Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp

1"""SVD decomposition functions."""

2import numpy

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

5# Local imports.

6from ._misc import LinAlgError, _datacopied

7from .lapack import get_lapack_funcs, _compute_lwork

8from ._decomp import _asarray_validated

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

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

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

15 """

16 Singular Value Decomposition.

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``.

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'``.

47 .. versionadded:: 0.18

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`.

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

63 Raises

64 ------

65 LinAlgError

66 If SVD computation does not converge.

68 See Also

69 --------

70 svdvals : Compute singular values of a matrix.

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

73 Examples

74 --------

75 >>> import numpy as np

76 >>> from scipy import linalg

77 >>> rng = np.random.default_rng()

78 >>> m, n = 9, 6

79 >>> a = rng.standard_normal((m, n)) + 1.j*rng.standard_normal((m, n))

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

81 >>> U.shape, s.shape, Vh.shape

82 ((9, 9), (6,), (6, 6))

84 Reconstruct the original matrix from the decomposition:

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

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

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

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

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

91 True

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

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

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

97 >>> U.shape, s.shape, Vh.shape

98 ((9, 6), (6,), (6, 6))

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

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

101 True

102

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

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

105 True

106

107 """

108 a1 = _asarray_validated(a, check_finite=check_finite)

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

110 raise ValueError('expected matrix')

111 m, n = a1.shape

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

113

114 if not isinstance(lapack_driver, str):

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

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

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

118 % (lapack_driver,))

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

120 gesXd, gesXd_lwork = get_lapack_funcs(funcs, (a1,), ilp64='preferred')

121

122 # compute optimal lwork

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

124 compute_uv=compute_uv, full_matrices=full_matrices)

125

126 # perform decomposition

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

128 full_matrices=full_matrices, overwrite_a=overwrite_a)

129

130 if info > 0:

131 raise LinAlgError("SVD did not converge")

132 if info < 0:

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

134 % -info)

135 if compute_uv:

136 return u, s, v

137 else:

138 return s

139

140

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

142 """

143 Compute singular values of a matrix.

144

145 Parameters

146 ----------

147 a : (M, N) array_like

148 Matrix to decompose.

149 overwrite_a : bool, optional

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

151 Default is False.

152 check_finite : bool, optional

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

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

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

156

157 Returns

158 -------

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

160 The singular values, sorted in decreasing order.

161

162 Raises

163 ------

164 LinAlgError

165 If SVD computation does not converge.

166

167 See Also

168 --------

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

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

171

172 Notes

173 -----

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

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

176 empty sequence:

177

178 >>> import numpy as np

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

180 >>> from scipy.linalg import svdvals

181 >>> svdvals(a)

182 array([], dtype=float64)

183

184 Examples

185 --------

186 >>> import numpy as np

187 >>> from scipy.linalg import svdvals

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

189 ... [2.0, 3.0],

190 ... [1.0, 1.0],

191 ... [0.0, 2.0],

192 ... [1.0, 0.0]])

193 >>> svdvals(m)

194 array([ 4.28091555, 1.63516424])

195

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

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

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

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

200

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

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

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

204 4.2809152422538475

205

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

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

208

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

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

211 >>> svdvals(p)

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

213 8.15115104e-34])

214

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

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

217 `scipy.stats.ortho_group`.

218

219 >>> from scipy.stats import ortho_group

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

221 >>> svdvals(orth)

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

223

224 """

225 a = _asarray_validated(a, check_finite=check_finite)

226 if a.size:

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

228 check_finite=False)

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

230 raise ValueError('expected matrix')

231 else:

232 return numpy.empty(0)

233

234

235def diagsvd(s, M, N):

236 """

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

238

239 Parameters

240 ----------

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

242 Singular values

243 M : int

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

245 N : int

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

247

248 Returns

249 -------

250 S : (M, N) ndarray

251 The S-matrix in the singular value decomposition

252

253 See Also

254 --------

255 svd : Singular value decomposition of a matrix

256 svdvals : Compute singular values of a matrix.

257

258 Examples

259 --------

260 >>> import numpy as np

261 >>> from scipy.linalg import diagsvd

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

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

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

265 [0, 2, 0, 0],

266 [0, 0, 3, 0]])

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

268 array([[1, 0, 0],

269 [0, 2, 0],

270 [0, 0, 3],

271 [0, 0, 0]])

272

273 """

274 part = diag(s)

275 typ = part.dtype.char

276 MorN = len(s)

277 if MorN == M:

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

279 elif MorN == N:

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

281 else:

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

283

284

285# Orthonormal decomposition

286

287def orth(A, rcond=None):

288 """

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

290

291 Parameters

292 ----------

293 A : (M, N) array_like

294 Input array

295 rcond : float, optional

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

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

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

299

300 Returns

301 -------

302 Q : (M, K) ndarray

303 Orthonormal basis for the range of A.

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

305

306 See Also

307 --------

308 svd : Singular value decomposition of a matrix

309 null_space : Matrix null space

310

311 Examples

312 --------

313 >>> import numpy as np

314 >>> from scipy.linalg import orth

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

316 >>> orth(A)

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

318 [1., 0.]])

319 >>> orth(A.T)

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

321 [1., 0.],

322 [0., 0.]])

323

324 """

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

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

327 if rcond is None:

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

329 tol = numpy.amax(s) * rcond

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

331 Q = u[:, :num]

332 return Q

333

334

335def null_space(A, rcond=None):

336 """

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

338

339 Parameters

340 ----------

341 A : (M, N) array_like

342 Input array

343 rcond : float, optional

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

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

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

347

348 Returns

349 -------

350 Z : (N, K) ndarray

351 Orthonormal basis for the null space of A.

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

353

354 See Also

355 --------

356 svd : Singular value decomposition of a matrix

357 orth : Matrix range

358

359 Examples

360 --------

361 1-D null space:

362

363 >>> import numpy as np

364 >>> from scipy.linalg import null_space

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

366 >>> ns = null_space(A)

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

368 array([[ 0.70710678],

369 [-0.70710678]])

370

371 2-D null space:

372

373 >>> from numpy.random import default_rng

374 >>> rng = default_rng()

375 >>> B = rng.random((3, 5))

376 >>> Z = null_space(B)

377 >>> Z.shape

378 (5, 2)

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

380 True

381

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

383

384 >>> Z.T.dot(Z)

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

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

387

388 """

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

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

391 if rcond is None:

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

393 tol = numpy.amax(s) * rcond

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

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

396 return Q

397

398

399def subspace_angles(A, B):

400 r"""

401 Compute the subspace angles between two matrices.

402

403 Parameters

404 ----------

405 A : (M, N) array_like

406 The first input array.

407 B : (M, K) array_like

408 The second input array.

409

410 Returns

411 -------

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

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

414 descending order.

415

416 See Also

417 --------

418 orth

419 svd

420

421 Notes

422 -----

423 This computes the subspace angles according to the formula

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

425 use ``angles[0]``.

426

427 .. versionadded:: 1.0

428

429 References

430 ----------

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

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

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

434

435 Examples

436 --------

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

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

439

440 >>> import numpy as np

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

442 >>> rng = np.random.default_rng()

443 >>> H = hadamard(4)

444 >>> print(H)

445 [[ 1 1 1 1]

446 [ 1 -1 1 -1]

447 [ 1 1 -1 -1]

448 [ 1 -1 -1 1]]

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

450 array([ 90., 90.])

451

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

453

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

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

456

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

458

459 >>> x = rng.standard_normal((4, 3))

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

461 array([ 55.832]) # random

462 """

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

464

465 # 1. Compute orthonormal bases of column-spaces

466 A = _asarray_validated(A, check_finite=True)

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

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

469 QA = orth(A)

470 del A

471

472 B = _asarray_validated(B, check_finite=True)

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

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

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

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

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

478 QB = orth(B)

479 del B

480

481 # 2. Compute SVD for cosine

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

483 sigma = svdvals(QA_H_QB)

484

485 # 3. Compute matrix B

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

487 B = QB - dot(QA, QA_H_QB)

488 else:

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

490 del QA, QB, QA_H_QB

491

492 # 4. Compute SVD for sine

493 mask = sigma ** 2 >= 0.5

494 if mask.any():

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

496 else:

497 mu_arcsin = 0.

498

499 # 5. Compute the principal angles

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

501 # angle theta

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

503 return theta

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_svd.py: 14%

86 statements