A grayscale image is just a table of numbers, so we can treat it as a matrix and apply the SVD to it. The SVD sorts the content of the image by importance. A few large singular values carry the overall shapes, and many small ones carry fine detail and noise. If we keep only the first n singular values and their vectors, we get an approximation that needs far fewer numbers than the original. I'll go through the decomposition, the truncation, the two results for n = 30 and n = 50, and the storage count. The Matlab code that produced the plots is at the end.

• How is the image turned into a matrix and decomposed ?
• How is the image rebuilt from only n singular values ?
• What do the results for n = 30 and n = 50 show ?
• How much storage does the truncation save ?
• What does the Matlab code do ?

How is the image turned into a matrix and decomposed ?

Before any compression, we need the image in a form that the SVD can take. A grayscale image already has that form. Each pixel is one number between 0 and 255, and the pixels sit in rows and columns like the entries of a matrix.

In this application, I assign the whole pixel data of an black-and-white image to the matrix M

M = imgBW; % 511 x 513 image, 511 x 513 data

Then I did SVD as follows.

Figure 1. The SVD of the image matrix. M is split into two orthonormal bases, U and V, and one diagonal matrix of weights, S.

• U is a 511 x 511 orthogonal matrix : its columns are the left singular vectors. Each one is a vertical pattern, with one value per pixel row.
• S is a 511 x 513 matrix with the singular values on its diagonal : there are min(511, 513) = 511 of them, and they are real, non-negative and sorted from largest to smallest.
• V is a 513 x 513 orthogonal matrix : its columns are the right singular vectors. Each one is a horizontal pattern, with one value per pixel column.
• V^H becomes V^T for an image : the pixel values are real numbers, so the conjugate transpose and the plain transpose are the same thing here.

The decomposition itself does not compress anything. U, S and V together hold more numbers than M does. The saving comes from the ordering of S, which the next section uses.

How is the image rebuilt from only n singular values ?

The full product U S V^T gives back M exactly. The question is how much of it we can drop and still see the same picture. Because S is sorted, the natural choice is to keep the first n singular values and drop the rest.

Then take out only n colums from U and V matrix and take out (n x n) submatrix from S as shown below.

Uc = U(:,1:n);

Sc = S(1:n,1:n);

Vc = V(:,1:n);

Then reconstruct the image from these submatrix as follows.

CompressedImage = Uc * Sc * Vc'   (Note : '*' indicate the matrix multiplication(Inner Product))

Note the transpose on Vc. Uc is 511 x n, Sc is n x n and Vc is 513 x n. So Uc * Sc is 511 x n, and it can only be multiplied by the n x 513 matrix Vc'. The result is 511 x 513, the same size as the original image. The code in List 1 writes it the same way, Uc*Sc*Vc'.

There is another way to read the same product. Write u_i for column i of U, v_i for column i of V and σ_i for the i-th singular value. Then the compressed image M_n is a sum of n layers. Each layer σ_iu_iv_i^T is a rank-one matrix, an outer product of one column pattern and one row pattern. The first layer is the strongest, and every later layer adds a smaller correction.

M_n = Uc * Sc * Vc'

    = sigma_1 u_1 v_1^T + sigma_2 u_2 v_2^T + ... + sigma_n u_n v_n^T

This truncation is also the best possible approximation of its size. The Eckart-Young theorem says that no matrix of rank n is closer to M than M_n, measured either by the Frobenius norm or by the 2-norm. The remaining error is set by the dropped singular values alone. In the Frobenius norm it is the square root of σ_n+1² + ... + σ₅₁₁², and in the 2-norm it is σ_n+1.

• The truncated image has rank n : it is a sum of n rank-one matrices, so it can hold at most n independent row patterns.
• The error depends only on what is dropped : a fast drop in the singular values means a small error for a small n.
• Uc * Sc * Vc needs the transpose on Vc : without it the sizes 511 x n and 513 x n do not fit together, so the product is not defined.

What do the results for n = 30 and n = 50 show ?

Now let's look at what the truncation does to a real picture. Both results below use the same image and the same singular values. Only n changes, so any difference between the two compressed images comes from the 20 extra layers.

Following result shows you the Original Image and all the diagonal values of S matrix and CompressedImage. This is the case when n = 30; In this case, Uc is (511 x 30) matrix, Sc is (30 x 30) matrix, Vc is (513 x 30) matrix.

Figure 2. n = 30. The singular values fall steeply, so 30 of 511 layers already give a recognizable image.

• Left, Original : the 511 x 513 grayscale image M.
• Middle, Diag(S) : all 511 singular values, with the index on the horizontal axis and the value on the vertical axis in units of 10⁴. The first value is about 6.5 x 10⁴, the next ones are already below 1.5 x 10⁴, and the curve is close to zero long before index 100.
• Right, Compressed : the rank-30 image. The face, the hat and the background shapes are clear. Fine texture, such as the hair and the band of the hat, is smoothed. Faint horizontal and vertical stripes appear across the image, because every layer is built from one row pattern and one column pattern.

The first singular value is much larger than all the others for a simple reason. Every pixel of a photo is positive, so the image has a large average brightness. The first layer σ₁u₁v₁^T mostly carries that average level and its slow changes across the picture. The later layers add edges and texture, and each of them carries much less energy.

Following result shows you the Original Image and all the diagonal values of S matrix and CompressedImage. This is the case when n = 50; In this case, Uc is (511 x 50) matrix, Sc is (50 x 50) matrix, Vc is (513 x 50) matrix.

Figure 3. n = 50. The 20 extra layers add fine detail. The Diag(S) plot is the same as in Figure 2, because n does not change the decomposition.

• Middle, Diag(S) : this plot is identical to the one in Figure 2. The singular values belong to the image, and n only decides how many of them are used.
• Right, Compressed : the rank-50 image is closer to the original. The eyes, the hair and the band of the hat are sharper, and the stripes are weaker.
• Diminishing return : layers 31 to 50 have small singular values. So they change the error less than the first 30 layers did, even though they cost the same storage per layer.

How much storage does the truncation save ?

A compressed image is only useful if it takes less space than the original. So let's count the numbers that must be stored. The original needs 511 x 513 = 262143 pixel values. The truncated form needs Uc, Vc and the n singular values on the diagonal of Sc.

That makes n x 511 + n x 513 + n = n x 1025 numbers. The table below applies this count to several values of n. The last column compares it with the 262143 values of the original.

So the n = 30 image in Figure 2 uses about one eighth of the original numbers, and the n = 50 image in Figure 3 uses about one fifth. The break-even point is at n = 262143 / 1025, which is about 255.7. From n = 256 on, the truncated form is larger than the image itself.

• The cost grows linearly with n : each extra layer costs one column of U, one column of V and one singular value, which is 1025 numbers here.
• The count is not the file size : the stored numbers are real values, while a pixel is one byte. A practical codec such as JPEG also uses quantization and entropy coding, which this example does not.
• The example shows the idea of low-rank approximation : the same idea appears in PCA, in noise reduction and in MIMO channel analysis, where a few strong singular values carry most of the energy.

What does the Matlab code do ?

The listing below produced Figure 2. For Figure 3, only the line n = 30 changes to n = 50. The code is short, but a few details in it need attention before you reuse it.

< List 1 >

clear all;

img=imread('Lena.png');
imgBW = rgb2gray(img);

M = imgBW;

[U,S,V]=svd(double(M));
Sd = diag(S);

n = 30;
Uc = U(:,1:n);
Sc = S(1:n,1:n);
Vc = V(:,1:n);
CompressedImag = Uc*Sc*Vc';

subplot(1,3,1);
imshow(M);
title('Original');

subplot(1,3,2);
plot(Sd,'ro','MarkerFaceColor',[1 0 0],'MarkerSize',2);
hold on;
plot(Sd,'k-');
hold on;
xlim([0 length(Sd)]);
title('Diag(S)');

subplot(1,3,3);
imshow(uint8(CompressedImag));
title('Compressed');

• rgb2gray gives the matrix M : it converts the color image to one uint8 value per pixel. The later calls are imshow(M), which shows the original, and svd(double(M)), which needs floating point input.
• svd returns the full matrices : U is 511 x 511, S is 511 x 513 and V is 513 x 513. For a large image, svd(double(M),'econ') or svds with n values saves memory and time.
• Vc' is the conjugate transpose : for real data it is the same as the plain transpose Vc.'.
• uint8 rounds and clips the result : the rank-n image can have values slightly below 0 or above 255, and uint8 saturates them to 0 and 255 before imshow displays it.
• The variable name CompressedImag is used consistently : it is spelled without the final e in all three places, so the code runs as written.

See the SVD page for the decomposition itself, and the Rank page for the meaning of a rank-n matrix.