Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Schur decomposition functions.""" 

2import numpy 

3from numpy import asarray_chkfinite, single, asarray, array 

4from numpy.linalg import norm 

5 

6 

7# Local imports. 

8from .misc import LinAlgError, _datacopied 

9from .lapack import get_lapack_funcs 

10from .decomp import eigvals 

11 

12__all__ = ['schur', 'rsf2csf'] 

13 

14_double_precision = ['i', 'l', 'd'] 

15 

16 

17def schur(a, output='real', lwork=None, overwrite_a=False, sort=None, 

18 check_finite=True): 

19 """ 

20 Compute Schur decomposition of a matrix. 

21 

22 The Schur decomposition is:: 

23 

24 A = Z T Z^H 

25 

26 where Z is unitary and T is either upper-triangular, or for real 

27 Schur decomposition (output='real'), quasi-upper triangular. In 

28 the quasi-triangular form, 2x2 blocks describing complex-valued 

29 eigenvalue pairs may extrude from the diagonal. 

30 

31 Parameters 

32 ---------- 

33 a : (M, M) array_like 

34 Matrix to decompose 

35 output : {'real', 'complex'}, optional 

36 Construct the real or complex Schur decomposition (for real matrices). 

37 lwork : int, optional 

38 Work array size. If None or -1, it is automatically computed. 

39 overwrite_a : bool, optional 

40 Whether to overwrite data in a (may improve performance). 

41 sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional 

42 Specifies whether the upper eigenvalues should be sorted. A callable 

43 may be passed that, given a eigenvalue, returns a boolean denoting 

44 whether the eigenvalue should be sorted to the top-left (True). 

45 Alternatively, string parameters may be used:: 

46 

47 'lhp' Left-hand plane (x.real < 0.0) 

48 'rhp' Right-hand plane (x.real > 0.0) 

49 'iuc' Inside the unit circle (x*x.conjugate() <= 1.0) 

50 'ouc' Outside the unit circle (x*x.conjugate() > 1.0) 

51 

52 Defaults to None (no sorting). 

53 check_finite : bool, optional 

54 Whether to check that the input matrix contains only finite numbers. 

55 Disabling may give a performance gain, but may result in problems 

56 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

57 

58 Returns 

59 ------- 

60 T : (M, M) ndarray 

61 Schur form of A. It is real-valued for the real Schur decomposition. 

62 Z : (M, M) ndarray 

63 An unitary Schur transformation matrix for A. 

64 It is real-valued for the real Schur decomposition. 

65 sdim : int 

66 If and only if sorting was requested, a third return value will 

67 contain the number of eigenvalues satisfying the sort condition. 

68 

69 Raises 

70 ------ 

71 LinAlgError 

72 Error raised under three conditions: 

73 

74 1. The algorithm failed due to a failure of the QR algorithm to 

75 compute all eigenvalues. 

76 2. If eigenvalue sorting was requested, the eigenvalues could not be 

77 reordered due to a failure to separate eigenvalues, usually because 

78 of poor conditioning. 

79 3. If eigenvalue sorting was requested, roundoff errors caused the 

80 leading eigenvalues to no longer satisfy the sorting condition. 

81 

82 See also 

83 -------- 

84 rsf2csf : Convert real Schur form to complex Schur form 

85 

86 Examples 

87 -------- 

88 >>> from scipy.linalg import schur, eigvals 

89 >>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]]) 

90 >>> T, Z = schur(A) 

91 >>> T 

92 array([[ 2.65896708, 1.42440458, -1.92933439], 

93 [ 0. , -0.32948354, -0.49063704], 

94 [ 0. , 1.31178921, -0.32948354]]) 

95 >>> Z 

96 array([[0.72711591, -0.60156188, 0.33079564], 

97 [0.52839428, 0.79801892, 0.28976765], 

98 [0.43829436, 0.03590414, -0.89811411]]) 

99 

100 >>> T2, Z2 = schur(A, output='complex') 

101 >>> T2 

102 array([[ 2.65896708, -1.22839825+1.32378589j, 0.42590089+1.51937378j], 

103 [ 0. , -0.32948354+0.80225456j, -0.59877807+0.56192146j], 

104 [ 0. , 0. , -0.32948354-0.80225456j]]) 

105 >>> eigvals(T2) 

106 array([2.65896708, -0.32948354+0.80225456j, -0.32948354-0.80225456j]) 

107 

108 An arbitrary custom eig-sorting condition, having positive imaginary part, 

109 which is satisfied by only one eigenvalue 

110 

111 >>> T3, Z3, sdim = schur(A, output='complex', sort=lambda x: x.imag > 0) 

112 >>> sdim 

113 1 

114 

115 """ 

116 if output not in ['real', 'complex', 'r', 'c']: 

117 raise ValueError("argument must be 'real', or 'complex'") 

118 if check_finite: 

119 a1 = asarray_chkfinite(a) 

120 else: 

121 a1 = asarray(a) 

122 if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]): 

123 raise ValueError('expected square matrix') 

124 typ = a1.dtype.char 

125 if output in ['complex', 'c'] and typ not in ['F', 'D']: 

126 if typ in _double_precision: 

127 a1 = a1.astype('D') 

128 typ = 'D' 

129 else: 

130 a1 = a1.astype('F') 

131 typ = 'F' 

132 overwrite_a = overwrite_a or (_datacopied(a1, a)) 

133 gees, = get_lapack_funcs(('gees',), (a1,)) 

134 if lwork is None or lwork == -1: 

135 # get optimal work array 

136 result = gees(lambda x: None, a1, lwork=-1) 

137 lwork = result[-2][0].real.astype(numpy.int) 

138 

139 if sort is None: 

140 sort_t = 0 

141 sfunction = lambda x: None 

142 else: 

143 sort_t = 1 

144 if callable(sort): 

145 sfunction = sort 

146 elif sort == 'lhp': 

147 sfunction = lambda x: (x.real < 0.0) 

148 elif sort == 'rhp': 

149 sfunction = lambda x: (x.real >= 0.0) 

150 elif sort == 'iuc': 

151 sfunction = lambda x: (abs(x) <= 1.0) 

152 elif sort == 'ouc': 

153 sfunction = lambda x: (abs(x) > 1.0) 

154 else: 

155 raise ValueError("'sort' parameter must either be 'None', or a " 

156 "callable, or one of ('lhp','rhp','iuc','ouc')") 

157 

158 result = gees(sfunction, a1, lwork=lwork, overwrite_a=overwrite_a, 

159 sort_t=sort_t) 

160 

161 info = result[-1] 

162 if info < 0: 

163 raise ValueError('illegal value in {}-th argument of internal gees' 

164 ''.format(-info)) 

165 elif info == a1.shape[0] + 1: 

166 raise LinAlgError('Eigenvalues could not be separated for reordering.') 

167 elif info == a1.shape[0] + 2: 

168 raise LinAlgError('Leading eigenvalues do not satisfy sort condition.') 

169 elif info > 0: 

170 raise LinAlgError("Schur form not found. Possibly ill-conditioned.") 

171 

172 if sort_t == 0: 

173 return result[0], result[-3] 

174 else: 

175 return result[0], result[-3], result[1] 

176 

177 

178eps = numpy.finfo(float).eps 

179feps = numpy.finfo(single).eps 

180 

181_array_kind = {'b': 0, 'h': 0, 'B': 0, 'i': 0, 'l': 0, 

182 'f': 0, 'd': 0, 'F': 1, 'D': 1} 

183_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1} 

184_array_type = [['f', 'd'], ['F', 'D']] 

185 

186 

187def _commonType(*arrays): 

188 kind = 0 

189 precision = 0 

190 for a in arrays: 

191 t = a.dtype.char 

192 kind = max(kind, _array_kind[t]) 

193 precision = max(precision, _array_precision[t]) 

194 return _array_type[kind][precision] 

195 

196 

197def _castCopy(type, *arrays): 

198 cast_arrays = () 

199 for a in arrays: 

200 if a.dtype.char == type: 

201 cast_arrays = cast_arrays + (a.copy(),) 

202 else: 

203 cast_arrays = cast_arrays + (a.astype(type),) 

204 if len(cast_arrays) == 1: 

205 return cast_arrays[0] 

206 else: 

207 return cast_arrays 

208 

209 

210def rsf2csf(T, Z, check_finite=True): 

211 """ 

212 Convert real Schur form to complex Schur form. 

213 

214 Convert a quasi-diagonal real-valued Schur form to the upper-triangular 

215 complex-valued Schur form. 

216 

217 Parameters 

218 ---------- 

219 T : (M, M) array_like 

220 Real Schur form of the original array 

221 Z : (M, M) array_like 

222 Schur transformation matrix 

223 check_finite : bool, optional 

224 Whether to check that the input arrays contain only finite numbers. 

225 Disabling may give a performance gain, but may result in problems 

226 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

227 

228 Returns 

229 ------- 

230 T : (M, M) ndarray 

231 Complex Schur form of the original array 

232 Z : (M, M) ndarray 

233 Schur transformation matrix corresponding to the complex form 

234 

235 See Also 

236 -------- 

237 schur : Schur decomposition of an array 

238 

239 Examples 

240 -------- 

241 >>> from scipy.linalg import schur, rsf2csf 

242 >>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]]) 

243 >>> T, Z = schur(A) 

244 >>> T 

245 array([[ 2.65896708, 1.42440458, -1.92933439], 

246 [ 0. , -0.32948354, -0.49063704], 

247 [ 0. , 1.31178921, -0.32948354]]) 

248 >>> Z 

249 array([[0.72711591, -0.60156188, 0.33079564], 

250 [0.52839428, 0.79801892, 0.28976765], 

251 [0.43829436, 0.03590414, -0.89811411]]) 

252 >>> T2 , Z2 = rsf2csf(T, Z) 

253 >>> T2 

254 array([[2.65896708+0.j, -1.64592781+0.743164187j, -1.21516887+1.00660462j], 

255 [0.+0.j , -0.32948354+8.02254558e-01j, -0.82115218-2.77555756e-17j], 

256 [0.+0.j , 0.+0.j, -0.32948354-0.802254558j]]) 

257 >>> Z2 

258 array([[0.72711591+0.j, 0.28220393-0.31385693j, 0.51319638-0.17258824j], 

259 [0.52839428+0.j, 0.24720268+0.41635578j, -0.68079517-0.15118243j], 

260 [0.43829436+0.j, -0.76618703+0.01873251j, -0.03063006+0.46857912j]]) 

261 

262 """ 

263 if check_finite: 

264 Z, T = map(asarray_chkfinite, (Z, T)) 

265 else: 

266 Z, T = map(asarray, (Z, T)) 

267 

268 for ind, X in enumerate([Z, T]): 

269 if X.ndim != 2 or X.shape[0] != X.shape[1]: 

270 raise ValueError("Input '{}' must be square.".format('ZT'[ind])) 

271 

272 if T.shape[0] != Z.shape[0]: 

273 raise ValueError("Input array shapes must match: Z: {} vs. T: {}" 

274 "".format(Z.shape, T.shape)) 

275 N = T.shape[0] 

276 t = _commonType(Z, T, array([3.0], 'F')) 

277 Z, T = _castCopy(t, Z, T) 

278 

279 for m in range(N-1, 0, -1): 

280 if abs(T[m, m-1]) > eps*(abs(T[m-1, m-1]) + abs(T[m, m])): 

281 mu = eigvals(T[m-1:m+1, m-1:m+1]) - T[m, m] 

282 r = norm([mu[0], T[m, m-1]]) 

283 c = mu[0] / r 

284 s = T[m, m-1] / r 

285 G = array([[c.conj(), s], [-s, c]], dtype=t) 

286 

287 T[m-1:m+1, m-1:] = G.dot(T[m-1:m+1, m-1:]) 

288 T[:m+1, m-1:m+1] = T[:m+1, m-1:m+1].dot(G.conj().T) 

289 Z[:, m-1:m+1] = Z[:, m-1:m+1].dot(G.conj().T) 

290 

291 T[m, m-1] = 0.0 

292 return T, Z