Singular Value Decomposition

Recall that in a previous post on eigendecomposition, An $n\times n$ symmetric matrix can be decomposed into $n$ matrices with the same shape ($n\times n$).

$$\mathbf{A}={\lambda }_1{\mathbf{u}}_1{\mathbf{u}}_1^{\mathrm{\top }}+{\lambda }_2{\mathbf{u}}_2{\mathbf{u}}_2^{\mathrm{\top }}+\dots +{\lambda }_n{\mathbf{u}}_n{\mathbf{u}}_n^{\mathrm{\top }},$$

where ${\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_n$ are eigenvectors of $\mathbf{A}$.

Written compactly,

$\mathbf{A}=\mathbf{P}\mathbf{D}{\mathbf{P}}^{\mathrm{\top }}$,

where $\mathbf{P}$ consists of ${\mathbf{u}}_i$ as its column vectors and $\mathbf{D}$ a diagonal matrix with $\mathbf{A}$’s eigenvalues on its diagonal.

More generally, any $m\times n$ matrix $\mathbf{A}$ can be docomposed into $r$ matrices of the same shape $m\times n$, where $r$ is the rank of $\mathbf{A}$. Why should we want to decompose a matrix? Similar to what decomposing a symmetric matrix does, we can approximate a matrix by the sum of its first $k$ components. And why do we want to approximate a matrix? I will answer this question in the next post. In this post, let’s focus on how to decompose a matrix, symmetric or otherwise.

Let $\mathbf{A}$ be an $m\times n$ matrix and $\mathrm{rank}(\mathbf{A})=r$. It can be shown that the number of the non-zero singular values of $\mathbf{A}$ is also its rank, $r$. Since all of those $r$ singular values are positive, we can label them in descending order as

$${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_n$$

where

$${\sigma }_{r+1}={\sigma }_{r+2}=\cdots ={\sigma }_n=0$$

We know that each signular value ${\sigma }_i$ is the square root of ${\lambda }_i$, which are eigenvalues of ${\mathbf{A}}^{\mathrm{\top }}\mathbf{A}$ and correspond to eigenvectors ${\mathbf{v}}_i$ of ${\mathbf{A}}^{\mathrm{\top }}\mathbf{A}$ in the same order. Now we can write the singular value decomposition of $\mathbf{A}$ as:

$$\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\mathrm{\top }}$$

Next, let’s unpack this equation, starting with the item in the middle, $\mathbf{\Sigma }$. $\mathbf{\Sigma }$ (pronouced “sigma”) is an $m\times n$ diagonal matrix, with ${\sigma }_i$ in its diagonal:

$${\mathbf{\Sigma }}_{m\times n}=\left[\begin{array}{ccccccc}{\sigma }_1& 0& \dots & 0& 0& \dots & 0\\ 0& {\sigma }_2& \dots & 0& 0& \dots & 0\\ ⋮& ⋮& \dots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & {\sigma }_r& 0& \dots & 0\\ 0& 0& \dots & 0& {\sigma }_{r+1}& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 0& 0& \dots & {\sigma }_n\end{array}\right]=\left[\begin{array}{ccccccc}{\sigma }_1& 0& \dots & 0& 0& \dots & 0\\ 0& {\sigma }_2& \dots & 0& 0& \dots & 0\\ ⋮& ⋮& \dots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & {\sigma }_r& 0& \dots & 0\\ 0& 0& \dots & 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 0& 0& \dots & 0\end{array}\right]$$

In practice, to construct $\mathbf{\Sigma }$, we can fill an $r\times r$ diagonal matrix with all the non-zero singular values of $\mathbf{A}$, ${\sigma }_1,{\sigma }_2,{\sigma }_r$. Then pad the rest with $0$s to make it an $m\times n$ matrix.

Next to $\mathbf{\Sigma }$ is $\mathbf{V}$, an $n\times n$ matrix consisting of column vectors ${\mathbf{v}}_i$, eigenvectors of the symmetric matrix ${\mathbf{A}}^{\mathrm{\top }}\mathbf{A}$.

$$\mathbf{V}=\left[\begin{array}{cccc}{\mathbf{v}}_1& {\mathbf{v}}_2& \cdots & {\mathbf{v}}_n\end{array}\right]$$

$\mathbf{V}$ is an orthogonal matrix, or orthonormal matrix, because its columns, ${\mathbf{v}}_i$, are orthogonal and normalized. In addition, a set of orthogonal and normalized vectors, such as the set $⟨{\mathbf{v}}_i⟩$ is called an orthonomal set.

Finally, $\mathbf{U}$ is an $m\times m$ orthogonal matrix. To understand how $\mathbf{U}$ is constructed, consider the following statement:

$⟨\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_r⟩$ is an orthogonal basis that spans $\mathrm{Col}(\mathbf{A})$.

Proof. Because ${\mathbf{v}}_i$ and ${\mathbf{v}}_j$ are orthogonal for $i\ne j$,

$$(\mathbf{A}{\mathbf{v}}_i{)}^{\mathrm{\top }}(\mathbf{A}{\mathbf{v}}_j)={\mathbf{v}}_i^{\mathrm{\top }}{\mathbf{A}}^{\mathrm{\top }}\mathbf{A}{\mathbf{v}}_j={\mathbf{v}}_i^{\mathrm{\top }}({\lambda }_j{\mathbf{v}}_j)=0.$$

Therefore, $\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_n$ are orthogonal to each other.

In addition, $\Vert \mathbf{A}{\mathbf{v}}_i\Vert ={\sigma }_i$, where ${\sigma }_i$ are singular values of $\mathbf{A}$. Recall that $\mathbf{A}{\mathbf{v}}_i\ne \mathbf{0}$ when $1\le i\le r$ and $\mathbf{A}{\mathbf{v}}_i=\mathbf{0}$ for $i>r$. So $\mathbf{A}{\mathbf{v}}_1,\dots ,\mathbf{A}{\mathbf{v}}_r$ are orthogonal and all non-zero, and thus a linearly independent set, all of which are in $\mathrm{col}(A)$.

$⟨\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_r⟩$ also spans $\mathrm{col}(A)$. To see why this is true, take any vector $\mathbf{y}$ in $\mathrm{col}(A)$, $\mathbf{y}=\mathbf{A}\mathbf{x}$, where $\mathbf{x}$ is an $n\times 1$ vector. Given that $⟨{\mathbf{v}}_i,\dots ,{\mathbf{v}}_n⟩$ is a basis for ${\mathbb{R}}^n$, we can write $$\mathbf{x}=c_1{\mathbf{v}}_1+\cdots +c_n{\mathbf{v}}_n,$$ so

$$\begin{array}{rlr}\mathbf{y}=\mathbf{A}\mathbf{x}& =c_1\mathbf{A}{\mathbf{v}}_1+\cdots +c_r\mathbf{A}{\mathbf{v}}_r+\cdots +c_n\mathbf{A}{\mathbf{v}}_n& =c_1\mathbf{A}{\mathbf{v}}_1+\cdots +c_r\mathbf{A}{\mathbf{v}}_r.\end{array}$$

In other words, any vector $\mathbf{y}$ in $\mathrm{col}(A)$ can be writen in terms of $⟨\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_r⟩$. Therefore, $⟨\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_r⟩$ is an orthogonal basis for $\mathrm{col}(A)$.

We can normalize vectors $\mathbf{A}{\mathbf{v}}_i$ ($i=1,\dots ,r$) to obtain an orthonormal basis

by dividing them by their length:

$${\mathbf{u}}_i=\frac{\mathbf{A}{\mathbf{v}}_i}{\Vert \mathbf{A}{\mathbf{v}}_i\Vert }=\frac{\mathbf{A}{\mathbf{v}}_i}{{\sigma }_i},1\le i\le r$$

We now have an orthonormal basis $⟨{\mathbf{u}}_i,\dots ,{\mathbf{u}}_r⟩$, where $r$ is the rank of $\mathbf{A}$. In the Singular Value Decomposition equation, $\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\mathrm{\top }}$, $\mathrm{\Sigma }$ is an $m\times n$ matrix. Therefore, $\mathbf{U}$ needs to be a $m\times m$ matrix. In case $r<m$, we need to add additional ortthonormal vectors $⟨{\mathbf{u}}_{r+1}\dots {\mathbf{u}}_m⟩$ to the set so that they span ${\mathbb{R}}^m$. One method to find these $m-r$ vectors is Gram-Schmidt Process. We will introduce it in another post.

Now we have successfully constructed all three components: $\mathrm{\Sigma }$, $\mathbf{V}$, and $\mathbf{U}$. And we can decompose $\mathbf{A}$ as:

$$\begin{array}{rl}\mathbf{A}& =\left[\begin{array}{cccc}{\mathbf{u}}_1& {\mathbf{u}}_2& \dots & {\mathbf{u}}_m\end{array}\right]\left[\begin{array}{ccccccc}{\sigma }_1& 0& \dots & 0& 0& \dots & 0\\ 0& {\sigma }_2& \dots & 0& 0& \dots & 0\\ ⋮& ⋮& \dots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & {\sigma }_r& 0& \dots & 0\\ 0& 0& \dots & 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 0& 0& \dots & 0\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_1^{\mathrm{\top }}\\ {\mathbf{v}}_2^{\mathrm{\top }}\\ ⋮\\ {\mathbf{v}}_n^{\mathrm{\top }}\end{array}\right]\\ & =\left[\begin{array}{cccc}{\mathbf{u}}_1& {\mathbf{u}}_2& \dots & {\mathbf{u}}_m\end{array}\right]\left[\begin{array}{c}{\sigma }_1{\mathbf{v}}_1^{\mathrm{\top }}\\ {\sigma }_2{\mathbf{v}}_2^{\mathrm{\top }}\\ ⋮\\ {\sigma }_r{\mathbf{v}}_r^{\mathrm{\top }}\\ 0\\ ⋮\\ 0\end{array}\right]\\ & ={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^{\mathrm{\top }}+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^{\mathrm{\top }}+\dots +{\sigma }_r{\mathbf{u}}_r{\mathbf{v}}_r^{\mathrm{\top }}\end{array}$$

Given that ${\mathbf{u}}_i$ is an $m$-dimensional column vector, and ${\mathbf{v}}_i$ is an $n$-dimensional column vector, and that ${\sigma }_i$ is a scalar, each ${\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\mathrm{\top }}$ is an $m\times n$ matrix. Therefore, the singular value decomposition equation decomposes matrix $\mathbf{A}$ into $r$ matrices of the same shape.

Closely inspecting the equation $$\mathbf{A}={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^{\mathrm{\top }}+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^{\mathrm{\top }}+\dots +{\sigma }_r{\mathbf{u}}_r{\mathbf{v}}_r^{\mathrm{\top }},$$ we can gain further insights. Multiply both sides of the equation with an $n\times 1$ vector $\mathbf{x}$ and we get:

$$\mathbf{A}\mathbf{x}={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^{\mathrm{\top }}\mathbf{x}+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^{\mathrm{\top }}\mathbf{x}+\dots +{\sigma }_r{\mathbf{u}}_r{\mathbf{v}}_r^{\mathrm{\top }}\mathbf{x},$$

Earlier, we have also shown that $⟨\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\dots ,\mathbf{A}{\mathbf{v}}_r⟩$ is an orthogonal basis that spans $\mathrm{Col}(\mathbf{A})$. Therefore, $\mathbf{A}\mathbf{x}$ can be written as

$$\mathbf{A}\mathbf{x}=a_1{\mathbf{u}}_1+a_2{\mathbf{u}}_2+\dots +a_r{\mathbf{u}}_r$$ Taken togethers, we will have $$a_i={\sigma }_i{\mathbf{v}}_i^{\mathrm{\top }}\mathbf{x}$$

${\mathbf{v}}_i^{\mathrm{\top }}\mathbf{x}$ gives the scalar projection of $\mathbf{x}$ onto ${\mathbf{v}}_i$, which is then multiplied by the $i$-th singular value ${\sigma }_i$. Thus, each component in the decomposition of $\mathbf{A}\mathbf{x}$ is the result of:

• Project $\mathbf{x}$ onto ${\mathbf{v}}_i$,
• multiply the scalar projection with the singular value ${\sigma }_i$,
• and finally scale ${\mathbf{u}}_i$ with the compound scalar resulted from the previous step.

Recall that in the eigendecomposition equation where $\mathbf{A}$ is a symmetric matrix, $$\mathbf{A}\mathbf{x}={\lambda }_1{\mathbf{u}}_1{\mathbf{u}}_1^{\mathrm{\top }}\mathbf{x}+{\lambda }_2{\mathbf{u}}_2{\mathbf{u}}_2^{\mathrm{\top }}\mathbf{x}+\cdots +{\lambda }_n{\mathbf{u}}_n{\mathbf{u}}_n^{\mathrm{\top }}\mathbf{x}$$ $\mathbf{x}$ is projected onto each eigenvector ${\mathbf{u}}_i$ of $\mathbf{A}$, then scaled by the corresponding eigenvalue ${\lambda }_i$.

In the singular value decomposition, $\mathbf{x}$ is projected onto ${\mathbf{u}}_i$, then scaled by the product of singular value ${\sigma }_i$ and scalar projection of $\mathbf{x}$ onto ${\mathbf{v}}_i$, where ${\mathbf{v}}_i$ are eigenvectors of ${\mathbf{A}}^{\mathrm{\top }}\mathbf{A}$, and ${\mathbf{u}}_i$ are unit vectors along the direction of $\mathbf{A}{\mathbf{v}}_i$.

The R package matlib has a function SVD() that can calculate the singular value decomposition of matrix $\mathbf{A}$. Let’s see an example:

mat_a <- matrix(c(4, 8, 1, 3, 3, -2), nrow = 2)
svd_a <- matlib::SVD(mat_a)
svd_a

$d
[1] 9.492982 3.589331

$U
          [,1]       [,2]
[1,] 0.4121068  0.9111355
[2,] 0.9111355 -0.4121068

$V
            [,1]        [,2]
[1,]  0.94148621  0.09686702
[2,]  0.33135146 -0.09059764
[3,] -0.06172462  0.99116540

Note that the last component, $\mathbf{V}$, is a $3\times 2$ matrix, whereas in the original signular value decomposition equation, $\mathbf{V}$ is expected to be a $3\times 3$ matrix, given that $\mathbf{A}$ is a $2\times 3$ matrix. Did the function matlib::SVD() miss something? Consider the following:

$$\begin{array}{rl}\left[\begin{array}{ccc}4& 1& 3\\ 8& 3& -2\end{array}\right]& =\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\mathrm{\top }}\\ & =\left[\begin{array}{cc}0.412& 0.911\\ 0.911& -0.412\end{array}\right]\left[\begin{array}{ccc}9.493& 0& 0\\ 0& 3.589& 0\end{array}\right]\left[\begin{array}{ccc}0.941& 0.331& -0.062\\ 0.097& -0.091& 0.991\\ ?& ?& ?\end{array}\right]\end{array}$$

I have used ? as a place holder for elements in ${\mathbf{V}}^{\mathrm{\top }}$ that were missing from the previous output. It is easy to observe that what’s replaced with ? would be multiplied by 0 from the last column of $\mathrm{\Sigma }$ and would therefore contribute nothing to the final result. This is why matlib::SVD() parsimoniously printed only the first two columns of $\mathbf{V}$.

In the next and final post of this series, we will walk through an example where singular value decoposition solves a real life problem.