Singular Value Decomposition¶

The last chapter was a lot of fun! This one will be a bit less heavy. We will see another way to decompose matrices: the Singular Value Decomposition or SVD. Since the beginning of this series, I emphasized the fact that you can see matrices as linear transformation in space. With the SVD, you decompose a matrix in three other matrices. You can see these new matrices as sub-transformations of the space. Instead of doing the transformation in one movement, we decompose it in three movements. As a bonus, we will apply the SVD to image processing. We will see the effect of SVD on an image of Lucy the goose (it is just a goose named Lucy...) so keep on reading!

The unit circle and its transformation by a matrix.

We saw in 2.7 that the eigendecomposition can be done only for square matrices. The way to go to decompose other types of matrices that can’t be decomposed with eigendecomposition is to use Singular Value Decomposition (SVD).

We will decompose $\boldsymbol{A}$ into 3 matrices (instead of two with eigendecomposition):

The singular value decomposition

The matrices $\boldsymbol{U}$, $\boldsymbol{D}$, and $\boldsymbol{V}$ have the following properties:

• $\boldsymbol{U}$ and $\boldsymbol{V}$ are orthogonal matrices ($\boldsymbol{U}^\text{T}=\boldsymbol{U}^{-1}$ and $\boldsymbol{V}^\text{T}=\boldsymbol{V}^{-1}$; see previous sections for more details about orthogonal matrices)

• $\boldsymbol{D}$ is a diagonal matrix (all 0 except the diagonal ; see previous sections). However $\boldsymbol{D}$ is not necessarily square.

The columns of $\boldsymbol{U}$ are called the left-singular vectors of $\boldsymbol{A}$ while the columns of $\boldsymbol{V}$ are the right-singular vectors of $\boldsymbol{A}$. The values along the diagonal of $\boldsymbol{D}$ are the singular values of $\boldsymbol{A}$.

Here are the dimensions of the factorization:

The dimensions of the singular value decomposition

The diagonal matrix of singular values is not square but have the shape of $\boldsymbol{A}$. Look at the example provided in the Numpy doc to see that they create a matrix of zeros with the same shape as $\boldsymbol{A}$ and fill it with the singular values:

smat = np.zeros((9, 6), dtype=complex)
smat[:6, :6] = np.diag(s)

Intuition¶

I think that the intuition behind the singular value decomposition needs some explanations about the idea of matrix transformation. For that reason, here are several examples showing how the space can be transformed by 2D square matrices. Hopefully, this will lead to a better understanding of this statement: $\boldsymbol{A}$ is a matrix that can be seen as a linear transformation. This transformation can be decomposed in three sub-transformations: 1. rotation, 2. re-scaling, 3. rotation. These three steps correspond to the three matrices $\boldsymbol{U}$, $\boldsymbol{D}$, and $\boldsymbol{V}$.

$\boldsymbol{A}$ is a matrix that can be seen as a linear transformation. This transformation can be decomposed in three sub-transformations: 1. rotation, 2. re-scaling, 3. rotation. These three steps correspond to the three matrices $\boldsymbol{U}$, $\boldsymbol{D}$, and $\boldsymbol{V}$.

You can look at this animation from the Wikipedia article on the SVD. If you scroll down the page you will see each step.

Every matrix can be seen as a linear transformation¶

You can see a matrix as a specific linear transformation. When you apply this matrix to a vector or to another matrix you will apply this linear transformation to it.

Example 1.¶

We will modify the vector:

by applying the matrix:

We will have:

We see that applying the matrix:

just doubled each coordinate of our vector. Here are the graphical representation of $\boldsymbol{v}$ and its transformation $\boldsymbol{w}$:

Applying the matrix on the vector multiplied each coordinate by two

You can look at other examples of simple transformations on vectors and unit circle in this video.

Example 2.¶

To represent the linear transformation associated with matrices we can also draw the unit circle and see how a matrix can transform it (see the BONUS in previous section). The unit circle represents the coordinates of every unit vectors (vector of length 1).

The unit circle

It is then possible to apply a matrix to all these unit vectors to see the kind of deformation it will produce.

Again, let’s apply the matrix:

to the unit circle:

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

We can see that the matrix doubled the size of the circle. But in some transformations, the change applied to the $x$ coordinate is different from the change applied to the $y$ coordinate. Let’s see what it means graphically.

Example 3.¶

We will apply the matrix:

to the unit circle:

This gives the following new circle:

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

We can check that with the equations associated with this matrix transformation. Let’s say that the coordinates of the new circle (after transformation) are $x'$ and $y'$. The relation between the old coordinates ($x$, $y$) and the new coordinates ($x'$, $y'$) is:

We also know that the equation of the unit circle is $x^2+y^2=1$ (the norm of the unit vectors is 1). By replacement we end up with:

We can check that this equation corresponds to our transformed circle. Let’s start by drawing the old circle. Its equation is:

So far so good!

Coding tip: You can see the trick to plot a circle here: you create the $x$ variable, then $y$ is defined from $x$. This means that for each $x$, the corresponding $y$ value is calculated (and thus $y$ has the same shape as $x$). Since the result of the square root can be negative or positive (for instance, 4 can be the result of 2² but also of $(-2)^2$) we need to plot both solutions ($y$ and $-y$ in plt.plot). Note also that a lot of values are needed if we want the connection between the two demi-spheres. See also some discussion here.

Now let’s add the circle obtained after matrix transformation. We saw that it is defined with

This shows that our transformation was correct.

Note that these examples used diagonal matrices (all zeros except the diagonal). The general rule is that the transformation associated with diagonal matrices imply only a rescaling of each coordinate without rotation. This is a first element to understand the SVD. Look again at the decomposition

We saw that the matrix $\boldsymbol{D}$ is a diagonal matrix. And we saw also that it corresponds to a rescaling without rotation.

Example 4. rotation matrix¶

Matrices that are not diagonal can produce a rotation (see more details here). Since it is easier to think about angles when we talk about rotation, we will use a matrix of the form

This matrix will rotate our vectors or matrices counterclockwise through an angle $\theta$. Our new coordinates will be

Let’s rotate some vectors through an angle of $\theta = 45^\circ$.

Let’s start with the vector $\boldsymbol{u}$ of coordinates $x=0$ and $y=1$ and the vector $\boldsymbol{v}$ of coordinates $x=1$ and $y=0$. The vectors $\boldsymbol{u'}$ $\boldsymbol{v'}$ are the rotated vectors.

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

First, let’s plot $\boldsymbol{u}$ and $\boldsymbol{v}$.

They are the basis vectors of our space. We will calculate the transformation of these vectors:

We will now plot these new vectors to check that they are well our basis vectors rotated through an angle of $45^\circ$.

Coding tip: the numpy functions sin and cos take input in radians. We can convert our angle from degrees to radians with the function np.radians().

We can also transform a circle. We will take a rescaled circle (the one from the example 3.) to be able to see the effect of the rotation.

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

We can see that the circle has been rotated by an angle of $45^\circ$. We have chosen the length of the vectors from the rescaling weight from example 3 (factor 3 and 2) to match the circle.

Summary¶

I hope that you got how vectors and matrices can be transformed by rotating or scaling matrices. The SVD can be seen as the decomposition of one complex transformation in 3 simpler transformations (a rotation, a scaling and another rotation).

Note that we took only square matrices. The SVD can be done even with non square matrices but it is harder to represent transformation associated with non square matrices. For instance, a 3 by 2 matrix will map a 2D space to a 3D space.

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

The three transformations¶

Now that the link between matrices and linear transformation is clearer we can check that a transformation associated with a matrix can be decomposed with the help of the SVD.

But first let’s create a function that takes a 2D matrix as an input and draw the unit circle transformation when we apply this matrix to it. It will be useful to visualize the transformations.

We can use it to check that the three transformations given by the SVD are equivalent to the transformation done with the original matrix. We will also draw each step of the SVD to see the independant effect of the first rotation, the scaling and the second rotation.

We will use the matrix:

and plot the unit circle and its transformation by $\boldsymbol{A}$:

This is what we get when we apply the matrix $\boldsymbol{A}$ to the unit circle and the basis vectors. We can see that the two base vectors are not necessarily rotated the same way. This is related to the sign of the determinant of the matrix (see 2.11).

Let’s now compute the SVD of $\boldsymbol{A}$:

We can now look at the sub-transformations by looking at the effect of the matrices $\boldsymbol{U}$, $\boldsymbol{D}$ and $\boldsymbol{V}$ in the reverse order. Note that it returns the right singular vector already transposed (see the doc).

Just to be sure, you can compare this last step with the transformation by $\boldsymbol{A}$. Fortunately, you will see that the result is the same:

Singular values interpretation¶

The singular values are ordered by descending order. They correspond to a new set of features (that are a linear combination of the original features) with the first feature explaining most of the variance. For instance from the last example we can visualize these new features. The major axis of the elipse will be the first left singular vector ($u_1$) and its norm will be the first singular value ($\sigma_1$).

They are the major ($\sigma_1u_1$) and minor ($\sigma_2u_2$) axes of the elipse. We can see that the feature corresponding to this major axis is associated with more variance (the range of value on this axis is bigger than the other). See 2.12 for more details about the variance explained.

SVD and eigendecomposition¶

Now that we understand the kind of decomposition done with the SVD, we want to know how the sub-transformations are found.

The matrices $\boldsymbol{U}$, $\boldsymbol{D}$ and $\boldsymbol{V}$ can be found by transforming $\boldsymbol{A}$ in a square matrix and by computing the eigenvectors of this square matrix. The square matrix can be obtain by multiplying the matrix $\boldsymbol{A}$ by its transpose in one way or the other:

• $\boldsymbol{U}$ corresponds to the eigenvectors of $\boldsymbol{AA}^\text{T}$

• $\boldsymbol{V}$ corresponds to the eigenvectors of $\boldsymbol{A^\text{T}A}$

• $\boldsymbol{D}$ corresponds to the eigenvalues $\boldsymbol{AA}^\text{T}$ or $\boldsymbol{A^\text{T}A}$ which are the same.

Let’s take an example of a non square matrix:

The singular value decomposition can be done with the linalg.svd() function from Numpy (note that np.linalg.eig(A) works only on square matrices and will give an error for A).

The left-singular values¶

The left-singular values of $\boldsymbol{A}$ correspond to the eigenvectors of $\boldsymbol{AA}^\text{T}$.

Example 5.¶

Note that the sign difference comes from the fact that eigenvectors are not unique. The linalg functions from Numpy return the normalized eigenvectors. Scaling by -1 doesn’t change their direction or the fact that they are unit vectors.

Left singular vectors of A:

Eigenvectors of AA_transpose:

The right-singular values¶

The right-singular values of $\boldsymbol{A}$ correspond to the eigenvectors of $\boldsymbol{A}^\text{T}\boldsymbol{A}$.

Example 6.¶

Right singular vectors of A:

Eigenvectors of A_transposeA:

The nonzero singular values¶

The nonzero singular values of $\boldsymbol{A}$ are the square roots of the eigenvalues of $\boldsymbol{A}^\text{T}\boldsymbol{A}$ and $\boldsymbol{AA}^\text{T}$.

Example 7.¶

Eigenvalues of A_transposeA:

Eigenvalues of AA_transpose:

Square root of the eigenvalues:

BONUS: Apply the SVD on images¶

In this example, we will use the SVD to extract the more important features from the image. It is nice to see the effect of the SVD on something very visual. The code is inspired/taken from this blog post.

Let’s start by loading an image in python and convert it to a Numpy array. We will convert it to grayscale to have one dimension per pixel. The shape of the matrix corresponds to the dimension of the image filled with intensity values: 1 cell per pixel.

We will see how to test the effect of SVD on Lucy the goose! Let’s start to extract the left singular vectors, the singular values and the right singular vectors:

Let’s check the shapes of our matrices:

Remember that $\boldsymbol{D}$ are the singular values that need to be put into a diagonal matrix. Also, $\boldsymbol{V}$ doesn’t need to be transposed (see above).

The singular vectors and singular values are ordered with the first ones corresponding to the more variance explained. For this reason, using just the first few singular vectors and singular values will provide the reconstruction of the principal elements of the image.

We can reconstruct an image from a certain number of singular values. For instance for 2 singular values we will have:

Another representation of the effect of the matrix: each coordinate of the unit circle was multiplied by two

In this example, we have reconstructed the 669px by 1000px image from two singular values.

It is hard to see Lucy with only two singular values and singular vectors. But we already see something!

We will now draw the reconstruction using different number of singular values.

Whaou! Even with 50 components, the quality of the image is not bad!

Conclusion¶

I like this chapter on the SVD because it uses what we have learned so far in a concrete application. The next chapter on the pseudo-inverse is quite cool as well so keep on reading! We will see how to find a near-solution of a system of equation that minimizes the error and at the end we will see an example that uses the pseudo-inverse to find the best fit line of a set of data points.

References¶

Drawing a circle with Matplotlib¶

• https://stackoverflow.com/questions/32092899/plot-equation-showing-a-circle

Rotation matrix¶

• Rotation matrix)

Basis vectors¶

• Basis (linear algebra)

Linear transformation¶

• Aran Glancy - Linear transformation and matrices

SVD¶

• Singular Value Decomposition - Wikipedia

• Professor-svd

• Intoli - PCA and SVD

Numpy¶

• Numpy SVD doc

Image processing¶

• Frank Cleary - SVD of an image