Math and Art - A Visual Representation of Singular Value Decomposition in Images

We're approaching the end of the fall 2022 semester at the university and by this time I'm ready for the change in scenery. It's dark outside and I'm ready for a good couple of weeks of hibernation. I don't get the luxury of relaxing yet though. It's my first semester back as a master's student, which means I need to make sure I finish strong to make a solid foundation for my GPA. Additionally, I was lucky enough to teach a class of undergraduates Linear Algebra, so I can't slack off for my student's sake.

It's one of my slower days, but as I was designing my student's quiz for singular value decomposition (SVD), I began to wonder just how much information I can cut from an image before it becomes noticeable to the average person. So of course, I open up MATLAB to see if I can manipulate a simple image from my device.

Importing Data to MATLAB

I started by importing an image directly into the console. This can be done with an imread('some-file.jpg'), though the command immediately starts spewing the file contents to the terminal, even if you assign the results to a variable. MATLAB does this a lot, but adding a semicolon (;) to the end of the statement will suppress most output. Upon inspection, the file contents from a JPEG were easily broken up into a tensor. (Vectors are 1D arrays, matrices are 2D arrays, tensors are 3D arrays... I'm not sure if anything comes after...) One dimension was pixel height, another for pixel length, and the last for color channel.

There were two weird things here, the first was that most of linear algebra (at least the amount I have been teaching) doesn't use tensors. I really only paused to ponder this for a moment before just blazing ahead, deciding it was a future me problem. The other is that imread gives back uint8 type numbers, which MATLAB was eager to tell me doesn't work with most of the matrix functions. The numbers read in were fairly small, so converting the whole tensor with single(...) seemed to work just fine. With that fixed, I went about trying to decompose the image with svd(...).

SVD

If you're unfamiliar with SVD, it's a way to break down (and "diagonalize") any matrix into a set of eigenvectors and eigenvalues. Normally you can only do diagonalization like this if the matrix is square because the process of finding eigenvalues involves using the determinant. Determinants don't exist if the matrix isn't square.

If I've lost you at this point, it may be worth checking out 3Blue1Brown's video on determinants. He's also got one on eigenstuff if that is a new concept for you as well.

Square Things Up!

So if a matrix isn't square, it doesn't have a determinant and won't have eigenvalues or eigenvectors, much less a diagonalized form. The way SVD gets around this is instead finding a way to make square matrices out of just the original, non-square matrix. If, for example, we were working with a matrix $A$, we know we will always get a square matrix if we compute either $AA^T$ or $A^TA$. SVD uses both of these to find "singular values" instead of eigenvalues and "singular vectors" instead of eigenvectors.

More specifically, the matrix $U$ from the equation above is made from the left singular vectors, which are the orthonormal eigenvectors of $AA^T$. The matrix $V$ is then made from the right singular vectors which come from the orthonormal set of eigenvectors of $A^TA$. Both the $U$ and $V$ matrices basically contain a normalized set of information from the original matrix that describes where each of the values in $A$ is the highest and lowest. It's $\Sigma$ that holds the information for the real magnitudes.

$\Sigma$ is a matrix that has the "singular values" along its main diagonal. These values are simply the square roots of the eigenvalues found for both $AA^T$ and $A^TA$. Yes, they will be the same every time, with the one caveat that when $A$ is not square, the larger of the two square matrices will just have extra zeros in the set of eigenvalues.

What does this have to do with images?

Well, when we're lucky (or have a really big set of data that isn't very diverse) SVD can allow us to reduce the amount of space we use to save information.

Below is a simple example of a full SVD.

As you can see, this is not very efficient for matrices of full rank. Even its reduced form isn't very efficient.

But small sets of data (and small matrices) don't compress well because there's not a lot of information that is repeated. Compression methods can actually make things worse if you aren't paying attention enough to check that the output is actually smaller. But lucky for us, image data can often have a lot of redundant information. It's unlikely that every pixel in a given image has a color that's wholely unique to the other pixels of that image. SVD can take advantage of this.

Approximation

The examples above show two ways to represent the same matrix in its exact form. There are many situations where retaining all the information in a system is important, and in that case, SVD will be more helpful when the original data is not full rank. However, we have a subjective threshold for how much information we can omit from an image and be unbothered by the difference. SVD also gives us a really concrete way to start omitting large percentages of data that are less meaningful in the context of the original data.

This is because SVD requires us to rank the singular values in $\Sigma$ by order of importance. This also means that $U$ and $V$ must follow suit, since their data is tied directly to those singular values. Singular values are considered more important when they have a larger absolute value. Singular values with the most magnitude are placed at the very top-left corner of $\Sigma$ and progressively smaller values are placed along the main diagonal behind them, ending with any zeros that were found along the way.

There are direct statistical applications to the eigenstuff and singular stuff that come from the SVD process, and I may expand on this later, but for now it's sufficient to know that there is less "information" stored in the smaller singular values and corresponding singular vectors. You might think of it as those later values as only adjusting one pixel by a very small color value in comparison to the larger values being responsible for covering a canvas with a complex mixture of a couple of initial pigments. As you add more and more of the singular values and corresponding vectors, it's like taking finer and finer brushes to a painting. At some point, we can simply choose to stop painting.

Image Reconstruction

Understanding this background, I had decided to ask my students on their quiz a short-answer type question revolving around the approximate SVD process. Nothing much, just making sure they understand the reason for ordering the SVD matrices the way we do. (The alternative was actually making them compute an SVD decomposition, and I'm a fan of maintaining a positive relationship with students.) This is what got me curious and in a squirrel moment I decided I was going to make MATLAB paint me some pretty pictures.

As I mentioned above, importing the data from a JPEG left me with a tensor, which is not the kind of math I had been teaching my students, nor was it any kind of math I had attempted before. Thankfully, MATLAB saw the madness I was attempting and gave me pagesvd(...) to solve the problem that I had three RGB matrix layers that I wanted to decompose. This basically let me treat each color channel as its own matrix and kept the math tidy.

I then had the computer compute the reduced SVD by passing the "econ" optional parameter and simply tried to re-combine the tensors using the SDV formula. MATLAB kindly told me again that tensors are not matrices and can't be multiplied as such, suggesting pagemtimes(...). This call conveniently had options that reminded me of CUDA for allowing one or both tensors to be transposed as needed.

The last part was to present the re-constructed image back again, which is done with imshow(...). If, like me, you needed to convert data in order to get the multiplications and SVD functions to work, you'll need to revert to the old data type to get what you expected again. For some reason the single and double data types don't play well with the image library, but once that was sorted out, the data was a perfect replica of the original image. That's not what I was there for, though.

How bad can it get?

I wanted to see what a less-than-perfect copy would look like. With the "econ" option on for the SVD call, I was still given back 1787 singular values in my image for each color channel, which was full rank in my case. I decided to try a fairly low percentage of singular values for my first approximation, 100. I ran the calculation and... the reconstructed image was actually fine! I cut down almost 91% of the original data from the imported image and I got out an image that was maybe a little blurry in some spots and mildly grainy in others. Think more "low-quality camera" and less "Oh, a computer made this". I was really impressed!

Knowing that I could get such a good reconstruction with only 9% of the original data, I decided to see what the extreme would do. I took a single feature from each of the color channels to see what it would build.

I thought this image was quite lovely. Nothing about it is sharp or jarring and I was anticipating much less nuance in the color gradients. This is the image that inspired this post. I love art that's born out of an elegant piece of mathematical theory. It's also not inaccurate to say that this piece is "the most important"! ... Well, it's the most important part of the original image I took it from, but that's good enough for me.

You can head over to my SVD Gallery to see the original image revealed as more features are added to the reconstruction.