linalg.decomp_svd

SVD decomposition functions.

Module Contents

Functions

svd(a,full_matrices=True,compute_uv=True,overwrite_a=False,check_finite=True,lapack_driver=”gesdd”) Singular Value Decomposition.
svdvals(a,overwrite_a=False,check_finite=True) Compute singular values of a matrix.
diagsvd(s,M,N) Construct the sigma matrix in SVD from singular values and size M, N.
orth(A) Construct an orthonormal basis for the range of A using SVD
subspace_angles(A,B) r
svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver="gesdd")

Singular Value Decomposition.

Factorizes the matrix a into two unitary matrices U and Vh, and a 1-D array s of singular values (real, non-negative) such that a == U @ S @ Vh, where S is a suitably shaped matrix of zeros with main diagonal s.

a : (M, N) array_like
Matrix to decompose.
full_matrices : bool, optional
If True (default), U and Vh are of shape (M, M), (N, N). If False, the shapes are (M, K) and (K, N), where K = min(M, N).
compute_uv : bool, optional
Whether to compute also U and Vh in addition to s. Default is True.
overwrite_a : bool, optional
Whether to overwrite a; may improve performance. Default is False.
check_finite : bool, optional
Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
lapack_driver : {‘gesdd’, ‘gesvd’}, optional

Whether to use the more efficient divide-and-conquer approach ('gesdd') or general rectangular approach ('gesvd') to compute the SVD. MATLAB and Octave use the 'gesvd' approach. Default is 'gesdd'.

New in version 0.18.

U : ndarray
Unitary matrix having left singular vectors as columns. Of shape (M, M) or (M, K), depending on full_matrices.
s : ndarray
The singular values, sorted in non-increasing order. Of shape (K,), with K = min(M, N).
Vh : ndarray
Unitary matrix having right singular vectors as rows. Of shape (N, N) or (K, N) depending on full_matrices.

For compute_uv=False, only s is returned.

LinAlgError
If SVD computation does not converge.

svdvals : Compute singular values of a matrix. diagsvd : Construct the Sigma matrix, given the vector s.

>>> from scipy import linalg
>>> m, n = 9, 6
>>> a = np.random.randn(m, n) + 1.j*np.random.randn(m, n)
>>> U, s, Vh = linalg.svd(a)
>>> U.shape,  s.shape, Vh.shape
((9, 9), (6,), (6, 6))

Reconstruct the original matrix from the decomposition:

>>> sigma = np.zeros((m, n))
>>> for i in range(min(m, n)):
...     sigma[i, i] = s[i]
>>> a1 = np.dot(U, np.dot(sigma, Vh))
>>> np.allclose(a, a1)
True

Alternatively, use full_matrices=False (notice that the shape of U is then (m, n) instead of (m, m)):

>>> U, s, Vh = linalg.svd(a, full_matrices=False)
>>> U.shape, s.shape, Vh.shape
((9, 6), (6,), (6, 6))
>>> S = np.diag(s)
>>> np.allclose(a, np.dot(U, np.dot(S, Vh)))
True
>>> s2 = linalg.svd(a, compute_uv=False)
>>> np.allclose(s, s2)
True
svdvals(a, overwrite_a=False, check_finite=True)

Compute singular values of a matrix.

a : (M, N) array_like
Matrix to decompose.
overwrite_a : bool, optional
Whether to overwrite a; may improve performance. Default is False.
check_finite : bool, optional
Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
s : (min(M, N),) ndarray
The singular values, sorted in decreasing order.
LinAlgError
If SVD computation does not converge.

svdvals(a) only differs from svd(a, compute_uv=False) by its handling of the edge case of empty a, where it returns an empty sequence:

>>> a = np.empty((0, 2))
>>> from scipy.linalg import svdvals
>>> svdvals(a)
array([], dtype=float64)

svd : Compute the full singular value decomposition of a matrix. diagsvd : Construct the Sigma matrix, given the vector s.

>>> from scipy.linalg import svdvals
>>> m = np.array([[1.0, 0.0],
...               [2.0, 3.0],
...               [1.0, 1.0],
...               [0.0, 2.0],
...               [1.0, 0.0]])
>>> svdvals(m)
array([ 4.28091555,  1.63516424])

We can verify the maximum singular value of m by computing the maximum length of m.dot(u) over all the unit vectors u in the (x,y) plane. We approximate “all” the unit vectors with a large sample. Because of linearity, we only need the unit vectors with angles in [0, pi].

>>> t = np.linspace(0, np.pi, 2000)
>>> u = np.array([np.cos(t), np.sin(t)])
>>> np.linalg.norm(m.dot(u), axis=0).max()
4.2809152422538475

p is a projection matrix with rank 1. With exact arithmetic, its singular values would be [1, 0, 0, 0].

>>> v = np.array([0.1, 0.3, 0.9, 0.3])
>>> p = np.outer(v, v)
>>> svdvals(p)
array([  1.00000000e+00,   2.02021698e-17,   1.56692500e-17,
         8.15115104e-34])

The singular values of an orthogonal matrix are all 1. Here we create a random orthogonal matrix by using the rvs() method of scipy.stats.ortho_group.

>>> from scipy.stats import ortho_group
>>> np.random.seed(123)
>>> orth = ortho_group.rvs(4)
>>> svdvals(orth)
array([ 1.,  1.,  1.,  1.])
diagsvd(s, M, N)

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

s : (M,) or (N,) array_like
Singular values
M : int
Size of the matrix whose singular values are s.
N : int
Size of the matrix whose singular values are s.
S : (M, N) ndarray
The S-matrix in the singular value decomposition
orth(A)

Construct an orthonormal basis for the range of A using SVD

A : (M, N) array_like
Input array
Q : (M, K) ndarray
Orthonormal basis for the range of A. K = effective rank of A, as determined by automatic cutoff

svd : Singular value decomposition of a matrix

subspace_angles(A, B)

r Compute the subspace angles between two matrices.

A : (M, N) array_like
The first input array.
B : (M, K) array_like
The second input array.
angles : ndarray, shape (min(N, K),)
The subspace angles between the column spaces of A and B.

orth svd

This computes the subspace angles according to the formula provided in [1]. For equivalence with MATLAB and Octave behavior, use angles[0].

New in version 1.0.

[1]Knyazev A, Argentati M (2002) Principal Angles between Subspaces in an A-Based Scalar Product: Algorithms and Perturbation Estimates. SIAM J. Sci. Comput. 23:2008-2040.

A Hadamard matrix, which has orthogonal columns, so we expect that the suspace angle to be :

>>> from scipy.linalg import hadamard, subspace_angles
>>> H = hadamard(4)
>>> print(H)
[[ 1  1  1  1]
 [ 1 -1  1 -1]
 [ 1  1 -1 -1]
 [ 1 -1 -1  1]]
>>> np.rad2deg(subspace_angles(H[:, :2], H[:, 2:]))
array([ 90.,  90.])

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

>>> subspace_angles(H[:, :2], H[:, :2]) <= 2 * np.finfo(float).eps
array([ True,  True], dtype=bool)

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

>>> x = np.random.RandomState(0).randn(4, 3)
>>> np.rad2deg(subspace_angles(x[:, :2], x[:, [2]]))
array([ 55.832])