Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/

1"""Matrix equation solver routines"""

2# Author: Jeffrey Armstrong <jeff@approximatrix.com>

3# February 24, 2012

5# Modified: Chad Fulton <ChadFulton@gmail.com>

6# June 19, 2014

8# Modified: Ilhan Polat <ilhanpolat@gmail.com>

9# September 13, 2016

11import warnings

12import numpy as np

13from numpy.linalg import inv, LinAlgError, norm, cond, svd

15from ._basic import solve, solve_triangular, matrix_balance

16from .lapack import get_lapack_funcs

17from ._decomp_schur import schur

18from ._decomp_lu import lu

19from ._decomp_qr import qr

20from ._decomp_qz import ordqz

21from ._decomp import _asarray_validated

22from ._special_matrices import kron, block_diag

24__all__ = ['solve_sylvester',

25 'solve_continuous_lyapunov', 'solve_discrete_lyapunov',

26 'solve_lyapunov',

27 'solve_continuous_are', 'solve_discrete_are']

30def solve_sylvester(a, b, q):

31 """

32 Computes a solution (X) to the Sylvester equation :math:`AX + XB = Q`.

34 Parameters

35 ----------

36 a : (M, M) array_like

37 Leading matrix of the Sylvester equation

38 b : (N, N) array_like

39 Trailing matrix of the Sylvester equation

40 q : (M, N) array_like

41 Right-hand side

43 Returns

44 -------

45 x : (M, N) ndarray

46 The solution to the Sylvester equation.

48 Raises

49 ------

50 LinAlgError

51 If solution was not found

53 Notes

54 -----

55 Computes a solution to the Sylvester matrix equation via the Bartels-

56 Stewart algorithm. The A and B matrices first undergo Schur

57 decompositions. The resulting matrices are used to construct an

58 alternative Sylvester equation (``RY + YS^T = F``) where the R and S

59 matrices are in quasi-triangular form (or, when R, S or F are complex,

60 triangular form). The simplified equation is then solved using

61 ``*TRSYL`` from LAPACK directly.

63 .. versionadded:: 0.11.0

65 Examples

66 --------

67 Given `a`, `b`, and `q` solve for `x`:

69 >>> import numpy as np

70 >>> from scipy import linalg

71 >>> a = np.array([[-3, -2, 0], [-1, -1, 3], [3, -5, -1]])

72 >>> b = np.array([[1]])

73 >>> q = np.array([[1],[2],[3]])

74 >>> x = linalg.solve_sylvester(a, b, q)

75 >>> x

76 array([[ 0.0625],

77 [-0.5625],

78 [ 0.6875]])

79 >>> np.allclose(a.dot(x) + x.dot(b), q)

80 True

82 """

84 # Compute the Schur decomposition form of a

85 r, u = schur(a, output='real')

87 # Compute the Schur decomposition of b

88 s, v = schur(b.conj().transpose(), output='real')

90 # Construct f = u'*q*v

91 f = np.dot(np.dot(u.conj().transpose(), q), v)

93 # Call the Sylvester equation solver

94 trsyl, = get_lapack_funcs(('trsyl',), (r, s, f))

95 if trsyl is None:

96 raise RuntimeError('LAPACK implementation does not contain a proper '

97 'Sylvester equation solver (TRSYL)')

98 y, scale, info = trsyl(r, s, f, tranb='C')

100 y = scale*y

101

102 if info < 0:

103 raise LinAlgError("Illegal value encountered in "

104 "the %d term" % (-info,))

105

106 return np.dot(np.dot(u, y), v.conj().transpose())

107

108

109def solve_continuous_lyapunov(a, q):

110 """

111 Solves the continuous Lyapunov equation :math:`AX + XA^H = Q`.

112

113 Uses the Bartels-Stewart algorithm to find :math:`X`.

114

115 Parameters

116 ----------

117 a : array_like

118 A square matrix

119

120 q : array_like

121 Right-hand side square matrix

122

123 Returns

124 -------

125 x : ndarray

126 Solution to the continuous Lyapunov equation

127

128 See Also

129 --------

130 solve_discrete_lyapunov : computes the solution to the discrete-time

131 Lyapunov equation

132 solve_sylvester : computes the solution to the Sylvester equation

133

134 Notes

135 -----

136 The continuous Lyapunov equation is a special form of the Sylvester

137 equation, hence this solver relies on LAPACK routine ?TRSYL.

138

139 .. versionadded:: 0.11.0

140

141 Examples

142 --------

143 Given `a` and `q` solve for `x`:

144

145 >>> import numpy as np

146 >>> from scipy import linalg

147 >>> a = np.array([[-3, -2, 0], [-1, -1, 0], [0, -5, -1]])

148 >>> b = np.array([2, 4, -1])

149 >>> q = np.eye(3)

150 >>> x = linalg.solve_continuous_lyapunov(a, q)

151 >>> x

152 array([[ -0.75 , 0.875 , -3.75 ],

153 [ 0.875 , -1.375 , 5.3125],

154 [ -3.75 , 5.3125, -27.0625]])

155 >>> np.allclose(a.dot(x) + x.dot(a.T), q)

156 True

157 """

158

159 a = np.atleast_2d(_asarray_validated(a, check_finite=True))

160 q = np.atleast_2d(_asarray_validated(q, check_finite=True))

161

162 r_or_c = float

163

164 for ind, _ in enumerate((a, q)):

165 if np.iscomplexobj(_):

166 r_or_c = complex

167

168 if not np.equal(*_.shape):

169 raise ValueError("Matrix {} should be square.".format("aq"[ind]))

170

171 # Shape consistency check

172 if a.shape != q.shape:

173 raise ValueError("Matrix a and q should have the same shape.")

174

175 # Compute the Schur decomposition form of a

176 r, u = schur(a, output='real')

177

178 # Construct f = u'*q*u

179 f = u.conj().T.dot(q.dot(u))

180

181 # Call the Sylvester equation solver

182 trsyl = get_lapack_funcs('trsyl', (r, f))

183

184 dtype_string = 'T' if r_or_c == float else 'C'

185 y, scale, info = trsyl(r, r, f, tranb=dtype_string)

186

187 if info < 0:

188 raise ValueError('?TRSYL exited with the internal error '

189 '"illegal value in argument number {}.". See '

190 'LAPACK documentation for the ?TRSYL error codes.'

191 ''.format(-info))

192 elif info == 1:

193 warnings.warn('Input "a" has an eigenvalue pair whose sum is '

194 'very close to or exactly zero. The solution is '

195 'obtained via perturbing the coefficients.',

196 RuntimeWarning)

197 y *= scale

198

199 return u.dot(y).dot(u.conj().T)

200

201

202# For backwards compatibility, keep the old name

203solve_lyapunov = solve_continuous_lyapunov

204

205

206def _solve_discrete_lyapunov_direct(a, q):

207 """

208 Solves the discrete Lyapunov equation directly.

209

210 This function is called by the `solve_discrete_lyapunov` function with

211 `method=direct`. It is not supposed to be called directly.

212 """

213

214 lhs = kron(a, a.conj())

215 lhs = np.eye(lhs.shape[0]) - lhs

216 x = solve(lhs, q.flatten())

217

218 return np.reshape(x, q.shape)

219

220

221def _solve_discrete_lyapunov_bilinear(a, q):

222 """

223 Solves the discrete Lyapunov equation using a bilinear transformation.

224

225 This function is called by the `solve_discrete_lyapunov` function with

226 `method=bilinear`. It is not supposed to be called directly.

227 """

228 eye = np.eye(a.shape[0])

229 aH = a.conj().transpose()

230 aHI_inv = inv(aH + eye)

231 b = np.dot(aH - eye, aHI_inv)

232 c = 2*np.dot(np.dot(inv(a + eye), q), aHI_inv)

233 return solve_lyapunov(b.conj().transpose(), -c)

234

235

236def solve_discrete_lyapunov(a, q, method=None):

237 """

238 Solves the discrete Lyapunov equation :math:`AXA^H - X + Q = 0`.

239

240 Parameters

241 ----------

242 a, q : (M, M) array_like

243 Square matrices corresponding to A and Q in the equation

244 above respectively. Must have the same shape.

245

246 method : {'direct', 'bilinear'}, optional

247 Type of solver.

248

249 If not given, chosen to be ``direct`` if ``M`` is less than 10 and

250 ``bilinear`` otherwise.

251

252 Returns

253 -------

254 x : ndarray

255 Solution to the discrete Lyapunov equation

256

257 See Also

258 --------

259 solve_continuous_lyapunov : computes the solution to the continuous-time

260 Lyapunov equation

261

262 Notes

263 -----

264 This section describes the available solvers that can be selected by the

265 'method' parameter. The default method is *direct* if ``M`` is less than 10

266 and ``bilinear`` otherwise.

267

268 Method *direct* uses a direct analytical solution to the discrete Lyapunov

269 equation. The algorithm is given in, for example, [1]_. However, it requires

270 the linear solution of a system with dimension :math:`M^2` so that

271 performance degrades rapidly for even moderately sized matrices.

272

273 Method *bilinear* uses a bilinear transformation to convert the discrete

274 Lyapunov equation to a continuous Lyapunov equation :math:`(BX+XB'=-C)`

275 where :math:`B=(A-I)(A+I)^{-1}` and

276 :math:`C=2(A' + I)^{-1} Q (A + I)^{-1}`. The continuous equation can be

277 efficiently solved since it is a special case of a Sylvester equation.

278 The transformation algorithm is from Popov (1964) as described in [2]_.

279

280 .. versionadded:: 0.11.0

281

282 References

283 ----------

284 .. [1] Hamilton, James D. Time Series Analysis, Princeton: Princeton

285 University Press, 1994. 265. Print.

286 http://doc1.lbfl.li/aca/FLMF037168.pdf

287 .. [2] Gajic, Z., and M.T.J. Qureshi. 2008.

288 Lyapunov Matrix Equation in System Stability and Control.

289 Dover Books on Engineering Series. Dover Publications.

290

291 Examples

292 --------

293 Given `a` and `q` solve for `x`:

294

295 >>> import numpy as np

296 >>> from scipy import linalg

297 >>> a = np.array([[0.2, 0.5],[0.7, -0.9]])

298 >>> q = np.eye(2)

299 >>> x = linalg.solve_discrete_lyapunov(a, q)

300 >>> x

301 array([[ 0.70872893, 1.43518822],

302 [ 1.43518822, -2.4266315 ]])

303 >>> np.allclose(a.dot(x).dot(a.T)-x, -q)

304 True

305

306 """

307 a = np.asarray(a)

308 q = np.asarray(q)

309 if method is None:

310 # Select automatically based on size of matrices

311 if a.shape[0] >= 10:

312 method = 'bilinear'

313 else:

314 method = 'direct'

315

316 meth = method.lower()

317

318 if meth == 'direct':

319 x = _solve_discrete_lyapunov_direct(a, q)

320 elif meth == 'bilinear':

321 x = _solve_discrete_lyapunov_bilinear(a, q)

322 else:

323 raise ValueError('Unknown solver %s' % method)

324

325 return x

326

327

328def solve_continuous_are(a, b, q, r, e=None, s=None, balanced=True):

329 r"""

330 Solves the continuous-time algebraic Riccati equation (CARE).

331

332 The CARE is defined as

333

334 .. math::

335

336 X A + A^H X - X B R^{-1} B^H X + Q = 0

337

338 The limitations for a solution to exist are :

339

340 * All eigenvalues of :math:`A` on the right half plane, should be

341 controllable.

342

343 * The associated hamiltonian pencil (See Notes), should have

344 eigenvalues sufficiently away from the imaginary axis.

345

346 Moreover, if ``e`` or ``s`` is not precisely ``None``, then the

347 generalized version of CARE

348

349 .. math::

350

351 E^HXA + A^HXE - (E^HXB + S) R^{-1} (B^HXE + S^H) + Q = 0

352

353 is solved. When omitted, ``e`` is assumed to be the identity and ``s``

354 is assumed to be the zero matrix with sizes compatible with ``a`` and

355 ``b``, respectively.

356

357 Parameters

358 ----------

359 a : (M, M) array_like

360 Square matrix

361 b : (M, N) array_like

362 Input

363 q : (M, M) array_like

364 Input

365 r : (N, N) array_like

366 Nonsingular square matrix

367 e : (M, M) array_like, optional

368 Nonsingular square matrix

369 s : (M, N) array_like, optional

370 Input

371 balanced : bool, optional

372 The boolean that indicates whether a balancing step is performed

373 on the data. The default is set to True.

374

375 Returns

376 -------

377 x : (M, M) ndarray

378 Solution to the continuous-time algebraic Riccati equation.

379

380 Raises

381 ------

382 LinAlgError

383 For cases where the stable subspace of the pencil could not be

384 isolated. See Notes section and the references for details.

385

386 See Also

387 --------

388 solve_discrete_are : Solves the discrete-time algebraic Riccati equation

389

390 Notes

391 -----

392 The equation is solved by forming the extended hamiltonian matrix pencil,

393 as described in [1]_, :math:`H - \lambda J` given by the block matrices ::

394

395 [ A 0 B ] [ E 0 0 ]

396 [-Q -A^H -S ] - \lambda * [ 0 E^H 0 ]

397 [ S^H B^H R ] [ 0 0 0 ]

398

399 and using a QZ decomposition method.

400

401 In this algorithm, the fail conditions are linked to the symmetry

402 of the product :math:`U_2 U_1^{-1}` and condition number of

403 :math:`U_1`. Here, :math:`U` is the 2m-by-m matrix that holds the

404 eigenvectors spanning the stable subspace with 2-m rows and partitioned

405 into two m-row matrices. See [1]_ and [2]_ for more details.

406

407 In order to improve the QZ decomposition accuracy, the pencil goes

408 through a balancing step where the sum of absolute values of

409 :math:`H` and :math:`J` entries (after removing the diagonal entries of

410 the sum) is balanced following the recipe given in [3]_.

411

412 .. versionadded:: 0.11.0

413

414 References

415 ----------

416 .. [1] P. van Dooren , "A Generalized Eigenvalue Approach For Solving

417 Riccati Equations.", SIAM Journal on Scientific and Statistical

418 Computing, Vol.2(2), :doi:`10.1137/0902010`

419

420 .. [2] A.J. Laub, "A Schur Method for Solving Algebraic Riccati

421 Equations.", Massachusetts Institute of Technology. Laboratory for

422 Information and Decision Systems. LIDS-R ; 859. Available online :

423 http://hdl.handle.net/1721.1/1301

424

425 .. [3] P. Benner, "Symplectic Balancing of Hamiltonian Matrices", 2001,

426 SIAM J. Sci. Comput., 2001, Vol.22(5), :doi:`10.1137/S1064827500367993`

427

428 Examples

429 --------

430 Given `a`, `b`, `q`, and `r` solve for `x`:

431

432 >>> import numpy as np

433 >>> from scipy import linalg

434 >>> a = np.array([[4, 3], [-4.5, -3.5]])

435 >>> b = np.array([[1], [-1]])

436 >>> q = np.array([[9, 6], [6, 4.]])

437 >>> r = 1

438 >>> x = linalg.solve_continuous_are(a, b, q, r)

439 >>> x

440 array([[ 21.72792206, 14.48528137],

441 [ 14.48528137, 9.65685425]])

442 >>> np.allclose(a.T.dot(x) + x.dot(a)-x.dot(b).dot(b.T).dot(x), -q)

443 True

444

445 """

446

447 # Validate input arguments

448 a, b, q, r, e, s, m, n, r_or_c, gen_are = _are_validate_args(

449 a, b, q, r, e, s, 'care')

450

451 H = np.empty((2*m+n, 2*m+n), dtype=r_or_c)

452 H[:m, :m] = a

453 H[:m, m:2*m] = 0.

454 H[:m, 2*m:] = b

455 H[m:2*m, :m] = -q

456 H[m:2*m, m:2*m] = -a.conj().T

457 H[m:2*m, 2*m:] = 0. if s is None else -s

458 H[2*m:, :m] = 0. if s is None else s.conj().T

459 H[2*m:, m:2*m] = b.conj().T

460 H[2*m:, 2*m:] = r

461

462 if gen_are and e is not None:

463 J = block_diag(e, e.conj().T, np.zeros_like(r, dtype=r_or_c))

464 else:

465 J = block_diag(np.eye(2*m), np.zeros_like(r, dtype=r_or_c))

466

467 if balanced:

468 # xGEBAL does not remove the diagonals before scaling. Also

469 # to avoid destroying the Symplectic structure, we follow Ref.3

470 M = np.abs(H) + np.abs(J)

471 M[np.diag_indices_from(M)] = 0.

472 _, (sca, _) = matrix_balance(M, separate=1, permute=0)

473 # do we need to bother?

474 if not np.allclose(sca, np.ones_like(sca)):

475 # Now impose diag(D,inv(D)) from Benner where D is

476 # square root of s_i/s_(n+i) for i=0,....

477 sca = np.log2(sca)

478 # NOTE: Py3 uses "Bankers Rounding: round to the nearest even" !!

479 s = np.round((sca[m:2*m] - sca[:m])/2)

480 sca = 2 ** np.r_[s, -s, sca[2*m:]]

481 # Elementwise multiplication via broadcasting.

482 elwisescale = sca[:, None] * np.reciprocal(sca)

483 H *= elwisescale

484 J *= elwisescale

485

486 # Deflate the pencil to 2m x 2m ala Ref.1, eq.(55)

487 q, r = qr(H[:, -n:])

488 H = q[:, n:].conj().T.dot(H[:, :2*m])

489 J = q[:2*m, n:].conj().T.dot(J[:2*m, :2*m])

490

491 # Decide on which output type is needed for QZ

492 out_str = 'real' if r_or_c == float else 'complex'

493

494 _, _, _, _, _, u = ordqz(H, J, sort='lhp', overwrite_a=True,

495 overwrite_b=True, check_finite=False,

496 output=out_str)

497

498 # Get the relevant parts of the stable subspace basis

499 if e is not None:

500 u, _ = qr(np.vstack((e.dot(u[:m, :m]), u[m:, :m])))

501 u00 = u[:m, :m]

502 u10 = u[m:, :m]

503

504 # Solve via back-substituion after checking the condition of u00

505 up, ul, uu = lu(u00)

506 if 1/cond(uu) < np.spacing(1.):

507 raise LinAlgError('Failed to find a finite solution.')

508

509 # Exploit the triangular structure

510 x = solve_triangular(ul.conj().T,

511 solve_triangular(uu.conj().T,

512 u10.conj().T,

513 lower=True),

514 unit_diagonal=True,

515 ).conj().T.dot(up.conj().T)

516 if balanced:

517 x *= sca[:m, None] * sca[:m]

518

519 # Check the deviation from symmetry for lack of success

520 # See proof of Thm.5 item 3 in [2]

521 u_sym = u00.conj().T.dot(u10)

522 n_u_sym = norm(u_sym, 1)

523 u_sym = u_sym - u_sym.conj().T

524 sym_threshold = np.max([np.spacing(1000.), 0.1*n_u_sym])

525

526 if norm(u_sym, 1) > sym_threshold:

527 raise LinAlgError('The associated Hamiltonian pencil has eigenvalues '

528 'too close to the imaginary axis')

529

530 return (x + x.conj().T)/2

531

532

533def solve_discrete_are(a, b, q, r, e=None, s=None, balanced=True):

534 r"""

535 Solves the discrete-time algebraic Riccati equation (DARE).

536

537 The DARE is defined as

538

539 .. math::

540

541 A^HXA - X - (A^HXB) (R + B^HXB)^{-1} (B^HXA) + Q = 0

542

543 The limitations for a solution to exist are :

544

545 * All eigenvalues of :math:`A` outside the unit disc, should be

546 controllable.

547

548 * The associated symplectic pencil (See Notes), should have

549 eigenvalues sufficiently away from the unit circle.

550

551 Moreover, if ``e`` and ``s`` are not both precisely ``None``, then the

552 generalized version of DARE

553

554 .. math::

555

556 A^HXA - E^HXE - (A^HXB+S) (R+B^HXB)^{-1} (B^HXA+S^H) + Q = 0

557

558 is solved. When omitted, ``e`` is assumed to be the identity and ``s``

559 is assumed to be the zero matrix.

560

561 Parameters

562 ----------

563 a : (M, M) array_like

564 Square matrix

565 b : (M, N) array_like

566 Input

567 q : (M, M) array_like

568 Input

569 r : (N, N) array_like

570 Square matrix

571 e : (M, M) array_like, optional

572 Nonsingular square matrix

573 s : (M, N) array_like, optional

574 Input

575 balanced : bool

576 The boolean that indicates whether a balancing step is performed

577 on the data. The default is set to True.

578

579 Returns

580 -------

581 x : (M, M) ndarray

582 Solution to the discrete algebraic Riccati equation.

583

584 Raises

585 ------

586 LinAlgError

587 For cases where the stable subspace of the pencil could not be

588 isolated. See Notes section and the references for details.

589

590 See Also

591 --------

592 solve_continuous_are : Solves the continuous algebraic Riccati equation

593

594 Notes

595 -----

596 The equation is solved by forming the extended symplectic matrix pencil,

597 as described in [1]_, :math:`H - \lambda J` given by the block matrices ::

598

599 [ A 0 B ] [ E 0 B ]

600 [ -Q E^H -S ] - \lambda * [ 0 A^H 0 ]

601 [ S^H 0 R ] [ 0 -B^H 0 ]

602

603 and using a QZ decomposition method.

604

605 In this algorithm, the fail conditions are linked to the symmetry

606 of the product :math:`U_2 U_1^{-1}` and condition number of

607 :math:`U_1`. Here, :math:`U` is the 2m-by-m matrix that holds the

608 eigenvectors spanning the stable subspace with 2-m rows and partitioned

609 into two m-row matrices. See [1]_ and [2]_ for more details.

610

611 In order to improve the QZ decomposition accuracy, the pencil goes

612 through a balancing step where the sum of absolute values of

613 :math:`H` and :math:`J` rows/cols (after removing the diagonal entries)

614 is balanced following the recipe given in [3]_. If the data has small

615 numerical noise, balancing may amplify their effects and some clean up

616 is required.

617

618 .. versionadded:: 0.11.0

619

620 References

621 ----------

622 .. [1] P. van Dooren , "A Generalized Eigenvalue Approach For Solving

623 Riccati Equations.", SIAM Journal on Scientific and Statistical

624 Computing, Vol.2(2), :doi:`10.1137/0902010`

625

626 .. [2] A.J. Laub, "A Schur Method for Solving Algebraic Riccati

627 Equations.", Massachusetts Institute of Technology. Laboratory for

628 Information and Decision Systems. LIDS-R ; 859. Available online :

629 http://hdl.handle.net/1721.1/1301

630

631 .. [3] P. Benner, "Symplectic Balancing of Hamiltonian Matrices", 2001,

632 SIAM J. Sci. Comput., 2001, Vol.22(5), :doi:`10.1137/S1064827500367993`

633

634 Examples

635 --------

636 Given `a`, `b`, `q`, and `r` solve for `x`:

637

638 >>> import numpy as np

639 >>> from scipy import linalg as la

640 >>> a = np.array([[0, 1], [0, -1]])

641 >>> b = np.array([[1, 0], [2, 1]])

642 >>> q = np.array([[-4, -4], [-4, 7]])

643 >>> r = np.array([[9, 3], [3, 1]])

644 >>> x = la.solve_discrete_are(a, b, q, r)

645 >>> x

646 array([[-4., -4.],

647 [-4., 7.]])

648 >>> R = la.solve(r + b.T.dot(x).dot(b), b.T.dot(x).dot(a))

649 >>> np.allclose(a.T.dot(x).dot(a) - x - a.T.dot(x).dot(b).dot(R), -q)

650 True

651

652 """

653

654 # Validate input arguments

655 a, b, q, r, e, s, m, n, r_or_c, gen_are = _are_validate_args(

656 a, b, q, r, e, s, 'dare')

657

658 # Form the matrix pencil

659 H = np.zeros((2*m+n, 2*m+n), dtype=r_or_c)

660 H[:m, :m] = a

661 H[:m, 2*m:] = b

662 H[m:2*m, :m] = -q

663 H[m:2*m, m:2*m] = np.eye(m) if e is None else e.conj().T

664 H[m:2*m, 2*m:] = 0. if s is None else -s

665 H[2*m:, :m] = 0. if s is None else s.conj().T

666 H[2*m:, 2*m:] = r

667

668 J = np.zeros_like(H, dtype=r_or_c)

669 J[:m, :m] = np.eye(m) if e is None else e

670 J[m:2*m, m:2*m] = a.conj().T

671 J[2*m:, m:2*m] = -b.conj().T

672

673 if balanced:

674 # xGEBAL does not remove the diagonals before scaling. Also

675 # to avoid destroying the Symplectic structure, we follow Ref.3

676 M = np.abs(H) + np.abs(J)

677 M[np.diag_indices_from(M)] = 0.

678 _, (sca, _) = matrix_balance(M, separate=1, permute=0)

679 # do we need to bother?

680 if not np.allclose(sca, np.ones_like(sca)):

681 # Now impose diag(D,inv(D)) from Benner where D is

682 # square root of s_i/s_(n+i) for i=0,....

683 sca = np.log2(sca)

684 # NOTE: Py3 uses "Bankers Rounding: round to the nearest even" !!

685 s = np.round((sca[m:2*m] - sca[:m])/2)

686 sca = 2 ** np.r_[s, -s, sca[2*m:]]

687 # Elementwise multiplication via broadcasting.

688 elwisescale = sca[:, None] * np.reciprocal(sca)

689 H *= elwisescale

690 J *= elwisescale

691

692 # Deflate the pencil by the R column ala Ref.1

693 q_of_qr, _ = qr(H[:, -n:])

694 H = q_of_qr[:, n:].conj().T.dot(H[:, :2*m])

695 J = q_of_qr[:, n:].conj().T.dot(J[:, :2*m])

696

697 # Decide on which output type is needed for QZ

698 out_str = 'real' if r_or_c == float else 'complex'

699

700 _, _, _, _, _, u = ordqz(H, J, sort='iuc',

701 overwrite_a=True,

702 overwrite_b=True,

703 check_finite=False,

704 output=out_str)

705

706 # Get the relevant parts of the stable subspace basis

707 if e is not None:

708 u, _ = qr(np.vstack((e.dot(u[:m, :m]), u[m:, :m])))

709 u00 = u[:m, :m]

710 u10 = u[m:, :m]

711

712 # Solve via back-substituion after checking the condition of u00

713 up, ul, uu = lu(u00)

714

715 if 1/cond(uu) < np.spacing(1.):

716 raise LinAlgError('Failed to find a finite solution.')

717

718 # Exploit the triangular structure

719 x = solve_triangular(ul.conj().T,

720 solve_triangular(uu.conj().T,

721 u10.conj().T,

722 lower=True),

723 unit_diagonal=True,

724 ).conj().T.dot(up.conj().T)

725 if balanced:

726 x *= sca[:m, None] * sca[:m]

727

728 # Check the deviation from symmetry for lack of success

729 # See proof of Thm.5 item 3 in [2]

730 u_sym = u00.conj().T.dot(u10)

731 n_u_sym = norm(u_sym, 1)

732 u_sym = u_sym - u_sym.conj().T

733 sym_threshold = np.max([np.spacing(1000.), 0.1*n_u_sym])

734

735 if norm(u_sym, 1) > sym_threshold:

736 raise LinAlgError('The associated symplectic pencil has eigenvalues'

737 'too close to the unit circle')

738

739 return (x + x.conj().T)/2

740

741

742def _are_validate_args(a, b, q, r, e, s, eq_type='care'):

743 """

744 A helper function to validate the arguments supplied to the

745 Riccati equation solvers. Any discrepancy found in the input

746 matrices leads to a ``ValueError`` exception.

747

748 Essentially, it performs:

749

750 - a check whether the input is free of NaN and Infs

751 - a pass for the data through ``numpy.atleast_2d()``

752 - squareness check of the relevant arrays

753 - shape consistency check of the arrays

754 - singularity check of the relevant arrays

755 - symmetricity check of the relevant matrices

756 - a check whether the regular or the generalized version is asked.

757

758 This function is used by ``solve_continuous_are`` and

759 ``solve_discrete_are``.

760

761 Parameters

762 ----------

763 a, b, q, r, e, s : array_like

764 Input data

765 eq_type : str

766 Accepted arguments are 'care' and 'dare'.

767

768 Returns

769 -------

770 a, b, q, r, e, s : ndarray

771 Regularized input data

772 m, n : int

773 shape of the problem

774 r_or_c : type

775 Data type of the problem, returns float or complex

776 gen_or_not : bool

777 Type of the equation, True for generalized and False for regular ARE.

778

779 """

780

781 if not eq_type.lower() in ('dare', 'care'):

782 raise ValueError("Equation type unknown. "

783 "Only 'care' and 'dare' is understood")

784

785 a = np.atleast_2d(_asarray_validated(a, check_finite=True))

786 b = np.atleast_2d(_asarray_validated(b, check_finite=True))

787 q = np.atleast_2d(_asarray_validated(q, check_finite=True))

788 r = np.atleast_2d(_asarray_validated(r, check_finite=True))

789

790 # Get the correct data types otherwise NumPy complains

791 # about pushing complex numbers into real arrays.

792 r_or_c = complex if np.iscomplexobj(b) else float

793

794 for ind, mat in enumerate((a, q, r)):

795 if np.iscomplexobj(mat):

796 r_or_c = complex

797

798 if not np.equal(*mat.shape):

799 raise ValueError("Matrix {} should be square.".format("aqr"[ind]))

800

801 # Shape consistency checks

802 m, n = b.shape

803 if m != a.shape[0]:

804 raise ValueError("Matrix a and b should have the same number of rows.")

805 if m != q.shape[0]:

806 raise ValueError("Matrix a and q should have the same shape.")

807 if n != r.shape[0]:

808 raise ValueError("Matrix b and r should have the same number of cols.")

809

810 # Check if the data matrices q, r are (sufficiently) hermitian

811 for ind, mat in enumerate((q, r)):

812 if norm(mat - mat.conj().T, 1) > np.spacing(norm(mat, 1))*100:

813 raise ValueError("Matrix {} should be symmetric/hermitian."

814 "".format("qr"[ind]))

815

816 # Continuous time ARE should have a nonsingular r matrix.

817 if eq_type == 'care':

818 min_sv = svd(r, compute_uv=False)[-1]

819 if min_sv == 0. or min_sv < np.spacing(1.)*norm(r, 1):

820 raise ValueError('Matrix r is numerically singular.')

821

822 # Check if the generalized case is required with omitted arguments

823 # perform late shape checking etc.

824 generalized_case = e is not None or s is not None

825

826 if generalized_case:

827 if e is not None:

828 e = np.atleast_2d(_asarray_validated(e, check_finite=True))

829 if not np.equal(*e.shape):

830 raise ValueError("Matrix e should be square.")

831 if m != e.shape[0]:

832 raise ValueError("Matrix a and e should have the same shape.")

833 # numpy.linalg.cond doesn't check for exact zeros and

834 # emits a runtime warning. Hence the following manual check.

835 min_sv = svd(e, compute_uv=False)[-1]

836 if min_sv == 0. or min_sv < np.spacing(1.) * norm(e, 1):

837 raise ValueError('Matrix e is numerically singular.')

838 if np.iscomplexobj(e):

839 r_or_c = complex

840 if s is not None:

841 s = np.atleast_2d(_asarray_validated(s, check_finite=True))

842 if s.shape != b.shape:

843 raise ValueError("Matrix b and s should have the same shape.")

844 if np.iscomplexobj(s):

845 r_or_c = complex

846

847 return a, b, q, r, e, s, m, n, r_or_c, generalized_case

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_solvers.py: 10%

215 statements