Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-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 f'"illegal value in argument number {-info}.". See '

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

191 elif info == 1:

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

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

194 'obtained via perturbing the coefficients.',

195 RuntimeWarning, stacklevel=2)

196 y *= scale

197

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

199

200

201# For backwards compatibility, keep the old name

202solve_lyapunov = solve_continuous_lyapunov

203

204

205def _solve_discrete_lyapunov_direct(a, q):

206 """

207 Solves the discrete Lyapunov equation directly.

208

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

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

211 """

212

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

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

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

216

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

218

219

220def _solve_discrete_lyapunov_bilinear(a, q):

221 """

222 Solves the discrete Lyapunov equation using a bilinear transformation.

223

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

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

226 """

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

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

229 aHI_inv = inv(aH + eye)

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

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

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

233

234

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

236 """

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

238

239 Parameters

240 ----------

241 a, q : (M, M) array_like

242 Square matrices corresponding to A and Q in the equation

243 above respectively. Must have the same shape.

244

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

246 Type of solver.

247

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

249 ``bilinear`` otherwise.

250

251 Returns

252 -------

253 x : ndarray

254 Solution to the discrete Lyapunov equation

255

256 See Also

257 --------

258 solve_continuous_lyapunov : computes the solution to the continuous-time

259 Lyapunov equation

260

261 Notes

262 -----

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

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

265 and ``bilinear`` otherwise.

266

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

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

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

270 performance degrades rapidly for even moderately sized matrices.

271

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

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

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

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

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

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

278

279 .. versionadded:: 0.11.0

280

281 References

282 ----------

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

284 University Press, 1994. 265. Print.

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

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

287 Lyapunov Matrix Equation in System Stability and Control.

288 Dover Books on Engineering Series. Dover Publications.

289

290 Examples

291 --------

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

293

294 >>> import numpy as np

295 >>> from scipy import linalg

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

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

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

299 >>> x

300 array([[ 0.70872893, 1.43518822],

301 [ 1.43518822, -2.4266315 ]])

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

303 True

304

305 """

306 a = np.asarray(a)

307 q = np.asarray(q)

308 if method is None:

309 # Select automatically based on size of matrices

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

311 method = 'bilinear'

312 else:

313 method = 'direct'

314

315 meth = method.lower()

316

317 if meth == 'direct':

318 x = _solve_discrete_lyapunov_direct(a, q)

319 elif meth == 'bilinear':

320 x = _solve_discrete_lyapunov_bilinear(a, q)

321 else:

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

323

324 return x

325

326

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

328 r"""

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

330

331 The CARE is defined as

332

333 .. math::

334

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

336

337 The limitations for a solution to exist are :

338

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

340 controllable.

341

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

343 eigenvalues sufficiently away from the imaginary axis.

344

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

346 generalized version of CARE

347

348 .. math::

349

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

351

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

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

354 ``b``, respectively.

355

356 Parameters

357 ----------

358 a : (M, M) array_like

359 Square matrix

360 b : (M, N) array_like

361 Input

362 q : (M, M) array_like

363 Input

364 r : (N, N) array_like

365 Nonsingular square matrix

366 e : (M, M) array_like, optional

367 Nonsingular square matrix

368 s : (M, N) array_like, optional

369 Input

370 balanced : bool, optional

371 The boolean that indicates whether a balancing step is performed

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

373

374 Returns

375 -------

376 x : (M, M) ndarray

377 Solution to the continuous-time algebraic Riccati equation.

378

379 Raises

380 ------

381 LinAlgError

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

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

384

385 See Also

386 --------

387 solve_discrete_are : Solves the discrete-time algebraic Riccati equation

388

389 Notes

390 -----

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

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

393

394 [ A 0 B ] [ E 0 0 ]

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

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

397

398 and using a QZ decomposition method.

399

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

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

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

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

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

405

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

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

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

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

410

411 .. versionadded:: 0.11.0

412

413 References

414 ----------

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

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

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

418

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

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

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

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

423

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

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

426

427 Examples

428 --------

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

430

431 >>> import numpy as np

432 >>> from scipy import linalg

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

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

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

436 >>> r = 1

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

438 >>> x

439 array([[ 21.72792206, 14.48528137],

440 [ 14.48528137, 9.65685425]])

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

442 True

443

444 """

445

446 # Validate input arguments

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

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

449

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

451 H[:m, :m] = a

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

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

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

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

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

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

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

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

460

461 if gen_are and e is not None:

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

463 else:

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

465

466 if balanced:

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

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

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

470 np.fill_diagonal(M, 0.)

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

472 # do we need to bother?

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

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

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

476 sca = np.log2(sca)

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

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

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

480 # Elementwise multiplication via broadcasting.

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

482 H *= elwisescale

483 J *= elwisescale

484

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

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

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

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

489

490 # Decide on which output type is needed for QZ

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

492

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

494 overwrite_b=True, check_finite=False,

495 output=out_str)

496

497 # Get the relevant parts of the stable subspace basis

498 if e is not None:

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

500 u00 = u[:m, :m]

501 u10 = u[m:, :m]

502

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

504 up, ul, uu = lu(u00)

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

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

507

508 # Exploit the triangular structure

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

510 solve_triangular(uu.conj().T,

511 u10.conj().T,

512 lower=True),

513 unit_diagonal=True,

514 ).conj().T.dot(up.conj().T)

515 if balanced:

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

517

518 # Check the deviation from symmetry for lack of success

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

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

521 n_u_sym = norm(u_sym, 1)

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

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

524

525 if norm(u_sym, 1) > sym_threshold:

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

527 'too close to the imaginary axis')

528

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

530

531

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

533 r"""

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

535

536 The DARE is defined as

537

538 .. math::

539

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

541

542 The limitations for a solution to exist are :

543

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

545 controllable.

546

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

548 eigenvalues sufficiently away from the unit circle.

549

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

551 generalized version of DARE

552

553 .. math::

554

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

556

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

558 is assumed to be the zero matrix.

559

560 Parameters

561 ----------

562 a : (M, M) array_like

563 Square matrix

564 b : (M, N) array_like

565 Input

566 q : (M, M) array_like

567 Input

568 r : (N, N) array_like

569 Square matrix

570 e : (M, M) array_like, optional

571 Nonsingular square matrix

572 s : (M, N) array_like, optional

573 Input

574 balanced : bool

575 The boolean that indicates whether a balancing step is performed

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

577

578 Returns

579 -------

580 x : (M, M) ndarray

581 Solution to the discrete algebraic Riccati equation.

582

583 Raises

584 ------

585 LinAlgError

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

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

588

589 See Also

590 --------

591 solve_continuous_are : Solves the continuous algebraic Riccati equation

592

593 Notes

594 -----

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

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

597

598 [ A 0 B ] [ E 0 B ]

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

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

601

602 and using a QZ decomposition method.

603

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

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

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

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

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

609

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

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

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

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

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

615 is required.

616

617 .. versionadded:: 0.11.0

618

619 References

620 ----------

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

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

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

624

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

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

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

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

629

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

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

632

633 Examples

634 --------

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

636

637 >>> import numpy as np

638 >>> from scipy import linalg as la

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

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

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

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

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

644 >>> x

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

646 [-4., 7.]])

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

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

649 True

650

651 """

652

653 # Validate input arguments

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

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

656

657 # Form the matrix pencil

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

659 H[:m, :m] = a

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

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

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

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

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

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

666

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

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

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

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

671

672 if balanced:

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

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

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

676 np.fill_diagonal(M, 0.)

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

678 # do we need to bother?

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

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

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

682 sca = np.log2(sca)

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

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

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

686 # Elementwise multiplication via broadcasting.

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

688 H *= elwisescale

689 J *= elwisescale

690

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

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

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

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

695

696 # Decide on which output type is needed for QZ

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

698

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

700 overwrite_a=True,

701 overwrite_b=True,

702 check_finite=False,

703 output=out_str)

704

705 # Get the relevant parts of the stable subspace basis

706 if e is not None:

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

708 u00 = u[:m, :m]

709 u10 = u[m:, :m]

710

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

712 up, ul, uu = lu(u00)

713

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

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

716

717 # Exploit the triangular structure

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

719 solve_triangular(uu.conj().T,

720 u10.conj().T,

721 lower=True),

722 unit_diagonal=True,

723 ).conj().T.dot(up.conj().T)

724 if balanced:

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

726

727 # Check the deviation from symmetry for lack of success

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

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

730 n_u_sym = norm(u_sym, 1)

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

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

733

734 if norm(u_sym, 1) > sym_threshold:

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

736 'too close to the unit circle')

737

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

739

740

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

742 """

743 A helper function to validate the arguments supplied to the

744 Riccati equation solvers. Any discrepancy found in the input

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

746

747 Essentially, it performs:

748

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

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

751 - squareness check of the relevant arrays

752 - shape consistency check of the arrays

753 - singularity check of the relevant arrays

754 - symmetricity check of the relevant matrices

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

756

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

758 ``solve_discrete_are``.

759

760 Parameters

761 ----------

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

763 Input data

764 eq_type : str

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

766

767 Returns

768 -------

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

770 Regularized input data

771 m, n : int

772 shape of the problem

773 r_or_c : type

774 Data type of the problem, returns float or complex

775 gen_or_not : bool

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

777

778 """

779

780 if eq_type.lower() not in ("dare", "care"):

781 raise ValueError("Equation type unknown. "

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

783

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

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

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

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

788

789 # Get the correct data types otherwise NumPy complains

790 # about pushing complex numbers into real arrays.

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

792

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

794 if np.iscomplexobj(mat):

795 r_or_c = complex

796

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

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

799

800 # Shape consistency checks

801 m, n = b.shape

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

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

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

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

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

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

808

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

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

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

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

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

814

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

816 if eq_type == 'care':

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

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

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

820

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

822 # perform late shape checking etc.

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

824

825 if generalized_case:

826 if e is not None:

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

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

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

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

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

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

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

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

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

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

837 if np.iscomplexobj(e):

838 r_or_c = complex

839 if s is not None:

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

841 if s.shape != b.shape:

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

843 if np.iscomplexobj(s):

844 r_or_c = complex

845

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

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

215 statements