Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_solvers.py: 10%
215 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-14 06:37 +0000
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
102 if info < 0:
103 raise LinAlgError("Illegal value encountered in "
104 "the %d term" % (-info,))
106 return np.dot(np.dot(u, y), v.conj().transpose())
109def solve_continuous_lyapunov(a, q):
110 """
111 Solves the continuous Lyapunov equation :math:`AX + XA^H = Q`.
113 Uses the Bartels-Stewart algorithm to find :math:`X`.
115 Parameters
116 ----------
117 a : array_like
118 A square matrix
120 q : array_like
121 Right-hand side square matrix
123 Returns
124 -------
125 x : ndarray
126 Solution to the continuous Lyapunov equation
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
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.
139 .. versionadded:: 0.11.0
141 Examples
142 --------
143 Given `a` and `q` solve for `x`:
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 """
159 a = np.atleast_2d(_asarray_validated(a, check_finite=True))
160 q = np.atleast_2d(_asarray_validated(q, check_finite=True))
162 r_or_c = float
164 for ind, _ in enumerate((a, q)):
165 if np.iscomplexobj(_):
166 r_or_c = complex
168 if not np.equal(*_.shape):
169 raise ValueError("Matrix {} should be square.".format("aq"[ind]))
171 # Shape consistency check
172 if a.shape != q.shape:
173 raise ValueError("Matrix a and q should have the same shape.")
175 # Compute the Schur decomposition form of a
176 r, u = schur(a, output='real')
178 # Construct f = u'*q*u
179 f = u.conj().T.dot(q.dot(u))
181 # Call the Sylvester equation solver
182 trsyl = get_lapack_funcs('trsyl', (r, f))
184 dtype_string = 'T' if r_or_c == float else 'C'
185 y, scale, info = trsyl(r, r, f, tranb=dtype_string)
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
198 return u.dot(y).dot(u.conj().T)
201# For backwards compatibility, keep the old name
202solve_lyapunov = solve_continuous_lyapunov
205def _solve_discrete_lyapunov_direct(a, q):
206 """
207 Solves the discrete Lyapunov equation directly.
209 This function is called by the `solve_discrete_lyapunov` function with
210 `method=direct`. It is not supposed to be called directly.
211 """
213 lhs = kron(a, a.conj())
214 lhs = np.eye(lhs.shape[0]) - lhs
215 x = solve(lhs, q.flatten())
217 return np.reshape(x, q.shape)
220def _solve_discrete_lyapunov_bilinear(a, q):
221 """
222 Solves the discrete Lyapunov equation using a bilinear transformation.
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)
235def solve_discrete_lyapunov(a, q, method=None):
236 """
237 Solves the discrete Lyapunov equation :math:`AXA^H - X + Q = 0`.
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.
245 method : {'direct', 'bilinear'}, optional
246 Type of solver.
248 If not given, chosen to be ``direct`` if ``M`` is less than 10 and
249 ``bilinear`` otherwise.
251 Returns
252 -------
253 x : ndarray
254 Solution to the discrete Lyapunov equation
256 See Also
257 --------
258 solve_continuous_lyapunov : computes the solution to the continuous-time
259 Lyapunov equation
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.
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.
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]_.
279 .. versionadded:: 0.11.0
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.
290 Examples
291 --------
292 Given `a` and `q` solve for `x`:
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
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'
315 meth = method.lower()
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)
324 return x
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).
331 The CARE is defined as
333 .. math::
335 X A + A^H X - X B R^{-1} B^H X + Q = 0
337 The limitations for a solution to exist are :
339 * All eigenvalues of :math:`A` on the right half plane, should be
340 controllable.
342 * The associated hamiltonian pencil (See Notes), should have
343 eigenvalues sufficiently away from the imaginary axis.
345 Moreover, if ``e`` or ``s`` is not precisely ``None``, then the
346 generalized version of CARE
348 .. math::
350 E^HXA + A^HXE - (E^HXB + S) R^{-1} (B^HXE + S^H) + Q = 0
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.
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.
374 Returns
375 -------
376 x : (M, M) ndarray
377 Solution to the continuous-time algebraic Riccati equation.
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.
385 See Also
386 --------
387 solve_discrete_are : Solves the discrete-time algebraic Riccati equation
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 ::
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 ]
398 and using a QZ decomposition method.
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.
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]_.
411 .. versionadded:: 0.11.0
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`
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
424 .. [3] P. Benner, "Symplectic Balancing of Hamiltonian Matrices", 2001,
425 SIAM J. Sci. Comput., 2001, Vol.22(5), :doi:`10.1137/S1064827500367993`
427 Examples
428 --------
429 Given `a`, `b`, `q`, and `r` solve for `x`:
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
444 """
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')
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
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))
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
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])
490 # Decide on which output type is needed for QZ
491 out_str = 'real' if r_or_c == float else 'complex'
493 _, _, _, _, _, u = ordqz(H, J, sort='lhp', overwrite_a=True,
494 overwrite_b=True, check_finite=False,
495 output=out_str)
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]
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.')
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]
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])
525 if norm(u_sym, 1) > sym_threshold:
526 raise LinAlgError('The associated Hamiltonian pencil has eigenvalues '
527 'too close to the imaginary axis')
529 return (x + x.conj().T)/2
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).
536 The DARE is defined as
538 .. math::
540 A^HXA - X - (A^HXB) (R + B^HXB)^{-1} (B^HXA) + Q = 0
542 The limitations for a solution to exist are :
544 * All eigenvalues of :math:`A` outside the unit disc, should be
545 controllable.
547 * The associated symplectic pencil (See Notes), should have
548 eigenvalues sufficiently away from the unit circle.
550 Moreover, if ``e`` and ``s`` are not both precisely ``None``, then the
551 generalized version of DARE
553 .. math::
555 A^HXA - E^HXE - (A^HXB+S) (R+B^HXB)^{-1} (B^HXA+S^H) + Q = 0
557 is solved. When omitted, ``e`` is assumed to be the identity and ``s``
558 is assumed to be the zero matrix.
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.
578 Returns
579 -------
580 x : (M, M) ndarray
581 Solution to the discrete algebraic Riccati equation.
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.
589 See Also
590 --------
591 solve_continuous_are : Solves the continuous algebraic Riccati equation
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 ::
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 ]
602 and using a QZ decomposition method.
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.
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.
617 .. versionadded:: 0.11.0
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`
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
630 .. [3] P. Benner, "Symplectic Balancing of Hamiltonian Matrices", 2001,
631 SIAM J. Sci. Comput., 2001, Vol.22(5), :doi:`10.1137/S1064827500367993`
633 Examples
634 --------
635 Given `a`, `b`, `q`, and `r` solve for `x`:
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
651 """
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')
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
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
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
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])
696 # Decide on which output type is needed for QZ
697 out_str = 'real' if r_or_c == float else 'complex'
699 _, _, _, _, _, u = ordqz(H, J, sort='iuc',
700 overwrite_a=True,
701 overwrite_b=True,
702 check_finite=False,
703 output=out_str)
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]
711 # Solve via back-substituion after checking the condition of u00
712 up, ul, uu = lu(u00)
714 if 1/cond(uu) < np.spacing(1.):
715 raise LinAlgError('Failed to find a finite solution.')
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]
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])
734 if norm(u_sym, 1) > sym_threshold:
735 raise LinAlgError('The associated symplectic pencil has eigenvalues '
736 'too close to the unit circle')
738 return (x + x.conj().T)/2
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.
747 Essentially, it performs:
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.
757 This function is used by ``solve_continuous_are`` and
758 ``solve_discrete_are``.
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'.
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.
778 """
780 if eq_type.lower() not in ("dare", "care"):
781 raise ValueError("Equation type unknown. "
782 "Only 'care' and 'dare' is understood")
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))
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
793 for ind, mat in enumerate((a, q, r)):
794 if np.iscomplexobj(mat):
795 r_or_c = complex
797 if not np.equal(*mat.shape):
798 raise ValueError("Matrix {} should be square.".format("aqr"[ind]))
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.")
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]))
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.')
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
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
846 return a, b, q, r, e, s, m, n, r_or_c, generalized_case