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"""SVD decomposition functions.""" 

2import numpy 

3from numpy import zeros, r_, diag, dot, arccos, arcsin, where, clip 

4 

5# Local imports. 

6from .misc import LinAlgError, _datacopied 

7from .lapack import get_lapack_funcs, _compute_lwork 

8from .decomp import _asarray_validated 

9 

10__all__ = ['svd', 'svdvals', 'diagsvd', 'orth', 'subspace_angles', 'null_space'] 

11 

12 

13def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, 

14 check_finite=True, lapack_driver='gesdd'): 

15 """ 

16 Singular Value Decomposition. 

17 

18 Factorizes the matrix `a` into two unitary matrices ``U`` and ``Vh``, and 

19 a 1-D array ``s`` of singular values (real, non-negative) such that 

20 ``a == U @ S @ Vh``, where ``S`` is a suitably shaped matrix of zeros with 

21 main diagonal ``s``. 

22 

23 Parameters 

24 ---------- 

25 a : (M, N) array_like 

26 Matrix to decompose. 

27 full_matrices : bool, optional 

28 If True (default), `U` and `Vh` are of shape ``(M, M)``, ``(N, N)``. 

29 If False, the shapes are ``(M, K)`` and ``(K, N)``, where 

30 ``K = min(M, N)``. 

31 compute_uv : bool, optional 

32 Whether to compute also ``U`` and ``Vh`` in addition to ``s``. 

33 Default is True. 

34 overwrite_a : bool, optional 

35 Whether to overwrite `a`; may improve performance. 

36 Default is False. 

37 check_finite : bool, optional 

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

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

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

41 lapack_driver : {'gesdd', 'gesvd'}, optional 

42 Whether to use the more efficient divide-and-conquer approach 

43 (``'gesdd'``) or general rectangular approach (``'gesvd'``) 

44 to compute the SVD. MATLAB and Octave use the ``'gesvd'`` approach. 

45 Default is ``'gesdd'``. 

46 

47 .. versionadded:: 0.18 

48 

49 Returns 

50 ------- 

51 U : ndarray 

52 Unitary matrix having left singular vectors as columns. 

53 Of shape ``(M, M)`` or ``(M, K)``, depending on `full_matrices`. 

54 s : ndarray 

55 The singular values, sorted in non-increasing order. 

56 Of shape (K,), with ``K = min(M, N)``. 

57 Vh : ndarray 

58 Unitary matrix having right singular vectors as rows. 

59 Of shape ``(N, N)`` or ``(K, N)`` depending on `full_matrices`. 

60 

61 For ``compute_uv=False``, only ``s`` is returned. 

62 

63 Raises 

64 ------ 

65 LinAlgError 

66 If SVD computation does not converge. 

67 

68 See also 

69 -------- 

70 svdvals : Compute singular values of a matrix. 

71 diagsvd : Construct the Sigma matrix, given the vector s. 

72 

73 Examples 

74 -------- 

75 >>> from scipy import linalg 

76 >>> m, n = 9, 6 

77 >>> a = np.random.randn(m, n) + 1.j*np.random.randn(m, n) 

78 >>> U, s, Vh = linalg.svd(a) 

79 >>> U.shape, s.shape, Vh.shape 

80 ((9, 9), (6,), (6, 6)) 

81 

82 Reconstruct the original matrix from the decomposition: 

83 

84 >>> sigma = np.zeros((m, n)) 

85 >>> for i in range(min(m, n)): 

86 ... sigma[i, i] = s[i] 

87 >>> a1 = np.dot(U, np.dot(sigma, Vh)) 

88 >>> np.allclose(a, a1) 

89 True 

90 

91 Alternatively, use ``full_matrices=False`` (notice that the shape of 

92 ``U`` is then ``(m, n)`` instead of ``(m, m)``): 

93 

94 >>> U, s, Vh = linalg.svd(a, full_matrices=False) 

95 >>> U.shape, s.shape, Vh.shape 

96 ((9, 6), (6,), (6, 6)) 

97 >>> S = np.diag(s) 

98 >>> np.allclose(a, np.dot(U, np.dot(S, Vh))) 

99 True 

100 

101 >>> s2 = linalg.svd(a, compute_uv=False) 

102 >>> np.allclose(s, s2) 

103 True 

104 

105 """ 

106 a1 = _asarray_validated(a, check_finite=check_finite) 

107 if len(a1.shape) != 2: 

108 raise ValueError('expected matrix') 

109 m, n = a1.shape 

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

111 

112 if not isinstance(lapack_driver, str): 

113 raise TypeError('lapack_driver must be a string') 

114 if lapack_driver not in ('gesdd', 'gesvd'): 

115 raise ValueError('lapack_driver must be "gesdd" or "gesvd", not "%s"' 

116 % (lapack_driver,)) 

117 funcs = (lapack_driver, lapack_driver + '_lwork') 

118 gesXd, gesXd_lwork = get_lapack_funcs(funcs, (a1,)) 

119 

120 # compute optimal lwork 

121 lwork = _compute_lwork(gesXd_lwork, a1.shape[0], a1.shape[1], 

122 compute_uv=compute_uv, full_matrices=full_matrices) 

123 

124 # perform decomposition 

125 u, s, v, info = gesXd(a1, compute_uv=compute_uv, lwork=lwork, 

126 full_matrices=full_matrices, overwrite_a=overwrite_a) 

127 

128 if info > 0: 

129 raise LinAlgError("SVD did not converge") 

130 if info < 0: 

131 raise ValueError('illegal value in %dth argument of internal gesdd' 

132 % -info) 

133 if compute_uv: 

134 return u, s, v 

135 else: 

136 return s 

137 

138 

139def svdvals(a, overwrite_a=False, check_finite=True): 

140 """ 

141 Compute singular values of a matrix. 

142 

143 Parameters 

144 ---------- 

145 a : (M, N) array_like 

146 Matrix to decompose. 

147 overwrite_a : bool, optional 

148 Whether to overwrite `a`; may improve performance. 

149 Default is False. 

150 check_finite : bool, optional 

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

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

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

154 

155 Returns 

156 ------- 

157 s : (min(M, N),) ndarray 

158 The singular values, sorted in decreasing order. 

159 

160 Raises 

161 ------ 

162 LinAlgError 

163 If SVD computation does not converge. 

164 

165 Notes 

166 ----- 

167 ``svdvals(a)`` only differs from ``svd(a, compute_uv=False)`` by its 

168 handling of the edge case of empty ``a``, where it returns an 

169 empty sequence: 

170 

171 >>> a = np.empty((0, 2)) 

172 >>> from scipy.linalg import svdvals 

173 >>> svdvals(a) 

174 array([], dtype=float64) 

175 

176 See Also 

177 -------- 

178 svd : Compute the full singular value decomposition of a matrix. 

179 diagsvd : Construct the Sigma matrix, given the vector s. 

180 

181 Examples 

182 -------- 

183 >>> from scipy.linalg import svdvals 

184 >>> m = np.array([[1.0, 0.0], 

185 ... [2.0, 3.0], 

186 ... [1.0, 1.0], 

187 ... [0.0, 2.0], 

188 ... [1.0, 0.0]]) 

189 >>> svdvals(m) 

190 array([ 4.28091555, 1.63516424]) 

191 

192 We can verify the maximum singular value of `m` by computing the maximum 

193 length of `m.dot(u)` over all the unit vectors `u` in the (x,y) plane. 

194 We approximate "all" the unit vectors with a large sample. Because 

195 of linearity, we only need the unit vectors with angles in [0, pi]. 

196 

197 >>> t = np.linspace(0, np.pi, 2000) 

198 >>> u = np.array([np.cos(t), np.sin(t)]) 

199 >>> np.linalg.norm(m.dot(u), axis=0).max() 

200 4.2809152422538475 

201 

202 `p` is a projection matrix with rank 1. With exact arithmetic, 

203 its singular values would be [1, 0, 0, 0]. 

204 

205 >>> v = np.array([0.1, 0.3, 0.9, 0.3]) 

206 >>> p = np.outer(v, v) 

207 >>> svdvals(p) 

208 array([ 1.00000000e+00, 2.02021698e-17, 1.56692500e-17, 

209 8.15115104e-34]) 

210 

211 The singular values of an orthogonal matrix are all 1. Here, we 

212 create a random orthogonal matrix by using the `rvs()` method of 

213 `scipy.stats.ortho_group`. 

214 

215 >>> from scipy.stats import ortho_group 

216 >>> np.random.seed(123) 

217 >>> orth = ortho_group.rvs(4) 

218 >>> svdvals(orth) 

219 array([ 1., 1., 1., 1.]) 

220 

221 """ 

222 a = _asarray_validated(a, check_finite=check_finite) 

223 if a.size: 

224 return svd(a, compute_uv=0, overwrite_a=overwrite_a, 

225 check_finite=False) 

226 elif len(a.shape) != 2: 

227 raise ValueError('expected matrix') 

228 else: 

229 return numpy.empty(0) 

230 

231 

232def diagsvd(s, M, N): 

233 """ 

234 Construct the sigma matrix in SVD from singular values and size M, N. 

235 

236 Parameters 

237 ---------- 

238 s : (M,) or (N,) array_like 

239 Singular values 

240 M : int 

241 Size of the matrix whose singular values are `s`. 

242 N : int 

243 Size of the matrix whose singular values are `s`. 

244 

245 Returns 

246 ------- 

247 S : (M, N) ndarray 

248 The S-matrix in the singular value decomposition 

249 

250 See Also 

251 -------- 

252 svd : Singular value decomposition of a matrix 

253 svdvals : Compute singular values of a matrix. 

254 

255 Examples 

256 -------- 

257 >>> from scipy.linalg import diagsvd 

258 >>> vals = np.array([1, 2, 3]) # The array representing the computed svd 

259 >>> diagsvd(vals, 3, 4) 

260 array([[1, 0, 0, 0], 

261 [0, 2, 0, 0], 

262 [0, 0, 3, 0]]) 

263 >>> diagsvd(vals, 4, 3) 

264 array([[1, 0, 0], 

265 [0, 2, 0], 

266 [0, 0, 3], 

267 [0, 0, 0]]) 

268 

269 """ 

270 part = diag(s) 

271 typ = part.dtype.char 

272 MorN = len(s) 

273 if MorN == M: 

274 return r_['-1', part, zeros((M, N-M), typ)] 

275 elif MorN == N: 

276 return r_[part, zeros((M-N, N), typ)] 

277 else: 

278 raise ValueError("Length of s must be M or N.") 

279 

280 

281# Orthonormal decomposition 

282 

283def orth(A, rcond=None): 

284 """ 

285 Construct an orthonormal basis for the range of A using SVD 

286 

287 Parameters 

288 ---------- 

289 A : (M, N) array_like 

290 Input array 

291 rcond : float, optional 

292 Relative condition number. Singular values ``s`` smaller than 

293 ``rcond * max(s)`` are considered zero. 

294 Default: floating point eps * max(M,N). 

295 

296 Returns 

297 ------- 

298 Q : (M, K) ndarray 

299 Orthonormal basis for the range of A. 

300 K = effective rank of A, as determined by rcond 

301 

302 See also 

303 -------- 

304 svd : Singular value decomposition of a matrix 

305 null_space : Matrix null space 

306 

307 Examples 

308 -------- 

309 >>> from scipy.linalg import orth 

310 >>> A = np.array([[2, 0, 0], [0, 5, 0]]) # rank 2 array 

311 >>> orth(A) 

312 array([[0., 1.], 

313 [1., 0.]]) 

314 >>> orth(A.T) 

315 array([[0., 1.], 

316 [1., 0.], 

317 [0., 0.]]) 

318 

319 """ 

320 u, s, vh = svd(A, full_matrices=False) 

321 M, N = u.shape[0], vh.shape[1] 

322 if rcond is None: 

323 rcond = numpy.finfo(s.dtype).eps * max(M, N) 

324 tol = numpy.amax(s) * rcond 

325 num = numpy.sum(s > tol, dtype=int) 

326 Q = u[:, :num] 

327 return Q 

328 

329 

330def null_space(A, rcond=None): 

331 """ 

332 Construct an orthonormal basis for the null space of A using SVD 

333 

334 Parameters 

335 ---------- 

336 A : (M, N) array_like 

337 Input array 

338 rcond : float, optional 

339 Relative condition number. Singular values ``s`` smaller than 

340 ``rcond * max(s)`` are considered zero. 

341 Default: floating point eps * max(M,N). 

342 

343 Returns 

344 ------- 

345 Z : (N, K) ndarray 

346 Orthonormal basis for the null space of A. 

347 K = dimension of effective null space, as determined by rcond 

348 

349 See also 

350 -------- 

351 svd : Singular value decomposition of a matrix 

352 orth : Matrix range 

353 

354 Examples 

355 -------- 

356 1-D null space: 

357 

358 >>> from scipy.linalg import null_space 

359 >>> A = np.array([[1, 1], [1, 1]]) 

360 >>> ns = null_space(A) 

361 >>> ns * np.sign(ns[0,0]) # Remove the sign ambiguity of the vector 

362 array([[ 0.70710678], 

363 [-0.70710678]]) 

364 

365 2-D null space: 

366 

367 >>> B = np.random.rand(3, 5) 

368 >>> Z = null_space(B) 

369 >>> Z.shape 

370 (5, 2) 

371 >>> np.allclose(B.dot(Z), 0) 

372 True 

373 

374 The basis vectors are orthonormal (up to rounding error): 

375 

376 >>> Z.T.dot(Z) 

377 array([[ 1.00000000e+00, 6.92087741e-17], 

378 [ 6.92087741e-17, 1.00000000e+00]]) 

379 

380 """ 

381 u, s, vh = svd(A, full_matrices=True) 

382 M, N = u.shape[0], vh.shape[1] 

383 if rcond is None: 

384 rcond = numpy.finfo(s.dtype).eps * max(M, N) 

385 tol = numpy.amax(s) * rcond 

386 num = numpy.sum(s > tol, dtype=int) 

387 Q = vh[num:,:].T.conj() 

388 return Q 

389 

390 

391def subspace_angles(A, B): 

392 r""" 

393 Compute the subspace angles between two matrices. 

394 

395 Parameters 

396 ---------- 

397 A : (M, N) array_like 

398 The first input array. 

399 B : (M, K) array_like 

400 The second input array. 

401 

402 Returns 

403 ------- 

404 angles : ndarray, shape (min(N, K),) 

405 The subspace angles between the column spaces of `A` and `B` in 

406 descending order. 

407 

408 See Also 

409 -------- 

410 orth 

411 svd 

412 

413 Notes 

414 ----- 

415 This computes the subspace angles according to the formula 

416 provided in [1]_. For equivalence with MATLAB and Octave behavior, 

417 use ``angles[0]``. 

418 

419 .. versionadded:: 1.0 

420 

421 References 

422 ---------- 

423 .. [1] Knyazev A, Argentati M (2002) Principal Angles between Subspaces 

424 in an A-Based Scalar Product: Algorithms and Perturbation 

425 Estimates. SIAM J. Sci. Comput. 23:2008-2040. 

426 

427 Examples 

428 -------- 

429 An Hadamard matrix, which has orthogonal columns, so we expect that 

430 the suspace angle to be :math:`\frac{\pi}{2}`: 

431 

432 >>> from scipy.linalg import hadamard, subspace_angles 

433 >>> H = hadamard(4) 

434 >>> print(H) 

435 [[ 1 1 1 1] 

436 [ 1 -1 1 -1] 

437 [ 1 1 -1 -1] 

438 [ 1 -1 -1 1]] 

439 >>> np.rad2deg(subspace_angles(H[:, :2], H[:, 2:])) 

440 array([ 90., 90.]) 

441 

442 And the subspace angle of a matrix to itself should be zero: 

443 

444 >>> subspace_angles(H[:, :2], H[:, :2]) <= 2 * np.finfo(float).eps 

445 array([ True, True], dtype=bool) 

446 

447 The angles between non-orthogonal subspaces are in between these extremes: 

448 

449 >>> x = np.random.RandomState(0).randn(4, 3) 

450 >>> np.rad2deg(subspace_angles(x[:, :2], x[:, [2]])) 

451 array([ 55.832]) 

452 """ 

453 # Steps here omit the U and V calculation steps from the paper 

454 

455 # 1. Compute orthonormal bases of column-spaces 

456 A = _asarray_validated(A, check_finite=True) 

457 if len(A.shape) != 2: 

458 raise ValueError('expected 2D array, got shape %s' % (A.shape,)) 

459 QA = orth(A) 

460 del A 

461 

462 B = _asarray_validated(B, check_finite=True) 

463 if len(B.shape) != 2: 

464 raise ValueError('expected 2D array, got shape %s' % (B.shape,)) 

465 if len(B) != len(QA): 

466 raise ValueError('A and B must have the same number of rows, got ' 

467 '%s and %s' % (QA.shape[0], B.shape[0])) 

468 QB = orth(B) 

469 del B 

470 

471 # 2. Compute SVD for cosine 

472 QA_H_QB = dot(QA.T.conj(), QB) 

473 sigma = svdvals(QA_H_QB) 

474 

475 # 3. Compute matrix B 

476 if QA.shape[1] >= QB.shape[1]: 

477 B = QB - dot(QA, QA_H_QB) 

478 else: 

479 B = QA - dot(QB, QA_H_QB.T.conj()) 

480 del QA, QB, QA_H_QB 

481 

482 # 4. Compute SVD for sine 

483 mask = sigma ** 2 >= 0.5 

484 if mask.any(): 

485 mu_arcsin = arcsin(clip(svdvals(B, overwrite_a=True), -1., 1.)) 

486 else: 

487 mu_arcsin = 0. 

488 

489 # 5. Compute the principal angles 

490 # with reverse ordering of sigma because smallest sigma belongs to largest 

491 # angle theta 

492 theta = where(mask, mu_arcsin, arccos(clip(sigma[::-1], -1., 1.))) 

493 return theta