Singular value decomposition
clear all
format short
set(0, 'defaultaxesfontsize', 12)
set(0, 'defaultlinelinewidth', 1.5)
set(0, 'defaultFunctionLinelinewidth', 1.5)
set(0, 'defaultscattermarkerfacecolor', 'flat')
gcf;
set(gcf, 'Position', [0 0 600 350], 'Theme', 'light')
addpath ../FNC_matlab7.3. Singular value decomposition ¶ We now introduce another factorization that is as fundamental as the EVD .
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\real^{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 S \mathbf{S} S S ∗ = S T \mathbf{S}^* = \mathbf{S}^T S ∗ = S T 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 , where we have used the fact that V \mathbf{V} V 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 ⎦ ⎤ . The proof also shows that 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 σ \sigma σ A \mathbf{A} A σ \sigma σ − σ -\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 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.
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)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?
✍ In order to prove (7.3.15) Theorem 7.3.3 ∥ A + ∥ 2 = σ n − 1 \twonorm{\mathbf{A}^+} = \sigma_n^{-1} ∥ A + ∥ 2 = σ n − 1 m ≥ n m\ge n m ≥ n A \mathbf{A} A n n n
✍ Let A ∈ R m × n \mathbf{A}\in\real^{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\real^{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) Suppose σ k \sigma_k σ k A \mathbf{A} A x \mathbf{x} x y \mathbf{y} y u k \mathbf{u}_k u k v k \mathbf{v}_k v k σ k \sigma_k σ k − σ k -\sigma_k − σ k C \mathbf{C} C