Home › Blogs › jpmay's blog

Stupid SVD Tricks Part 1

jpmay's picture
Posted on 2008-03-18 17:40 By jpmay (Maple Rating 2 114)
Categories:

Of all the ways to decompose a numerical (floating point) matrix, my favorite is the singular value decomposition (SVD).  There are a lot of applications of the SVD (see my dissertation for one related to polynomial algebra) but my favorite ones are probably two applications related to image processing.

The first one I want to talk about comes from the cover of James Demmel's book "Applied Numerical Linear Algebra": image compression.  This example gives a really cool intuitive understanding of the Rank of Matrix and is also nice excuse to play with Maple's ImageTools package.

So, the first thing you need a test image. I used the classic image compression benchmark of a Mandrill.



Read this in with:


mandrill:=ImageTools:-Read("4.2.03.tiff");


The result is a 512x512x3 array.  In order to do something with this, we need to make it into a matrix so, call


manmat:=convert(ArrayTools:-Reshape(mandrill, 512*3, 512), Matrix);


Now we can compute a singular value decomposition of the image:


(U, S, V) := LinearAlgebra:-SingularValues(manmat, output = ['U', 'S', 'Vt']);

Now we can zero-out small singular values and multiply things back together to create low-rank approximations of the matrix that are also compressed versions of the image.
Rank 32 will give us 1/8 of the data (64 dimension 512 vectors: 32 rows of U, 32 columns of V, and the 32 corresponding singular values) but still a pretty good image:


rank32approx:=MatrixMatrixMultiply(`.`(U, DiagonalMatrix(S[1..32], 3*512, 512)), V, outputoptions = [order = C_order]);


This reshape it back to an image and display:


Preview((Reshape(rank32approx, 512, 512, 3)));


Taking things down to rank 8, is leaving only 1/32 of the data, but it is amazing how what is left resembles the original image:


rank8approx:=MatrixMatrixMultiply(`.`(U, DiagonalMatrix(S[1..8], 3*512, 512)), V, outputoptions = [order = C_order]);
Preview((Reshape(rank8approx, 512, 512, 3)));


To look at more images in order of descending rank, take a look at my worksheet:
Download 5480_SVD-face-colour-improved.mw
View file details

Next time: eigenfaces

Comments

JacquesC's picture

Natural compression?

Comment on 2008-03-19 21:40 from JacquesC (Maple Rating 4 1655)

Instead of using an arbitrary number (like 32 or 8) for the rank, what about looking through the resulting singular values instead for the most natural point where a 'jump' occurs?  In almost all cases, there is such a jump, sometimes more than one.  These all seem to correspond to natural breakpoints in the amount of information contained in the "rest".  A nice piece of Maple code to find those natural points would be a convenient addition to this fun and informative post.

jpmay's picture

Good Suggestion

Comment on 2008-03-20 00:17 from jpmay (Maple Rating 2 114)

It is actually not true that all matrices have a clear natural gap.  In approximate algebra applications you can run across matrices where the singular values are pretty much evenly spaced.

That said, this image has a very large ("natural") gap between the third and fourth singular values:
# S is the vector of singular values from the worksheet

Srat := [seq(S[i]/S[i+1], i = 1 .. 511)]:
SSrat:=sort(Srat);
evalf[5](SSrat[-8..-1]);
  [1.0875, 1.1044, 1.1639, 1.1741, 1.2214, 1.2647, 1.2871, 4.1747]
member(SSrat[-1], Srat, 'loc'): loc;
  3

And this is the largest gap by a pretty large margin.  The next five biggest gaps are at 4, 6, 7, 5, and 1 respectively. This is pretty interesting.  This suggests you get a lot of information when adding in the first 7 or so singular vector pairs but that it tapers off quite a bit after that.  I think you can kind of see this in rank 8 vs. rank 32 pictures.

I'll post the rank 3 image tomorow for comparison.
 

Tip to older computer (e.g. first gen g4 powerbook) users:  run the worksheet on an 256x256 image instead.

acer's picture

sometimes difficult

Comment on 2008-03-20 02:49 from acer (Maple Rating 4 1107)

I'm not sure whether the "information jumps" to which Jacques alluded are better measured as (the ratio of) absolute deltas or actual entries.

Eg,

plots[pointplot]([seq([i, (S[i]-S[i+1])/(S[i+1]-S[i+2])], i = 1 .. 50)]);

That seems to bring out jump points that match OK with the reciprocal log plot in another comment below.

Like so much else with plots, it's sometimes far easier to see a quality than to measure and detect it with code in a reliable way.

acer

JacquesC's picture

The gap

Comment on 2008-03-20 09:06 from JacquesC (Maple Rating 4 1655)

Note that I was quite careful to say "almost always" when talking about finding gaps.  Yes, one can construct matrices with no gaps.  However, it seems that for naturally occuring phenomenon with real information in them, there is usually a gap [counter-examples here would be very interesting].

So yes, it would be very interesting to see rank-3 vs rank-4, as it seems that that is where a lot of information "finally" comes in, and that going up to 7 gives you the majority of the information, and the rest just fills-in the details.

acer's picture

interesting

Comment on 2008-03-19 23:25 from acer (Maple Rating 4 1107)

It does seem like there are jumps, at elements 2,4,8,15,18, 21,26,28,..

plots[pointplot]([seq([i, 1/abs(ln(S[i]))], i = 1 .. 30)]);

acer

jpmay's picture

More about Error

Comment on 2008-03-20 14:33 from jpmay (Maple Rating 2 114)

If you want to compute the pixel-by-pixel error between the original image and the low-rank "compressed" image you can compute the matrix norm of the differences between the two images:
LinearAlgebra:-MatrixNorm(faceM-face8, Frobenius); # face8 is the rank=4 image
   78.987

A well known theorem says that this norm is equal to the euclidean norm of the vector S[5..512], that is: sqrt(S[5]^2+S[6]^2+..+S[511]^2+S[512]^2)
Thus, it is pretty easy to plot a graph of how the pixel-by-pixel error decreases as you increase the rank of the compressed image:

5480_rankvserror.jpeg

And the curve stays pretty smooth, decreasing to 0 error at rank 512.  This plot confirms the information from the gaps in the singular values: most of  the information about the image is in the first 3 singular vectors, and after rank 7, addtional increases in rank don't buy you much.  But, that rank 3 approximation is still pretty far from the original image:

 

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}