Singular value decomposition

Singular value decomposition#

The SVD of a \(m \times n\) matrix \(A\) is

\[ A = U \Sigma V^\star \]

We use numpy.linalg.svd or scipy.linalg.svd to compute this. Below we use numpy.

%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt

Construct a random \(3 \times 2\) matrix.

A = 2 * np.random.rand(3,2) - 1

Reduced SVD#

The reduced SVD is given by

U,S,Vh = np.linalg.svd(A, full_matrices=False)

The function returns \(\textrm{Vh} = V^\star\).

The singular vectors and values are given by

\[\begin{align*} u_j &= \textrm{U[:,j]}, & 0 \le j < n \\ v_j &= \textrm{Vh[j,:]}, & 0 \le j < n \\ \sigma_j &= \textrm{S[j]}, & 0 \le j < n \end{align*}\]
U,S,Vh = np.linalg.svd(A, full_matrices=False)
print('U = \n', U)
print('sigma = ', *S)
print('V = \n', Vh.T)

U = 
 [[ 0.56362633  0.81337542]
 [-0.14178377  0.2670426 ]
 [-0.81377068  0.51682557]]
sigma =  1.390363151754217 0.7636365942544742
V = 
 [[-0.54540601  0.83817199]
 [ 0.83817199  0.54540601]]

\(U\) is \(m \times n\) and \(V\) is \(n \times n\). We can verify that \(A = U \Sigma V^\star\) by computing the element-wise maximum difference

print(np.abs(A - U@np.diag(S)@Vh).max())

2.220446049250313e-16

Full SVD#

The full SVD is given by

U,S,Vh = np.linalg.svd(A, full_matrices=True)

The singular vectors and values are given by

\[\begin{align*} u_j &= \textrm{U[:,j]}, & 0 \le j < m \\ v_j &= \textrm{Vh[j,:]}, & 0 \le j < n \\ \sigma_j &= \textrm{S[j]}, & 0 \le j < n \end{align*}\]
U,S,Vh = np.linalg.svd(A, full_matrices=True)
print('U = \n', U)
print('sigma = ', *S)
print('V = \n', Vh.T)

U = 
 [[ 0.56362633  0.81337542  0.14403396]
 [-0.14178377  0.2670426  -0.95319757]
 [-0.81377068  0.51682557  0.26583567]]
sigma =  1.390363151754217 0.7636365942544742
V = 
 [[-0.54540601  0.83817199]
 [ 0.83817199  0.54540601]]

\(U\) is \(m \times m\) and \(V\) is \(n \times n\).

Rank deficient case#

Construct a random \(4 \times 3\) matrix whose last column is sum of first two columns. Its rank is 2.

A = 2 * np.random.rand(4,3) - 1
A[:,2] = A[:,0] + A[:,1]

Compute the reduced SVD

U,S,Vh = np.linalg.svd(A, full_matrices=False)
print('U = \n', U)
print('sigma = ', *S)
print('V = \n', Vh.T)

U = 
 [[-0.59953071 -0.34508729 -0.39074021]
 [ 0.09225911 -0.48661847  0.80308371]
 [ 0.41885745 -0.76715029 -0.39806246]
 [-0.67572895 -0.23579155  0.20958273]]
sigma =  3.037127455137475 0.41484360285223504 2.539493840585527e-17
V = 
 [[-0.49551087  0.64894964  0.57735027]
 [-0.31425144 -0.75359983  0.57735027]
 [-0.80976231 -0.10465018 -0.57735027]]

We notice that one out of three singular values is zero. The number of non-zero singular values corresponds to the rank of the matrix.

Random matrix#

m, n = 100, 50
A = 2 * np.random.rand(m,n) - 1
U,S,Vh = np.linalg.svd(A, full_matrices=False)

plt.subplot(121)
plt.plot(S,'o')
plt.xlabel('j')
plt.title('$\\sigma_j$')

plt.subplot(122)
plt.plot(S/S[0],'o')
plt.xlabel('j')
plt.title('$\\sigma_j/\\sigma_1$');
../_images/419dc8cc6783dfef89895be6775e7e00636585f62d70f85945907d687bc0e2f5.svg

All singular values are of very similar size.

Contents