Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_qz.py: 13%
119 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import warnings
3import numpy as np
4from numpy import asarray_chkfinite
5from ._misc import LinAlgError, _datacopied, LinAlgWarning
6from .lapack import get_lapack_funcs
9__all__ = ['qz', 'ordqz']
11_double_precision = ['i', 'l', 'd']
14def _select_function(sort):
15 if callable(sort):
16 # assume the user knows what they're doing
17 sfunction = sort
18 elif sort == 'lhp':
19 sfunction = _lhp
20 elif sort == 'rhp':
21 sfunction = _rhp
22 elif sort == 'iuc':
23 sfunction = _iuc
24 elif sort == 'ouc':
25 sfunction = _ouc
26 else:
27 raise ValueError("sort parameter must be None, a callable, or "
28 "one of ('lhp','rhp','iuc','ouc')")
30 return sfunction
33def _lhp(x, y):
34 out = np.empty_like(x, dtype=bool)
35 nonzero = (y != 0)
36 # handles (x, y) = (0, 0) too
37 out[~nonzero] = False
38 out[nonzero] = (np.real(x[nonzero]/y[nonzero]) < 0.0)
39 return out
42def _rhp(x, y):
43 out = np.empty_like(x, dtype=bool)
44 nonzero = (y != 0)
45 # handles (x, y) = (0, 0) too
46 out[~nonzero] = False
47 out[nonzero] = (np.real(x[nonzero]/y[nonzero]) > 0.0)
48 return out
51def _iuc(x, y):
52 out = np.empty_like(x, dtype=bool)
53 nonzero = (y != 0)
54 # handles (x, y) = (0, 0) too
55 out[~nonzero] = False
56 out[nonzero] = (abs(x[nonzero]/y[nonzero]) < 1.0)
57 return out
60def _ouc(x, y):
61 out = np.empty_like(x, dtype=bool)
62 xzero = (x == 0)
63 yzero = (y == 0)
64 out[xzero & yzero] = False
65 out[~xzero & yzero] = True
66 out[~yzero] = (abs(x[~yzero]/y[~yzero]) > 1.0)
67 return out
70def _qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False,
71 overwrite_b=False, check_finite=True):
72 if sort is not None:
73 # Disabled due to segfaults on win32, see ticket 1717.
74 raise ValueError("The 'sort' input of qz() has to be None and will be "
75 "removed in a future release. Use ordqz instead.")
77 if output not in ['real', 'complex', 'r', 'c']:
78 raise ValueError("argument must be 'real', or 'complex'")
80 if check_finite:
81 a1 = asarray_chkfinite(A)
82 b1 = asarray_chkfinite(B)
83 else:
84 a1 = np.asarray(A)
85 b1 = np.asarray(B)
87 a_m, a_n = a1.shape
88 b_m, b_n = b1.shape
89 if not (a_m == a_n == b_m == b_n):
90 raise ValueError("Array dimensions must be square and agree")
92 typa = a1.dtype.char
93 if output in ['complex', 'c'] and typa not in ['F', 'D']:
94 if typa in _double_precision:
95 a1 = a1.astype('D')
96 typa = 'D'
97 else:
98 a1 = a1.astype('F')
99 typa = 'F'
100 typb = b1.dtype.char
101 if output in ['complex', 'c'] and typb not in ['F', 'D']:
102 if typb in _double_precision:
103 b1 = b1.astype('D')
104 typb = 'D'
105 else:
106 b1 = b1.astype('F')
107 typb = 'F'
109 overwrite_a = overwrite_a or (_datacopied(a1, A))
110 overwrite_b = overwrite_b or (_datacopied(b1, B))
112 gges, = get_lapack_funcs(('gges',), (a1, b1))
114 if lwork is None or lwork == -1:
115 # get optimal work array size
116 result = gges(lambda x: None, a1, b1, lwork=-1)
117 lwork = result[-2][0].real.astype(np.int_)
119 sfunction = lambda x: None
120 result = gges(sfunction, a1, b1, lwork=lwork, overwrite_a=overwrite_a,
121 overwrite_b=overwrite_b, sort_t=0)
123 info = result[-1]
124 if info < 0:
125 raise ValueError("Illegal value in argument {} of gges".format(-info))
126 elif info > 0 and info <= a_n:
127 warnings.warn("The QZ iteration failed. (a,b) are not in Schur "
128 "form, but ALPHAR(j), ALPHAI(j), and BETA(j) should be "
129 "correct for J={},...,N".format(info-1), LinAlgWarning,
130 stacklevel=3)
131 elif info == a_n+1:
132 raise LinAlgError("Something other than QZ iteration failed")
133 elif info == a_n+2:
134 raise LinAlgError("After reordering, roundoff changed values of some "
135 "complex eigenvalues so that leading eigenvalues "
136 "in the Generalized Schur form no longer satisfy "
137 "sort=True. This could also be due to scaling.")
138 elif info == a_n+3:
139 raise LinAlgError("Reordering failed in <s,d,c,z>tgsen")
141 return result, gges.typecode
144def qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False,
145 overwrite_b=False, check_finite=True):
146 """
147 QZ decomposition for generalized eigenvalues of a pair of matrices.
149 The QZ, or generalized Schur, decomposition for a pair of n-by-n
150 matrices (A,B) is::
152 (A,B) = (Q @ AA @ Z*, Q @ BB @ Z*)
154 where AA, BB is in generalized Schur form if BB is upper-triangular
155 with non-negative diagonal and AA is upper-triangular, or for real QZ
156 decomposition (``output='real'``) block upper triangular with 1x1
157 and 2x2 blocks. In this case, the 1x1 blocks correspond to real
158 generalized eigenvalues and 2x2 blocks are 'standardized' by making
159 the corresponding elements of BB have the form::
161 [ a 0 ]
162 [ 0 b ]
164 and the pair of corresponding 2x2 blocks in AA and BB will have a complex
165 conjugate pair of generalized eigenvalues. If (``output='complex'``) or
166 A and B are complex matrices, Z' denotes the conjugate-transpose of Z.
167 Q and Z are unitary matrices.
169 Parameters
170 ----------
171 A : (N, N) array_like
172 2-D array to decompose
173 B : (N, N) array_like
174 2-D array to decompose
175 output : {'real', 'complex'}, optional
176 Construct the real or complex QZ decomposition for real matrices.
177 Default is 'real'.
178 lwork : int, optional
179 Work array size. If None or -1, it is automatically computed.
180 sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
181 NOTE: THIS INPUT IS DISABLED FOR NOW. Use ordqz instead.
183 Specifies whether the upper eigenvalues should be sorted. A callable
184 may be passed that, given a eigenvalue, returns a boolean denoting
185 whether the eigenvalue should be sorted to the top-left (True). For
186 real matrix pairs, the sort function takes three real arguments
187 (alphar, alphai, beta). The eigenvalue
188 ``x = (alphar + alphai*1j)/beta``. For complex matrix pairs or
189 output='complex', the sort function takes two complex arguments
190 (alpha, beta). The eigenvalue ``x = (alpha/beta)``. Alternatively,
191 string parameters may be used:
193 - 'lhp' Left-hand plane (x.real < 0.0)
194 - 'rhp' Right-hand plane (x.real > 0.0)
195 - 'iuc' Inside the unit circle (x*x.conjugate() < 1.0)
196 - 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
198 Defaults to None (no sorting).
199 overwrite_a : bool, optional
200 Whether to overwrite data in a (may improve performance)
201 overwrite_b : bool, optional
202 Whether to overwrite data in b (may improve performance)
203 check_finite : bool, optional
204 If true checks the elements of `A` and `B` are finite numbers. If
205 false does no checking and passes matrix through to
206 underlying algorithm.
208 Returns
209 -------
210 AA : (N, N) ndarray
211 Generalized Schur form of A.
212 BB : (N, N) ndarray
213 Generalized Schur form of B.
214 Q : (N, N) ndarray
215 The left Schur vectors.
216 Z : (N, N) ndarray
217 The right Schur vectors.
219 See Also
220 --------
221 ordqz
223 Notes
224 -----
225 Q is transposed versus the equivalent function in Matlab.
227 .. versionadded:: 0.11.0
229 Examples
230 --------
231 >>> import numpy as np
232 >>> from scipy.linalg import qz
234 >>> A = np.array([[1, 2, -1], [5, 5, 5], [2, 4, -8]])
235 >>> B = np.array([[1, 1, -3], [3, 1, -1], [5, 6, -2]])
237 Compute the decomposition. The QZ decomposition is not unique, so
238 depending on the underlying library that is used, there may be
239 differences in the signs of coefficients in the following output.
241 >>> AA, BB, Q, Z = qz(A, B)
242 >>> AA
243 array([[-1.36949157, -4.05459025, 7.44389431],
244 [ 0. , 7.65653432, 5.13476017],
245 [ 0. , -0.65978437, 2.4186015 ]]) # may vary
246 >>> BB
247 array([[ 1.71890633, -1.64723705, -0.72696385],
248 [ 0. , 8.6965692 , -0. ],
249 [ 0. , 0. , 2.27446233]]) # may vary
250 >>> Q
251 array([[-0.37048362, 0.1903278 , 0.90912992],
252 [-0.90073232, 0.16534124, -0.40167593],
253 [ 0.22676676, 0.96769706, -0.11017818]]) # may vary
254 >>> Z
255 array([[-0.67660785, 0.63528924, -0.37230283],
256 [ 0.70243299, 0.70853819, -0.06753907],
257 [ 0.22088393, -0.30721526, -0.92565062]]) # may vary
259 Verify the QZ decomposition. With real output, we only need the
260 transpose of ``Z`` in the following expressions.
262 >>> Q @ AA @ Z.T # Should be A
263 array([[ 1., 2., -1.],
264 [ 5., 5., 5.],
265 [ 2., 4., -8.]])
266 >>> Q @ BB @ Z.T # Should be B
267 array([[ 1., 1., -3.],
268 [ 3., 1., -1.],
269 [ 5., 6., -2.]])
271 Repeat the decomposition, but with ``output='complex'``.
273 >>> AA, BB, Q, Z = qz(A, B, output='complex')
275 For conciseness in the output, we use ``np.set_printoptions()`` to set
276 the output precision of NumPy arrays to 3 and display tiny values as 0.
278 >>> np.set_printoptions(precision=3, suppress=True)
279 >>> AA
280 array([[-1.369+0.j , 2.248+4.237j, 4.861-5.022j],
281 [ 0. +0.j , 7.037+2.922j, 0.794+4.932j],
282 [ 0. +0.j , 0. +0.j , 2.655-1.103j]]) # may vary
283 >>> BB
284 array([[ 1.719+0.j , -1.115+1.j , -0.763-0.646j],
285 [ 0. +0.j , 7.24 +0.j , -3.144+3.322j],
286 [ 0. +0.j , 0. +0.j , 2.732+0.j ]]) # may vary
287 >>> Q
288 array([[ 0.326+0.175j, -0.273-0.029j, -0.886-0.052j],
289 [ 0.794+0.426j, -0.093+0.134j, 0.402-0.02j ],
290 [-0.2 -0.107j, -0.816+0.482j, 0.151-0.167j]]) # may vary
291 >>> Z
292 array([[ 0.596+0.32j , -0.31 +0.414j, 0.393-0.347j],
293 [-0.619-0.332j, -0.479+0.314j, 0.154-0.393j],
294 [-0.195-0.104j, 0.576+0.27j , 0.715+0.187j]]) # may vary
296 With complex arrays, we must use ``Z.conj().T`` in the following
297 expressions to verify the decomposition.
299 >>> Q @ AA @ Z.conj().T # Should be A
300 array([[ 1.-0.j, 2.-0.j, -1.-0.j],
301 [ 5.+0.j, 5.+0.j, 5.-0.j],
302 [ 2.+0.j, 4.+0.j, -8.+0.j]])
303 >>> Q @ BB @ Z.conj().T # Should be B
304 array([[ 1.+0.j, 1.+0.j, -3.+0.j],
305 [ 3.-0.j, 1.-0.j, -1.+0.j],
306 [ 5.+0.j, 6.+0.j, -2.+0.j]])
308 """
309 # output for real
310 # AA, BB, sdim, alphar, alphai, beta, vsl, vsr, work, info
311 # output for complex
312 # AA, BB, sdim, alpha, beta, vsl, vsr, work, info
313 result, _ = _qz(A, B, output=output, lwork=lwork, sort=sort,
314 overwrite_a=overwrite_a, overwrite_b=overwrite_b,
315 check_finite=check_finite)
316 return result[0], result[1], result[-4], result[-3]
319def ordqz(A, B, sort='lhp', output='real', overwrite_a=False,
320 overwrite_b=False, check_finite=True):
321 """QZ decomposition for a pair of matrices with reordering.
323 Parameters
324 ----------
325 A : (N, N) array_like
326 2-D array to decompose
327 B : (N, N) array_like
328 2-D array to decompose
329 sort : {callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
330 Specifies whether the upper eigenvalues should be sorted. A
331 callable may be passed that, given an ordered pair ``(alpha,
332 beta)`` representing the eigenvalue ``x = (alpha/beta)``,
333 returns a boolean denoting whether the eigenvalue should be
334 sorted to the top-left (True). For the real matrix pairs
335 ``beta`` is real while ``alpha`` can be complex, and for
336 complex matrix pairs both ``alpha`` and ``beta`` can be
337 complex. The callable must be able to accept a NumPy
338 array. Alternatively, string parameters may be used:
340 - 'lhp' Left-hand plane (x.real < 0.0)
341 - 'rhp' Right-hand plane (x.real > 0.0)
342 - 'iuc' Inside the unit circle (x*x.conjugate() < 1.0)
343 - 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
345 With the predefined sorting functions, an infinite eigenvalue
346 (i.e., ``alpha != 0`` and ``beta = 0``) is considered to lie in
347 neither the left-hand nor the right-hand plane, but it is
348 considered to lie outside the unit circle. For the eigenvalue
349 ``(alpha, beta) = (0, 0)``, the predefined sorting functions
350 all return `False`.
351 output : str {'real','complex'}, optional
352 Construct the real or complex QZ decomposition for real matrices.
353 Default is 'real'.
354 overwrite_a : bool, optional
355 If True, the contents of A are overwritten.
356 overwrite_b : bool, optional
357 If True, the contents of B are overwritten.
358 check_finite : bool, optional
359 If true checks the elements of `A` and `B` are finite numbers. If
360 false does no checking and passes matrix through to
361 underlying algorithm.
363 Returns
364 -------
365 AA : (N, N) ndarray
366 Generalized Schur form of A.
367 BB : (N, N) ndarray
368 Generalized Schur form of B.
369 alpha : (N,) ndarray
370 alpha = alphar + alphai * 1j. See notes.
371 beta : (N,) ndarray
372 See notes.
373 Q : (N, N) ndarray
374 The left Schur vectors.
375 Z : (N, N) ndarray
376 The right Schur vectors.
378 See Also
379 --------
380 qz
382 Notes
383 -----
384 On exit, ``(ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N``, will be the
385 generalized eigenvalues. ``ALPHAR(j) + ALPHAI(j)*i`` and
386 ``BETA(j),j=1,...,N`` are the diagonals of the complex Schur form (S,T)
387 that would result if the 2-by-2 diagonal blocks of the real generalized
388 Schur form of (A,B) were further reduced to triangular form using complex
389 unitary transformations. If ALPHAI(j) is zero, then the jth eigenvalue is
390 real; if positive, then the ``j``th and ``(j+1)``st eigenvalues are a
391 complex conjugate pair, with ``ALPHAI(j+1)`` negative.
393 .. versionadded:: 0.17.0
395 Examples
396 --------
397 >>> import numpy as np
398 >>> from scipy.linalg import ordqz
399 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
400 >>> B = np.array([[0, 6, 0, 0], [5, 0, 2, 1], [5, 2, 6, 6], [4, 7, 7, 7]])
401 >>> AA, BB, alpha, beta, Q, Z = ordqz(A, B, sort='lhp')
403 Since we have sorted for left half plane eigenvalues, negatives come first
405 >>> (alpha/beta).real < 0
406 array([ True, True, False, False], dtype=bool)
408 """
409 (AA, BB, _, *ab, Q, Z, _, _), typ = _qz(A, B, output=output, sort=None,
410 overwrite_a=overwrite_a,
411 overwrite_b=overwrite_b,
412 check_finite=check_finite)
414 if typ == 's':
415 alpha, beta = ab[0] + ab[1]*np.complex64(1j), ab[2]
416 elif typ == 'd':
417 alpha, beta = ab[0] + ab[1]*1.j, ab[2]
418 else:
419 alpha, beta = ab
421 sfunction = _select_function(sort)
422 select = sfunction(alpha, beta)
424 tgsen = get_lapack_funcs('tgsen', (AA, BB))
425 # the real case needs 4n + 16 lwork
426 lwork = 4*AA.shape[0] + 16 if typ in 'sd' else 1
427 AAA, BBB, *ab, QQ, ZZ, _, _, _, _, info = tgsen(select, AA, BB, Q, Z,
428 ijob=0,
429 lwork=lwork, liwork=1)
431 # Once more for tgsen output
432 if typ == 's':
433 alpha, beta = ab[0] + ab[1]*np.complex64(1j), ab[2]
434 elif typ == 'd':
435 alpha, beta = ab[0] + ab[1]*1.j, ab[2]
436 else:
437 alpha, beta = ab
439 if info < 0:
440 raise ValueError(f"Illegal value in argument {-info} of tgsen")
441 elif info == 1:
442 raise ValueError("Reordering of (A, B) failed because the transformed"
443 " matrix pair (A, B) would be too far from "
444 "generalized Schur form; the problem is very "
445 "ill-conditioned. (A, B) may have been partially "
446 "reordered.")
448 return AAA, BBB, alpha, beta, QQ, ZZ