Mathematica Lab 11

MATH 220E Spring 2025

Singular Value Decomposition (SVD)

In this lab, we will use singular value decomposition (SVD) to compress a digital image.

Singular value decomposition is a generalization of diagonalization that can be applied to any matrix, even those that are not square. A singular value decomposition of a matrix A is a product A=U Σ (V^T), where:

U is an n×n matrix whose columns are orthogonal unit vectors.

Σ (the Greek letter capital sigma) is an n×m matrix with zero entries except on its diagonal. The diagonal starts in the upper-left corner, but if Σ is not a square matrix, its diagonal doesn’t go exactly to the lower-right corner. In this case, Σ could look like the following:

(({{1, 0, 0, 0, 0, 0, 0}, {0, 2, 0, 0, 0, 0, 0}, {0, 0, 3, 0, 0, 0, 0}, {0, 0, 0, 4, 0, 0, 0}}))

V is an m×m matrix whose columns are orthogonal unit vectors.

(Singular value decomposition is described in Section 8.4 of our textbook. You do not need to read that section to do this lab.)

Import a Grayscale Image

Begin by importing a smallish picture. The code below imports a sample picture that you may use.

pic=Import["https://mlwright.org/teaching/math220s25/files/table.jpg"]

If you see a picture after entering the line above, you can move to the next step.

Now, we want to convert the picture into data that is recognizable by Mathematica. The following line will do this.

picData=Import["https://mlwright.org/teaching/math220s25/files/table.jpg","Data"]

We have basically imported a big matrix containing grayscale color values. We can check the dimensions of our matrix.

Dimensions[picData]

The sample image is 420 pixels by 620 pixels. (Remember that a digital image is a grid of tiny colored squares called pixels.) Since each pixel has its own color value stored as a number, the image is stored as 420×620=260400 numbers.

By default, the grayscale color value are between 0 and 255. We need to rescale the data values to make each number between 0 and 1, stored as a decimal value. The following code cell does this, storing the result in the variable grayMatrix.

grayMatrix=picData/255.0;

We can display grayMatrix as an image to confirm that things are working so far.

Image[grayMatrix]

Optional: Use your own image

If you want to use your own picture, you may do so. (Choose a sharable picture that isn’t too big!) The easiest way to do this is to type Import[, then when the tool tip appears, click the down arrow symbol and choose File Browser. Navigate to the picture you want to use. Be sure to specify the “Data” parameter as above, then import the data from your picture.

If your import a color picture, then the Dimensions function will give you three numbers, such as {1200, 800, 3}. This means that you have three matrices, each 1200×800 in size. These three matrices give the RBG (red, green, and blue) color values for your image. We need to convert this color image data into grayscale image data. A good way to do this is to average the three color values for each pixel.

Mathematica provides the Mean function that computes averages. Suppose that your color image data is stored in the variable colorPicData. Then the following code cell will compute the mean of the RGB colors of each pixel, storing the results in myPicData. The Dimensions function will now confirm that you have a single matrix of values.

myPicData=Map[Mean,colorPicData,{2}]
Dimensions[myPicData]

Then divide myPicData by 255.0 and store the result in grayMatrix. To do this, un-comment and run the following code cell. (This code is commented so that people who aren’t using their own image won’t accidentally run it and mess things up in their notebook.)

(* grayMatrix=myPicData/255.0
Image[grayMatrix]*)

If you see a grayscale version of your picture, everything is working and you can move on.

SVG and Image compression

Now we’re going to get the singular value decomposition of the grayscale data and check the dimensions of the resulting matrices U, Σ, and V.

{U,Σ,V}=SingularValueDecomposition[grayMatrix];

Dimensions[U]
Dimensions[Σ]
Dimensions[V]

Did the dimensions work out to be what they should? Our original matrix is the product U Σ (V^T), so let’s make sure the image of the product looks the the original image.

Image[U.Σ.Transpose[V]]

We know that the Σ matrix has entries only on its diagonal, so the number of nonzero numbers on the diagonal of Σ must be the same as the number of pivots. Let’s see what this is:

MatrixRank[Σ]

For my image, out of the 420 singular values, all of them ended up being not zero. How many of yours were not zero?

One of the amazing things about SVD is that you can use a truncated SVD to approximate your original matrix with a matrix of much smaller rank. The original Σ above had rank 420 in my case (that is, of the 620 columns, 420 of them were linearly independent). We’ll see soon that we can use a matrix of rank 50 to reconstruct the image quite well. In other words, we are going to use less than ten percent of the data to reconstruct the image!

First we’ll get the truncated SVD, check the new dimensions, then display some compressed images. The SVD[k] function below computes the SVD of my grayscale matrix, truncated at rank k. This simply means is that singular values and corresponding vectors after the kth will be ignored.

SVD[k_]:=SingularValueDecomposition[grayMatrix,k];

We’ll compute SVD[5] below, then look at the (truncated) sizes of the U, Σ, and V matrices. These are the 1st, 2nd, and 3rd components of the return value from SVD[5], so we can look at them by running SVD[5][[i]] for i from 1 to 3, like this:

Dimensions[SVD[5][[1]]]  (*This should ONLY keep the first five columns of the original U matrix*)

Dimensions[SVD[5][[2]]] (*This should ONLY keep the first five diagonal entries of the original Σ matrix*)

Dimensions[SVD[5][[3]]] (*This should ONLY keep the first five rows entries of the original V matrix*)

MatrixRank[SVD[5][[2]]](*this gives the number of singular values, and is 5 since we truncated at k=5*)

Now that we have the truncated SVD components, we’ll multiply them to get an approximation to the original grayscale matrix. The function below will do the multiplying and then show the compressed image.  We’ll do this first for k=5, and then for some larger values of k.

image[k_]:=Image[SVD[k][[1]].SVD[k][[2]].Transpose[SVD[k][[3]]]];

image[5]

The rank 5 image should be a blur, with only a faint resemblance to the original image. This is because it’s generated from only five of the singular values. In my picture, that is only 5 out of 420 singular values, which is just over 1% of the singular values. However, we can store the rank 5 image with much less data than the original image. Adding up the total number of entries in the rank 5 matrices U, Σ, and V, we have to store only 420×5+5×5+620×5=5225 total numbers, which requires vastly less memory than the original image! (In fact, we can do slightly better, since we don’t need to store the off-diagonal entries of Σ, which are all zero.)

Let’s see what happens when we use more singular values!

image[25]

image[50]

image[120]

At rank 25 you can tell roughly what is in my image. By rank 50 the image is mostly clear. At rank 120 the image appears basically the same as the original!  Even at rank 120, we are using only one-third of the total number of singular values, and we can store the three matrices of the SVD in much less memory than the original image

Animation: increasing the rank

Next we’ll animate the images so that you can see the picture coming into focus as the rank increases. Run the following code, and pause the animation when you think your image basically the same as the original.

Animate[image[k],{k,2,MatrixRank[Σ],2,Appearance->"Labeled"},AnimationRate->5]

What is the rank (the k value displayed by the slider above the image)? It should be much less than the original rank. For example, my original image had 420 singular values, but I think my image is recognizable with just 20 singular values. After about rank 100, I can't see much difference in the clarity of the compressed image compared to the original.

Different applications of SVD use different thresholds for singular values. You might choose only singular values that are larger than 1, or maybe just the value larger than 10, or just the overall largest two values, etc. Let’s look at the singular values associated with your picture in a list:

list=SingularValueList[grayMatrix]

Let’s count the number of singular values that are greater than 1. The following code produces this count. The /; syntax is the Condition operator (it evaluates the pattern to its left to see if it satisfies the condition on its right).

Count[list,x_/;x>1]

I found 121 singular values greater than 1. Here is the image created with these singular values:

image[121]

For my image, I don’t even need 121 singular values to get an image that looks like the original. Maybe I only need the singular values larger than 2? Or even just the ones larger than 5? You may experiment more to find out!

Assignment

Experiment a bit and find an integer value n so that you simply need the singular values that are larger than n in order to recognize your picture.

Turn in a single document (not this notebook) that contains:

Your original grayscale picture and its rank.

A compressed but recognizable version of your picture and its rank.

A paragraph describing the experiments you did in order to find out whether you had to keep all singular values greater than 1, or greater than 2, etc., in order to obtain a clear picture.

You can simply copy and paste images into a new Mathematica file, Word document, or Google doc. Add written explanations, save your document as a PDF, and upload it to the Mathematica Lab 11 assignment on Moodle.