Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 231 lines (197 sloc) 8.441 kb
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
1 import warnings
2
3 import numpy as np
4 from numpy import asarray_chkfinite, single
5
6 from misc import LinAlgError, _datacopied
7 from lapack import get_lapack_funcs
8
9 __all__ = ['qz']
10
11 _double_precision = ['i','l','d']
12
13 def _select_function(sort, typ):
14 if typ in ['F','D']:
15 if callable(sort):
16 #assume the user knows what they're doing
17 sfunction = sort
18 elif sort == 'lhp':
19 sfunction = lambda x,y: (np.real(x/y) < 0.0)
20 elif sort == 'rhp':
21 sfunction = lambda x,y: (np.real(x/y) >= 0.0)
22 elif sort == 'iuc':
23 sfunction = lambda x,y: (abs(x/y) <= 1.0)
24 elif sort == 'ouc':
25 sfunction = lambda x,y: (abs(x/y) > 1.0)
26 else:
27 raise ValueError("sort parameter must be None, a callable, or "
28 "one of ('lhp','rhp','iuc','ouc')")
29 elif typ in ['f','d']:
30 if callable(sort):
31 #assume the user knows what they're doing
32 sfunction = sort
33 elif sort == 'lhp':
34 sfunction = lambda x,y,z: (np.real((x+y*1j)/z) < 0.0)
35 elif sort == 'rhp':
36 sfunction = lambda x,y,z: (np.real((x+y*1j)/z) >= 0.0)
37 elif sort == 'iuc':
38 sfunction = lambda x,y,z: (abs((x+y*1j)/z) <= 1.0)
39 elif sort == 'ouc':
40 sfunction = lambda x,y,z: (abs((x+y*1j)/z) > 1.0)
41 else:
42 raise ValueError("sort parameter must be None, a callable, or "
43 "one of ('lhp','rhp','iuc','ouc')")
44 else: # to avoid an error later
45 raise ValueError("dtype %s not understood" % typ)
46 return sfunction
47
48
49 def qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False,
50 overwrite_b=False):
51 """
52 QZ decompostion for generalized eigenvalues of a pair of matrices.
53
54 The QZ, or generalized Schur, decomposition for a pair of N x N
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
55 nonsymmetric matrices (A,B) is::
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
56
57 (A,B) = (Q*AA*Z', Q*BB*Z')
58
59 where AA, BB is in generalized Schur form if BB is upper-triangular
60 with non-negative diagonal and AA is upper-triangular, or for real QZ
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
61 decomposition (``output='real'``) block upper triangular with 1x1
62 and 2x2 blocks. In this case, the 1x1 blocks correspond to real
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
63 generalized eigenvalues and 2x2 blocks are 'standardized' by making
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
64 the corresponding elements of BB have the form::
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
65
66 [ a 0 ]
67 [ 0 b ]
68
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
69 and the pair of corresponding 2x2 blocks in AA and BB will have a complex
70 conjugate pair of generalized eigenvalues. If (``output='complex'``) or
71 A and B are complex matrices, Z' denotes the conjugate-transpose of Z.
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
72 Q and Z are unitary matrices.
73
74 Parameters
75 ----------
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
76 A : array_like, shape (N,N)
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
77 2-D array to decompose.
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
78 B : array_like, shape (N,N)
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
79 2-D array to decompose.
80 output : {'real','complex'}, optional
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
81 Construct the real or complex QZ decomposition for real matrices.
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
82 Default is 'real'.
83 lwork : int, optional
84 Work array size. If None or -1, it is automatically computed.
85 sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
e1c61f1 @rgommers MAINT: disable sort keyword of linalg.qz, due to Windows issues. See…
rgommers authored
86 NOTE: THIS INPUT IS DISABLED FOR NOW, IT DOESN'T WORK WELL ON WINDOWS.
87
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
88 Specifies whether the upper eigenvalues should be sorted. A callable
89 may be passed that, given a eigenvalue, returns a boolean denoting
90 whether the eigenvalue should be sorted to the top-left (True). For
91 real matrix pairs, the sort function takes three real arguments
92 (alphar, alphai, beta). The eigenvalue x = (alphar + alphai*1j)/beta.
93 For complex matrix pairs or output='complex', the sort function
94 takes two complex arguments (alpha, beta). The eigenvalue
95 x = (alpha/beta).
96 Alternatively, string parameters may be used:
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
97
98 - 'lhp' Left-hand plane (x.real < 0.0)
99 - 'rhp' Right-hand plane (x.real > 0.0)
100 - 'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
101 - 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
102
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
103 Defaults to None (no sorting).
104
105 Returns
106 -------
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
107 AA : ndarray, shape (N,N)
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
108 Generalized Schur form of A.
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
109 BB : ndarray, shape (N,N)
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
110 Generalized Schur form of B.
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
111 Q : ndarray, shape (N,N)
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
112 The left Schur vectors.
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
113 Z : ndarray, shape (N,N)
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
114 The right Schur vectors.
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
115 sdim : int, optional
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
116 If sorting was requested, a fifth return value will contain the
117 number of eigenvalues for which the sort condition was True.
118
119 Notes
120 -----
121 Q is transposed versus the equivalent function in Matlab.
2ff9a95 @jseabold DOC: Add release related documentation
jseabold authored
122
123 .. versionadded:: 0.11.0
d829bb0 @rgommers DOC: merge more wiki doc edits.
rgommers authored
124
5ca6e31 @rgommers DOC: fix typos and formatting in linalg.qz docstring, and add an exam…
rgommers authored
125 Examples
126 --------
127 >>> from scipy import linalg
128 >>> np.random.seed(1234)
129 >>> A = np.arange(9).reshape((3, 3))
130 >>> B = np.random.randn(3, 3)
131
132 >>> AA, BB, Q, Z = linalg.qz(A, B)
133 >>> AA
134 array([[-13.40928183, -4.62471562, 1.09215523],
135 [ 0. , 0. , 1.22805978],
136 [ 0. , 0. , 0.31973817]])
137 >>> BB
138 array([[ 0.33362547, -1.37393632, 0.02179805],
139 [ 0. , 1.68144922, 0.74683866],
140 [ 0. , 0. , 0.9258294 ]])
141 >>> Q
142 array([[ 0.14134727, -0.97562773, 0.16784365],
143 [ 0.49835904, -0.07636948, -0.86360059],
144 [ 0.85537081, 0.20571399, 0.47541828]])
145 >>> Z
146 array([[-0.24900855, -0.51772687, 0.81850696],
147 [-0.79813178, 0.58842606, 0.12938478],
148 [-0.54861681, -0.6210585 , -0.55973739]])
149
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
150 """
e1c61f1 @rgommers MAINT: disable sort keyword of linalg.qz, due to Windows issues. See…
rgommers authored
151 if sort is not None:
152 # Disabled due to segfaults on win32, see ticket 1717.
153 raise ValueError("The 'sort' input of qz() has to be None (will "
154 " change when this functionality is made more robust).")
155
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
156 if not output in ['real','complex','r','c']:
157 raise ValueError("argument must be 'real', or 'complex'")
158
159 a1 = asarray_chkfinite(A)
160 b1 = asarray_chkfinite(B)
161
162 a_m, a_n = a1.shape
163 b_m, b_n = b1.shape
164 try:
165 assert a_m == a_n == b_m == b_n
166 except AssertionError:
167 raise ValueError("Array dimensions must be square and agree")
168
169 typa = a1.dtype.char
170 if output in ['complex', 'c'] and typa not in ['F','D']:
171 if typa in _double_precision:
172 a1 = a1.astype('D')
173 typa = 'D'
174 else:
175 a1 = a1.astype('F')
176 typa = 'F'
177 typb = b1.dtype.char
178 if output in ['complex', 'c'] and typb not in ['F','D']:
179 if typb in _double_precision:
180 b1 = b1.astype('D')
181 typb = 'D'
182 else:
183 b1 = b1.astype('F')
184 typb = 'F'
185
186 overwrite_a = overwrite_a or (_datacopied(a1,A))
187 overwrite_b = overwrite_b or (_datacopied(b1,B))
188
189 gges, = get_lapack_funcs(('gges',), (a1,b1))
190
191 if lwork is None or lwork == -1:
192 # get optimal work array size
193 result = gges(lambda x: None, a1, b1, lwork=-1)
194 lwork = result[-2][0].real.astype(np.int)
195
196 if sort is None:
197 sort_t = 0
198 sfunction = lambda x : None
a394b8a @jseabold STY: Remove outdated comment
jseabold authored
199 else:
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
200 sort_t = 1
201 sfunction = _select_function(sort, typa)
202
203 result = gges(sfunction, a1, b1, lwork=lwork, overwrite_a=overwrite_a,
204 overwrite_b=overwrite_b, sort_t=sort_t)
205
206 info = result[-1]
207 if info < 0:
208 raise ValueError("Illegal value in argument %d of gges" % -info)
209 elif info > 0 and info <= a_n:
210 warnings.warn("The QZ iteration failed. (a,b) are not in Schur "
e1c61f1 @rgommers MAINT: disable sort keyword of linalg.qz, due to Windows issues. See…
rgommers authored
211 "form, but ALPHAR(j), ALPHAI(j), and BETA(j) should be correct "
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
212 "for J=%d,...,N" % info-1, UserWarning)
213 elif info == a_n+1:
214 raise LinAlgError("Something other than QZ iteration failed")
215 elif info == a_n+2:
e1c61f1 @rgommers MAINT: disable sort keyword of linalg.qz, due to Windows issues. See…
rgommers authored
216 raise LinAlgError("After reordering, roundoff changed values of some "
217 "complex eigenvalues so that leading eigenvalues in the "
218 "Generalized Schur form no longer satisfy sort=True. "
356c90b @jseabold ENH: Add QZ decomposition with tests
jseabold authored
219 "This could also be caused due to scaling.")
220 elif info == a_n+3:
221 raise LinAlgError("Reordering failed in <s,d,c,z>tgsen")
222
223 # output for real
224 #AA, BB, sdim, alphar, alphai, beta, vsl, vsr, work, info
225 # output for complex
226 #AA, BB, sdim, alphai, beta, vsl, vsr, work, info
227 if sort_t == 0:
228 return result[0], result[1], result[-4], result[-3]
229 else:
230 return result[0], result[1], result[-4], result[-3], result[2]
Something went wrong with that request. Please try again.