Singular value decomposition Tobin A. Driscoll Richard J. Braun
We now introduce another factorization that is as fundamental as the EVD .
Definition 7.3.1 (Singular value decomposition (
SVD ))
The singular value decomposition m × n m\times n m × n A \mathbf{A} A
A = U S V ∗ , \mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^*, A = US V ∗ , where U ∈ C m × m \mathbf{U}\in\mathbb{C}^{m\times m} U ∈ C m × m V ∈ C n × n \mathbf{V}\in\mathbb{C}^{n\times n} V ∈ C n × n S ∈ R m × n \mathbf{S}\in\mathbb{R}^{m\times n} S ∈ R m × n
The columns of U \mathbf{U} U V \mathbf{V} V left and right singular vectors , respectively. The diagonal elements of S \mathbf{S} S σ 1 , … , σ r \sigma_1,\ldots,\sigma_r σ 1 , … , σ r r = min { m , n } r=\min\{m,n\} r = min { m , n } singular values of A \mathbf{A} A
σ 1 ≥ σ 2 ≥ ⋯ ≥ σ r ≥ 0 , r = min { m , n } . \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r\ge 0, \qquad r=\min\{m,n\}. σ 1 ≥ σ 2 ≥ ⋯ ≥ σ r ≥ 0 , r = min { m , n } . We call σ 1 \sigma_1 σ 1 principal singular value and u 1 \mathbf{u}_{1} u 1 v 1 \mathbf{v}_{1} v 1 principal singular vectors .
Every m × n m\times n m × n SVD . The singular values of a matrix are unique, but the singular vectors are not. If the matrix is real, then U \mathbf{U} U V \mathbf{V} V (7.3.1)
It is easy to check that
[ 3 4 ] = ( 1 5 [ 3 − 4 4 3 ] ) ⏟ U ⋅ [ 5 0 ] ⏟ S ⋅ [ 1 ] ⏟ V T \begin{bmatrix} 3 \\ 4 \end{bmatrix}
= \underbrace{\left(\frac{1}{5} \begin{bmatrix}
3 & -4 \\ 4 & 3
\end{bmatrix}\,\right)}_{\mathbf{U}} \cdot
\underbrace{\begin{bmatrix}
5 \\ 0
\end{bmatrix}}_{\mathbf{S}} \cdot
\underbrace{\begin{bmatrix}
1
\end{bmatrix}}_{\mathbf{V}^T} [ 3 4 ] = U ( 5 1 [ 3 4 − 4 3 ] ) ⋅ S [ 5 0 ] ⋅ V T [ 1 ] meets all the requirements of an SVD . Hence, interpreted as a 2 × 1 2\times 1 2 × 1 [ 3 , 4 ] [3,4] [ 3 , 4 ]
Suppose A \mathbf{A} A A = U S V T \mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^T A = US V T SVD . Then A T = V S T U T \mathbf{A}^T=\mathbf{V}\mathbf{S}^T\mathbf{U}^T A T = V S T U T SVD for A T \mathbf{A}^T A T A \mathbf{A} A A T \mathbf{A}^T A T
7.3.1 Connections to the EVD ¶ The eigenvalues of A ∗ A \mathbf{A}^*\mathbf{A} A ∗ A min { m , n } \min\{m,n\} min { m , n } A \mathbf{A} A
Let A = U S V ∗ \mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^* A = US V ∗ m × n m\times n m × n B = A ∗ A \mathbf{B}=\mathbf{A}^*\mathbf{A} B = A ∗ A
B = ( V S ∗ U ∗ ) ( U S V ∗ ) = V S ∗ S V ∗ = V ( S T S ) V − 1 . \mathbf{B} = (\mathbf{V}\mathbf{S}^*\mathbf{U}^*) (\mathbf{U}\mathbf{S}\mathbf{V}^*) = \mathbf{V}\mathbf{S}^*\mathbf{S}\mathbf{V}^* = \mathbf{V}(\mathbf{S}^T\mathbf{S})\mathbf{V}^{-1}. B = ( V S ∗ U ∗ ) ( US V ∗ ) = V S ∗ S V ∗ = V ( S T S ) V − 1 . Note that S T S \mathbf{S}^T\mathbf{S} S T S n × n n \times n n × n m ≥ n m \ge n m ≥ n
S T S = [ σ 1 2 ⋱ σ n 2 ] . \mathbf{S}^T\mathbf{S} =
\begin{bmatrix}
\sigma_1^2 & & \\
& \ddots & \\
& & \sigma_n^2
\end{bmatrix}. S T S = ⎣ ⎡ σ 1 2 ⋱ σ n 2 ⎦ ⎤ . On the other hand, if m < n m<n m < n
S T S = [ σ 1 2 ⋱ σ m 2 0 ] . \mathbf{S}^T\mathbf{S} =
\begin{bmatrix}
\sigma_1^2 & & & \\
& \ddots & & \\
& & \sigma_m^2 & \\
& & & \boldsymbol{0}
\end{bmatrix}. S T S = ⎣ ⎡ σ 1 2 ⋱ σ m 2 0 ⎦ ⎤ . Except for some unimportant technicalities, the eigenvectors of A ∗ A \mathbf{A}^*\mathbf{A} A ∗ A A \mathbf{A} A A V = U S \mathbf{A}\mathbf{V} = \mathbf{U}\mathbf{S} AV = US
Another close connection between EVD and SVD comes via the ( m + n ) × ( m + n ) (m+n)\times (m+n) ( m + n ) × ( m + n )
C = [ 0 A ∗ A 0 ] . \mathbf{C} =
\begin{bmatrix}
0 & \mathbf{A}^* \\ \mathbf{A} & 0
\end{bmatrix}. C = [ 0 A A ∗ 0 ] . If σ is a singular value of B \mathbf{B} B − σ -\sigma − σ C \mathbf{C} C Exercise 7.3.11 SVD .
7.3.2 Interpreting the SVD ¶ Another way to write A = U S V ∗ \mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^* A = US V ∗
A V = U S . \mathbf{A}\mathbf{V}=\mathbf{U}\mathbf{S}. AV = US . Taken columnwise, this equation means
A v k = σ k u k , k = 1 , … , r = min { m , n } . \mathbf{A} \mathbf{v}_{k} = \sigma_k \mathbf{u}_{k}, \qquad k=1,\ldots,r=\min\{m,n\}. A v k = σ k u k , k = 1 , … , r = min { m , n } . In words, each right singular vector is mapped by A \mathbf{A} A
Both the SVD and the EVD describe a matrix in terms of some special vectors and a few scalars. Table 7.3.1 SVD sacrifices having the same basis in both source and image spaces—after all, they may not even have the same dimension—but as a result gains orthogonality in both spaces.
Table 7.3.1: Comparison of the EVD and SVD
EVD SVD exists for most square matrices exists for all rectangular and square matrices A x k = λ k x k \mathbf{A}\mathbf{x}_k = \lambda_k \mathbf{x}_k A x k = λ k x k A v k = σ k u k \mathbf{A} \mathbf{v}_k = \sigma_k \mathbf{u}_k A v k = σ k u k same basis for domain and range of A \mathbf{A} A two orthonormal bases may have poor conditioning perfectly conditioned
In The QR factorization we saw that a matrix has both full and thin forms of the QR factorization. A similar situation holds with the SVD .
Suppose A \mathbf{A} A m × n m\times n m × n m > n m > n m > n A = U S V ∗ \mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^* A = US V ∗ SVD . The last m − n m-n m − n S \mathbf{S} S S \mathbf{S} S
U S = [ u 1 ⋯ u n u n + 1 ⋯ u m ] [ σ 1 ⋱ σ n 0 ] = [ u 1 ⋯ u n ] [ σ 1 ⋱ σ n ] = U ^ S ^ , \begin{align*}
\mathbf{U} \mathbf{S} & =
\begin{bmatrix}
\mathbf{u}_1 & \cdots & \mathbf{u}_n & \mathbf{u}_{n+1} & \cdots & \mathbf{u}_m
\end{bmatrix}
\begin{bmatrix}
\sigma_1 & & \\
& \ddots & \\
& & \sigma_n \\
& & \\
& \boldsymbol{0} & \\
& &
\end{bmatrix} \\
&=
\begin{bmatrix}
\mathbf{u}_1 & \cdots & \mathbf{u}_n
\end{bmatrix}
\begin{bmatrix}
\sigma_1 & & \\
& \ddots & \\
& & \sigma_n
\end{bmatrix} = \hat{\mathbf{U}} \hat{\mathbf{S}},
\end{align*} US = [ u 1 ⋯ u n u n + 1 ⋯ u m ] ⎣ ⎡ σ 1 ⋱ 0 σ n ⎦ ⎤ = [ u 1 ⋯ u n ] ⎣ ⎡ σ 1 ⋱ σ n ⎦ ⎤ = U ^ S ^ , in which U ^ \hat{\mathbf{U}} U ^ m × n m\times n m × n S ^ \hat{\mathbf{S}} S ^ n × n n\times n n × n thin SVD
A = U ^ S ^ V ∗ , \mathbf{A}=\hat{\mathbf{U}}\hat{\mathbf{S}}\mathbf{V}^*, A = U ^ S ^ V ∗ , in which S ^ \hat{\mathbf{S}} S ^ U ^ \hat{\mathbf{U}} U ^ ONC but not square.
So, in sketch form, a full SVD of a matrix that is taller than it is wide looks like
= ⋅ ⋅ \rule{1cm}{2.4cm} \; \raisebox{11mm}{=} \; \rule{2.4cm}{2.4cm} \; \raisebox{11mm}{$\centerdot$} \; \rule{1cm}{2.4cm}\; \raisebox{11mm}{$\centerdot$} \; \rule[6mm]{1cm}{1cm}\quad = ⋅ ⋅ while a thin SVD looks like
= ⋅ ⋅ \rule{1cm}{2.4cm} \; \raisebox{11mm}{=} \; \rule{1cm}{2.4cm} \; \raisebox{11mm}{$\centerdot$} \; \rule[6mm]{1cm}{1cm} \; \raisebox{11mm}{$\centerdot$} \; \rule[6mm]{1cm}{1cm} = ⋅ ⋅ The thin form retains all the information about A \mathbf{A} A SVD ; the factorization is still an equality, not an approximation. It is computationally preferable when m ≫ n m \gg n m ≫ n SVD . For a matrix with more columns than rows, one can derive a thin form by taking the adjoint of the thin SVD of A ∗ \mathbf{A}^* A ∗
7.3.4 SVD and the 2-norm¶ The SVD is intimately connected to the 2-norm, as the following theorem describes.
The conclusion (7.3.14) m = n m=n m = n A \mathbf{A} A SVD is the usual means for computing the 2-norm and condition number of a matrix.
Example 7.3.4 (
SVD properties)
Example 7.3.4 We verify some of the fundamental SVD properties using standard Julia functions from LinearAlgebra.
A = [i^j for i = 1:5, j = 0:3]5×4 Matrix{Int64}:
1 1 1 1
1 2 4 8
1 3 9 27
1 4 16 64
1 5 25 125
To get only the singular values, use svdvals.
4-element Vector{Float64}:
146.69715365883005
5.738569780953702
0.9998486640841027
0.11928082685241923
Here is verification of the connections between the singular values, norm, and condition number.
@show opnorm(A, 2);
@show σ[1];opnorm(A, 2) = 146.69715365883005
σ[1] = 146.69715365883005
@show cond(A, 2);
@show σ[1] / σ[end];cond(A, 2) = 1229.846887633767
σ[1] / σ[end] = 1229.846887633767
To get singular vectors as well, use svd. The thin form of the factorization is the default.
U, σ, V = svd(A);
@show size(U);
@show size(V);We verify the orthogonality of the singular vectors as follows:
@show opnorm(U' * U - I);
@show opnorm(V' * V - I);opnorm(U' * U - I) = 4.947764483928556e-16
opnorm(V' * V - I) = Example 7.3.4 We verify some of the fundamental SVD properties using the built-in svd function.
A = vander(1:5);
A = A(:, 1:4)[U, S, V] = svd(A);
disp(sprintf("U is %d by %d. S is %d by %d. V is %d by %d.\n", size(U), size(S), size(V)))U is 5 by 5. S is 5 by 4. V is 4 by 4.
We verify the orthogonality of the singular vectors as follows:
norm(U' * U - eye(5))
norm(V' * V - eye(4))Here is verification of the connections between the singular values, norm, and condition number.
s = diag(S);
norm_A = norm(A)
sigma_max = s(1)cond_A = cond(A)
sigma_ratio = s(1) / s(end)Example 7.3.4 We verify some of the fundamental SVD properties using standard functions from numpy.linalg.
A = array([[(i + 1.0) ** j for j in range(4)] for i in range(5)])
set_printoptions(precision=4)
print(A)[[ 1. 1. 1. 1.]
[ 1. 2. 4. 8.]
[ 1. 3. 9. 27.]
[ 1. 4. 16. 64.]
[ 1. 5. 25. 125.]]
The factorization is obtained using svd from numpy.linalg.
from numpy.linalg import svd
U, sigma, Vh = svd(A)
print("singular values:")
print(sigma)singular values:
[1.4670e+02 5.7386e+00 9.9985e-01 1.1928e-01]
By default, the full factorization type is returned. This can be a memory hog if one of the dimensions of A \mathbf{A} A
print("size of U:", U.shape)
print("size of V:", Vh.T.shape)size of U: (5, 5)
size of V: (4, 4)
Both U \mathbf{U} U V \mathbf{V} V V ∗ \mathbf{V}^* V ∗ V \mathbf{V} V
print(f"should be near zero: {norm(U.T @ U - eye(5), 2):.2e}")
print(f"should be near zero: {norm(Vh @ Vh.T - eye(4), 2):.2e}")should be near zero: 5.41e-16
should be near zero: 7.48e-16
Next we test that we have the factorization promised by the SVD , using diagsvd to construct a rectangular diagonal matrix.
from scipy.linalg import diagsvd
S = diagsvd(sigma, 5, 4)
print(f"should be near zero: {norm(A - U @ S @ Vh, 2):.2e}")should be near zero: 6.55e-14
Here is verification of the connections between the singular values, norm, and condition number.
from numpy.linalg import cond
print("largest singular value:", sigma[0])
print("2-norm of the matrix: ", norm(A, 2))
print("singular value ratio:", sigma[0] / sigma[-1])
print("2-norm condition no.:", cond(A, 2))largest singular value: 146.69715365883002
2-norm of the matrix: 146.69715365883005
singular value ratio: 1229.8468876337672
2-norm condition no.: 1229.846887633767
For matrices that are much taller than they are wide, the thin SVD form is more memory-efficient, because U \mathbf{U} U
A = random.randn(1000, 10)
U, sigma, Vh = svd(A, full_matrices=False)
print("size of U:", U.shape)
print("size of V:", Vh.shape)size of U: (1000, 10)
size of V: (10, 10)
7.3.5 Exercises ¶ ✍ Each factorization below is algebraically correct. The notation I n \mathbf{I}_n I n n × n n\times n n × n SVD . If it is, write down σ 1 \sigma_1 σ 1 u 1 \mathbf{u}_1 u 1 v 1 \mathbf{v}_1 v 1
(a)
[ 0 0 0 − 1 ] = [ 0 1 1 0 ] [ 1 0 0 0 ] [ 0 1 − 1 0 ] \begin{bmatrix}
0 & 0 \\ 0 & -1
\end{bmatrix} = \begin{bmatrix}
0 & 1 \\ 1 & 0
\end{bmatrix} \begin{bmatrix}
1 & 0 \\ 0 & 0
\end{bmatrix} \begin{bmatrix}
0 & 1 \\ -1 & 0
\end{bmatrix} [ 0 0 0 − 1 ] = [ 0 1 1 0 ] [ 1 0 0 0 ] [ 0 − 1 1 0 ]
(b)
[ 0 0 0 − 1 ] = I 2 [ 0 0 0 − 1 ] I 2 \begin{bmatrix}
0 & 0 \\ 0 & -1
\end{bmatrix} =
\mathbf{I}_2 \begin{bmatrix}
0 & 0 \\ 0 & -1
\end{bmatrix}
\mathbf{I}_2 [ 0 0 0 − 1 ] = I 2 [ 0 0 0 − 1 ] I 2
(c)
[ 1 0 0 2 1 0 ] = [ α 0 − α 0 1 0 α 0 − α ] [ 2 0 0 2 0 0 ] [ 0 1 1 0 ] , α = 1 / 2 \begin{bmatrix}
1 & 0\\ 0 & \sqrt{2}\\ 1 & 0
\end{bmatrix} = \begin{bmatrix}
\alpha & 0 & -\alpha \\ 0 & 1 & 0 \\ \alpha & 0 & -\alpha
\end{bmatrix} \begin{bmatrix}
\sqrt{2} & 0 \\ 0 & \sqrt{2} \\ 0 & 0
\end{bmatrix} \begin{bmatrix}
0 & 1 \\ 1 & 0
\end{bmatrix}, \quad \alpha=1/\sqrt{2} ⎣ ⎡ 1 0 1 0 2 0 ⎦ ⎤ = ⎣ ⎡ α 0 α 0 1 0 − α 0 − α ⎦ ⎤ ⎣ ⎡ 2 0 0 0 2 0 ⎦ ⎤ [ 0 1 1 0 ] , α = 1/ 2
(d)
[ 2 2 − 1 1 0 0 ] = I 3 [ 2 0 0 2 0 0 ] [ α α − α α ] , α = 1 / 2 \begin{bmatrix}
\sqrt{2} & \sqrt{2}\\ -1 & 1\\ 0 & 0
\end{bmatrix} =
\mathbf{I}_3 \begin{bmatrix}
2 & 0 \\ 0 & \sqrt{2} \\ 0 & 0
\end{bmatrix} \begin{bmatrix}
\alpha & \alpha \\ -\alpha & \alpha
\end{bmatrix}, \quad \alpha=1/\sqrt{2} ⎣ ⎡ 2 − 1 0 2 1 0 ⎦ ⎤ = I 3 ⎣ ⎡ 2 0 0 0 2 0 ⎦ ⎤ [ α − α α α ] , α = 1/ 2
⌨ Let x be a vector of 1000 equally spaced points between 0 and 1. Suppose A n \mathbf{A}_n A n 1000 × n 1000\times n 1000 × n ( i , j ) (i,j) ( i , j ) x i j − 1 x_i^{j-1} x i j − 1 j = 1 , … , n j=1,\ldots,n j = 1 , … , n
(a) Print out the singular values of A 1 \mathbf{A}_1 A 1 A 2 \mathbf{A}_2 A 2 A 3 \mathbf{A}_3 A 3
(b) Make a log-linear plot of the singular values of A 40 \mathbf{A}_{40} A 40
(c) Repeat part (b) after converting the elements of x to single precision.
(d) Having seen the plot for part (c), which singular values in part (b) do you suspect may be incorrect?
✍ Let A ∈ R m × n \mathbf{A}\in\mathbb{R}^{m\times n} A ∈ R m × n m > n m>n m > n SVD A = U ^ S ^ V T \mathbf{A}=\hat{\mathbf{U}}\hat{\mathbf{S}}\mathbf{V}^T A = U ^ S ^ V T A A + \mathbf{A}\mathbf{A}^{+} A A + U ^ U ^ T \hat{\mathbf{U}}\hat{\mathbf{U}}^T U ^ U ^ T A + \mathbf{A}^{+} A + pseudoinverse A \mathbf{A} A (3.2.4)
✍ In (3.2.6) κ ( A ) = ∥ A ∥ ⋅ ∥ A + ∥ \kappa(\mathbf{A})=\|\mathbf{A}\|\cdot \|\mathbf{A}^{+}\| κ ( A ) = ∥ A ∥ ⋅ ∥ A + ∥ κ ( A ∗ A ) = κ ( A ) 2 \kappa(\mathbf{A}^*\mathbf{A})=\kappa(\mathbf{A})^2 κ ( A ∗ A ) = κ ( A ) 2 SVD .
✍ In this problem you will see how (7.3.14)
(a) Use the technique of Lagrange multipliers to show that among vectors that satisfy ∥ x ∥ 2 2 = 1 \|\mathbf{x}\|_2^2=1 ∥ x ∥ 2 2 = 1 ∥ A x ∥ 2 2 \|\mathbf{A}\mathbf{x}\|_2^2 ∥ Ax ∥ 2 2 A T A \mathbf{A}^T\mathbf{A} A T A B \mathbf{B} B x T B x \mathbf{x}^T\mathbf{B}\mathbf{x} x T Bx x \mathbf{x} x 2 B x 2\mathbf{B}\mathbf{x} 2 Bx
(b) Use the result of part (a) to prove (7.3.14)
✍ Suppose A ∈ R n × n \mathbf{A}\in\mathbb{R}^{n \times n} A ∈ R n × n C \mathbf{C} C (7.3.7)
(a) Suppose that v = [ x y ] \mathbf{v}=\begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix} v = [ x y ] C v = λ v \mathbf{C}\mathbf{v} = \lambda \mathbf{v} Cv = λ v x \mathbf{x} x y \mathbf{y} y
(b) By applying some substitutions, rewrite the equations from part (a) as one in which x \mathbf{x} x y \mathbf{y} y
(c) Substitute the SVD A = U S V T \mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^T A = US V T λ 2 = σ k 2 \lambda^2=\sigma_k^2 λ 2 = σ k 2 σ k \sigma_k σ k
(d) As a more advanced variation, modify the argument to show that λ = 0 \lambda=0 λ = 0 A \mathbf{A} A