49
$\begingroup$

I have read about singular value decomposition (SVD). In almost all textbooks it is mentioned that it factorizes the matrix into three matrices with given specification.

But what is the intuition behind splitting the matrix in such form? PCA and other algorithms for dimensionality reduction are intuitive in the sense that algorithm has nice visualization property but with SVD it is not the case.

$\endgroup$
  • 4
    $\begingroup$ You might want to start from the intuition of eigenvalue-eigenvector decomposition as SVD is an extension of it for all kinds of matrices, instead of just square ones. $\endgroup$ – JohnK Oct 15 '15 at 17:43
  • $\begingroup$ There are plenty of notes on internet and answers here on CV about SVD and its workings. $\endgroup$ – Vladislavs Dovgalecs Oct 15 '15 at 17:53
  • 2
    $\begingroup$ SVD can be thought as a compression/learning algorithm. It is a linear compressor decompressor. A matrix M can be represented by multiplication of SVD. S is the compressor V determines how much error you would like to have (lossy compression) and D is the decompressor. If you keep all diagonal values of V then you have a lossless compressor. If you start throwing away small singular values (zeroing them) then you cannot reconstruct the initial matrix exactly but will still be close. Here the term close is measured with Frobenius norm. $\endgroup$ – Cowboy Trader Oct 15 '15 at 18:32
  • 2
    $\begingroup$ @Cagdas if you do that please carefully define what you're taking "S" "V" and "D" to be mathematically. I've not seen the initials overloaded into the notation itself before (which has the singular values in it, for example?). It seems to be a likely source of confusion, $\endgroup$ – Glen_b Oct 15 '15 at 21:46
  • 3
    $\begingroup$ Do you know how to estimate PCA with SVD? If you do, then can you explain why you feel that something missing in your understanding of SVD? See this $\endgroup$ – Aksakal Oct 28 '15 at 16:52
62
$\begingroup$

Write the SVD of matrix $X$ (real, $n\times p$) as $$ X = U D V^T $$ where $U$ is $n\times p$, $D$ is diagonal $p\times p$ and $V^T$ is $p\times p$. In terms of the columns of the matrices $U$ and $V$ we can write $X=\sum_{i=1}^p d_i u_i v_i^T$. That shows $X$ written as a sum of $p$ rank-1 matrices. What does a rank-1 matrix look like? Let's see: $$ \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \begin{pmatrix} 4 & 5 & 6 \end{pmatrix} = \begin{pmatrix} 4 & 5 & 6 \\ 8 & 10 & 12 \\ 12 & 15 & 18 \end{pmatrix} $$ The rows are proportional, and the columns are proportional.

Think now about $X$ as containing the greyscale values of a black-and-white image, each entry in the matrix representing one pixel. For instance the following picture of a baboon:

image of a baboon

Then read this image into R and get the matrix part of the resulting structure, maybe using the library pixmap.


If you want a step-by-step guide as to how to reproduce the results, you can find the code here.


Calculate the SVD:

baboon.svd  <-  svd(bab) # May take some time

How can we think about this? We get the $512 \times 512$ baboon image represented as a sum of $512$ simple images, with each one only showing vertical and horizontal structure, i.e. it is an image of vertical and horizontal stripes! So, the SVD of the baboon represents the baboon image as a superposition of $512$ simple images, each one only showing horizontal/vertical stripes. Let us calculate a low-rank reconstruction of the image with $1$ and with $20$ components:

baboon.1  <-  sweep(baboon.svd$u[,1,drop=FALSE],2,baboon.svd$d[1],"*") %*%
                   t(baboon.svd$v[,1,drop=FALSE])

baboon.20 <-  sweep(baboon.svd$u[,1:20,drop=FALSE],2,baboon.svd$d[1:20],"*") %*%
                   t(baboon.svd$v[,1:20,drop=FALSE])

resulting in the following two images:

rank one and rank 20 reconstruction of baboon image

On the left we can easily see the vertical/horizontal stripes in the rank-1 image.

Let us finally look at the "residual image", the image reconstructed (as above, code not shown) from the $20$ rank-one images with the lowest singular values. Here it is:

image of residuals from rank 20 baboon reconstruction

Which is quite interesting: we see the parts of the original image that are difficult to represent as superposition of vertical/horizontal lines, mostly diagonal nose hair and some texture, and the eyes!

$\endgroup$
  • 10
    $\begingroup$ I think you meant low-rank reconstruction, not low range. Never mind. This is a very good illustration (+1). That's why it is a linear compressor decompressor. Image is approximated with lines. If you actually perform a similar autoencoder with a neural network with linear activation functions, you will actually see that it also allows lines with any slope not only vertical and horizontal lines, which makes it slightly more powerful than SVD. $\endgroup$ – Cowboy Trader Oct 28 '15 at 14:31
  • $\begingroup$ @kjetil-b-halvorsen Shouldn't the SVD $X = U \Sigma V^*$ for a $n \times p$ matrix $X$ result in $U$ being $n \times n$, $\Sigma$ being $n \times p$ and $V$ being $p \times p$? Typo? $\endgroup$ – Martin Krämer Jan 26 '17 at 17:06
  • 1
    $\begingroup$ See math.stackexchange.com/questions/92171/… for some other examples $\endgroup$ – kjetil b halvorsen May 9 '17 at 10:45
  • $\begingroup$ @kjetil-b-halvorsen I am interested in knowing how would the decription change if I would have used PCA for denosing application. I would appreciate if you could answer my question here stats.stackexchange.com/questions/412123/… $\endgroup$ – Dushyant Kumar Jun 8 at 21:31
  • $\begingroup$ @CowboyTrader interesting observation. My understanding of machine learning/neural network is pretty limited. So, I fail to understand that if one has a single noisy image and nothing else to train on, how would the neural network work? $\endgroup$ – Dushyant Kumar Jun 8 at 21:38
2
$\begingroup$

Let $A$ be a real $m \times n$ matrix. I'll assume that $m \geq n$ for simplicity. It's natural to ask in which direction $v$ does $A$ have the most impact (or the most explosiveness, or the most amplifying power). The answer is \begin{align} \tag{1}v_1 = \,\,& \arg \max_{v \in \mathbb R^n} \quad \| A v \|_2 \\ & \text{subject to } \, \|v\|_2 = 1. \end{align} A natural follow-up question is, after $v_1$, what is the next most explosive direction for $A$? The answer is \begin{align} v_2 = \,\,& \arg \max_{v \in \mathbb R^n} \quad \| A v \|_2 \\ & \text{subject to } \,\langle v_1, v \rangle = 0, \\ & \qquad \qquad \, \, \, \, \|v\|_2 = 1. \end{align} Continuing like this, we obtain an orthonormal basis $v_1, \ldots, v_n$ of $\mathbb R^n$. This special basis of $\mathbb R^n$ tells us the directions that are, in some sense, most important for understanding $A$.

Let $\sigma_i = \|A v_i \|_2$ (so $\sigma_i$ quantifies the explosive power of $A$ in the direction $v_i$). Suppose that unit vectors $u_i$ are defined so that $$ \tag{2} A v_i = \sigma_i u_i \quad \text{for } i = 1, \ldots, n. $$ The equations (2) can be expressed concisely using matrix notation as $$ \tag{3} A V = U \Sigma, $$ where $V$ is the $n \times n$ matrix whose $i$th column is $v_i$, $U$ is the $m \times n$ matrix whose $i$th column is $u_i$, and $\Sigma$ is the $n \times n$ diagonal matrix whose $i$th diagonal entry is $\sigma_i$. The matrix $V$ is orthogonal, so we can multiply both sides of (3) by $V^T$ to obtain $$ A = U \Sigma V^T. $$ It might appear that we have now derived the SVD of $A$ with almost zero effort. None of the steps so far have been difficult. However, a crucial piece of the picture is missing -- we do not yet know that $U$ is orthogonal.

Here is the crucial fact, the missing piece: it turns out that $A v_1$ is orthogonal to $A v_2$: $$ \tag{4} \langle A v_1, A v_2 \rangle = 0. $$ I claim that if this were not true, then $v_1$ would not be optimal for problem (1). Indeed, if (4) were not satisfied, then it would be possible to improve $v_1$ by perturbing it a bit in the direction $v_2$.

Suppose (for a contradiction) that (4) is not satisfied. If $v_1$ is perturbed slightly in the orthogonal direction $v_2$, the norm of $v_1$ does not change (or at least, the change in the norm of $v_1$ is negligible). When I walk on the surface of the earth, my distance from the center of the earth does not change. However, when $v_1$ is perturbed in the direction $v_2$, the vector $A v_1$ is perturbed in the non-orthogonal direction $A v_2$, and so the change in the norm of $A v_1$ is non-negligible. The norm of $A v_1$ can be increased by a non-negligible amount. This means that $v_1$ is not optimal for problem (1), which is a contradiction. I love this argument because: 1) the intuition is very clear; 2) the intuition can be converted directly into a rigorous proof.

A similar argument shows that $A v_3$ is orthogonal to both $A v_1$ and $A v_2$, and so on. The vectors $A v_1, \ldots, A v_n$ are pairwise orthogonal. This means that the unit vectors $u_1, \ldots, u_n$ can be chosen to be pairwise orthogonal, which means the matrix $U$ above is an orthogonal matrix. This completes our discovery of the SVD.


To convert the above intuitive argument into a rigorous proof, we must confront the fact that if $v_1$ is perturbed in the direction $v_2$, the perturbed vector $$ \tilde v_1 = v_1 + \epsilon v_2 $$ is not truly a unit vector. (Its norm is $\sqrt{1 + \epsilon^2}$.) To obtain a rigorous proof, define $$ \bar v_1(\epsilon) = \sqrt{1 - \epsilon^2} v_1 + \epsilon v_2. $$ The vector $\bar v_1(\epsilon)$ is truly a unit vector. But as you can easily show, if (4) is not satisfied, then for sufficiently small values of $\epsilon$ we have $$ f(\epsilon) = \| A \bar v_1(\epsilon) \|_2^2 > \| A v_1 \|_2^2 $$ (assuming that the sign of $\epsilon$ is chosen correctly). To show this, just check that $f'(0) \neq 0$. This means that $v_1$ is not optimal for problem (1), which is a contradiction.

(By the way, I recommend reading Qiaochu Yuan's explanation of the SVD here. In particular, take a look at "Key lemma # 1", which is what we discussed above. As Qiaochu says, key lemma # 1 is "the technical heart of singular value decomposition".)

$\endgroup$
0
$\begingroup$

Dude take an hour of your day and watch this lecture: https://www.youtube.com/watch?v=EokL7E6o1AE

This guy is super straight forward, it's important not to skip any of it because it all comes together in the end. Even if it might seem a little slow at the beginning, he is trying to pin down a critical point, which he does!

I'll sum it up for you, rather than just giving you the three matrices that everyone does (because that was confusing me when I read other descriptions). Where do those matrices come from and why do we set it up like that? The lecture nails it! Every matrix (ever in the history of everness) can be constructed from a base matrix with same dimensions, then rotate it, and stretch it (this is the fundamental theorem of linear algebra). Each of those three matrices people throw around represent an initial matrix (U), A scaling matrix (sigma), and a rotation matrix (V).

The scaling matrix shows you which rotation vectors are dominating, these are called the singular values. The decomposition is solving for U, sigma, and V.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.