Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/linalg/decomp_svd.py : 14%

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
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 >>> 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))
82 Reconstruct the original matrix from the decomposition:
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
91 Alternatively, use ``full_matrices=False`` (notice that the shape of
92 ``U`` is then ``(m, n)`` instead of ``(m, m)``):
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
101 >>> s2 = linalg.svd(a, compute_uv=False)
102 >>> np.allclose(s, s2)
103 True
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))
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,))
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)
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)
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
139def svdvals(a, overwrite_a=False, check_finite=True):
140 """
141 Compute singular values of a matrix.
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.
155 Returns
156 -------
157 s : (min(M, N),) ndarray
158 The singular values, sorted in decreasing order.
160 Raises
161 ------
162 LinAlgError
163 If SVD computation does not converge.
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:
171 >>> a = np.empty((0, 2))
172 >>> from scipy.linalg import svdvals
173 >>> svdvals(a)
174 array([], dtype=float64)
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.
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])
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].
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
202 `p` is a projection matrix with rank 1. With exact arithmetic,
203 its singular values would be [1, 0, 0, 0].
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])
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`.
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.])
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)
232def diagsvd(s, M, N):
233 """
234 Construct the sigma matrix in SVD from singular values and size M, N.
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`.
245 Returns
246 -------
247 S : (M, N) ndarray
248 The S-matrix in the singular value decomposition
250 See Also
251 --------
252 svd : Singular value decomposition of a matrix
253 svdvals : Compute singular values of a matrix.
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]])
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.")
281# Orthonormal decomposition
283def orth(A, rcond=None):
284 """
285 Construct an orthonormal basis for the range of A using SVD
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).
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
302 See also
303 --------
304 svd : Singular value decomposition of a matrix
305 null_space : Matrix null space
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.]])
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
330def null_space(A, rcond=None):
331 """
332 Construct an orthonormal basis for the null space of A using SVD
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).
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
349 See also
350 --------
351 svd : Singular value decomposition of a matrix
352 orth : Matrix range
354 Examples
355 --------
356 1-D null space:
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]])
365 2-D null space:
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
374 The basis vectors are orthonormal (up to rounding error):
376 >>> Z.T.dot(Z)
377 array([[ 1.00000000e+00, 6.92087741e-17],
378 [ 6.92087741e-17, 1.00000000e+00]])
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
391def subspace_angles(A, B):
392 r"""
393 Compute the subspace angles between two matrices.
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.
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.
408 See Also
409 --------
410 orth
411 svd
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]``.
419 .. versionadded:: 1.0
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.
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}`:
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.])
442 And the subspace angle of a matrix to itself should be zero:
444 >>> subspace_angles(H[:, :2], H[:, :2]) <= 2 * np.finfo(float).eps
445 array([ True, True], dtype=bool)
447 The angles between non-orthogonal subspaces are in between these extremes:
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
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
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
471 # 2. Compute SVD for cosine
472 QA_H_QB = dot(QA.T.conj(), QB)
473 sigma = svdvals(QA_H_QB)
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
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.
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