9
$\begingroup$

I am trying to do SVD by hand:

m<-matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U=eigen(m%*%t(m))$vector
V=eigen(t(m)%*%m)$vector
D=sqrt(diag(eigen(m%*%t(m))$values))

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d)

U1%*%D1%*%t(V1)
U%*%D%*%t(V)

But the last line does not return m back. Why? It seems to has something to do with signs of these eigenvectors... Or did I misunderstand the procedure?

$\endgroup$
  • $\begingroup$ See stats.stackexchange.com/search?q=svd+sign. $\endgroup$ – whuber Feb 3 '17 at 22:31
  • $\begingroup$ I am repeatedly told that sign doesn't matter in SVDs...like this $\endgroup$ – failedstatistician Feb 3 '17 at 22:33
  • $\begingroup$ @Amoeba Thank you for clarifying that. I was focusing on the English question rather than the code. Failedstatistician: see what D=diag(c(-1,1,1)*sqrt(eigen(m%*%t(m))$values)) does and bear in mind that the square root (as well as any normalized eigenvector) is defined only up to sign. For more insight, change m to m <- matrix(-2,1,1) and include ,1,1) at the end of each of the calls to diag. This is a $1\times 1$ example that creates the same problem--but it's so simple the nature of the problem will become completely obvious. $\endgroup$ – whuber Feb 3 '17 at 22:40
  • 1
    $\begingroup$ Got it. Thank you! Do you have a general rule of determining the vector c(-1, 1, 1)? Or how the signs of the two decomposition should be linked? $\endgroup$ – failedstatistician Feb 3 '17 at 23:00
  • 1
    $\begingroup$ Note that @whuber's trick with c(-1,1,1) does work, but D defined like that is not giving you singular values. Singular values must all be positive by definition. The question of how to link the signs of U and V is good, and I don't have an answer. Why don't you just do an SVD? :-) $\endgroup$ – amoeba Feb 3 '17 at 23:14
13
$\begingroup$

Analysis of the Problem

The SVD of a matrix is never unique. Let matrix $A$ have dimensions $n\times k$ and let its SVD be

$$A = U D V^\prime$$

for an $n\times p$ matrix $U$ with orthonormal columns, a diagonal $p\times p$ matrix $D$ with non-negative entries, and a $k\times p$ matrix $V$ with orthonormal columns.

Now choose, arbitrarily, any diagonal $p\times p$ matrix $S$ having $\pm 1$s on the diagonal, so that $S^2 = I$ is the $p\times p$ identity $I_p$. Then

$$A = U D V^\prime = U I D I V^\prime = U (S^2) D (S^2) V^\prime = (US) (SDS) (VS)^\prime$$

is also an SVD of $A$ because $$(US)^\prime(US) = S^\prime U^\prime U S = S^\prime I_p S = S^\prime S = S^2 = I_p$$ demonstrates $US$ has orthonormal columns and a similar calculation demonstrates $VS$ has orthonormal columns. Moreover, since $S$ and $D$ are diagonal, they commute, whence $$S D S = DS^2 = D$$ shows $D$ still has non-negative entries.

The method implemented in the code to find an SVD finds a $U$ that diagonalizes $$AA^\prime = (UDV^\prime)(UDV^\prime)^\prime = UDV^\prime V D^\prime U^\prime = UD^2 U^\prime$$ and, similarly, a $V$ that diagonalizes $$A^\prime A = VD^2V^\prime.$$ It proceeds to compute $D$ in terms of the eigenvalues found in $D^2$. The problem is this does not assure a consistent matching of the columns of $U$ with the columns of $V$.

A Solution

Instead, after finding such a $U$ and such a $V$, use them to compute

$$U^\prime A V = U^\prime (U D V^\prime) V = (U^\prime U) D (V^\prime V) = D$$

directly and efficiently. The diagonal values of this $D$ are not necessarily positive. (That is because there is nothing about the process of diagonalizing either $A^\prime A$ or $AA^\prime$ that will guarantee that, since those two processes were carried out separately.) Make them positive by choosing the entries along the diagonal of $S$ to equal the signs of the entries of $D$, so that $SD$ has all positive values. Compensate for this by right-multiplying $U$ by $S$:

$$A = U D V^\prime = (US) (SD) V^\prime.$$

That is an SVD.

Example

Let $n=p=k=1$ with $A=(-2)$. An SVD is

$$(-2) = (1)(2)(-1)$$

with $U=(1)$, $D=(2)$, and $V=(-1)$.

If you diagonalize $A^\prime A = (4)$ you would naturally choose $U=(1)$ and $D=(\sqrt{4})=(2)$. Likewise if you diagonalize $AA^\prime=(4)$ you would choose $V=(1)$. Unfortunately, $$UDV^\prime = (1)(2)(1) = (2) \ne A.$$ Instead, compute $$D=U^\prime A V = (1)^\prime (-2) (1) = (-2).$$ Because this is negative, set $S=(-1)$. This adjusts $U$ to $US = (1)(-1)=(-1)$ and $D$ to $SD = (-1)(-2)=(2)$. You have obtained $$A = (-1)(2)(1),$$ which is one of the two possible SVDs (but not the same as the original!).

Code

Here is modified code. Its output confirms

  1. The method recreates m correctly.
  2. $U$ and $V$ really are still orthonormal.
  3. But the result is not the same SVD returned by svd. (Both are equally valid.)
m <- matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U <- eigen(tcrossprod(m))$vector
V <- eigen(crossprod(m))$vector
D <- diag(zapsmall(diag(t(U) %*% m %*% V)))
s <- diag(sign(diag(D)))  # Find the signs of the eigenvalues
U <- U %*% s              # Adjust the columns of U
D <- s %*% D              # Fix up D.  (D <- abs(D) would be more efficient.)

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d,n,n)

zapsmall(U1 %*% D1 %*% t(V1)) # SVD
zapsmall(U %*% D %*% t(V))    # Hand-rolled SVD
zapsmall(crossprod(U))        # Check that U is orthonormal
zapsmall(tcrossprod(V))       # Check that V' is orthonormal
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1. This is very clear. I would only add that in practice it's enough to compute either U or V and then to obtain another matrix via multiplying with A. This way one performs only one (instead of two) eigendecompositions, and the signs will come out right. $\endgroup$ – amoeba Feb 3 '17 at 23:36
  • 2
    $\begingroup$ @Amoeba That is right: in the spirit of hand-computing an SVD, which evidently is an educational exercise, no attention is paid here to efficiency. $\endgroup$ – whuber Feb 3 '17 at 23:38
  • 2
    $\begingroup$ Thank you for your kind help! I think I understand this issue (finally). $\endgroup$ – failedstatistician Feb 4 '17 at 16:06
  • 3
    $\begingroup$ @Federico Thank you for that reminder. You are quite correct--I have implicitly assumed all eigenvalues are distinct, for indeed that is almost surely going to be the case in statistical applications and one gets out of the habit of considering the ambiguities with "degenerate" eigenspaces. $\endgroup$ – whuber Feb 8 '17 at 16:41
  • 3
    $\begingroup$ You are correct, this is only an edge case, and a tricky one indeed. In some sense, it is another manifestation of the same problem that you outline in your answer, that this method doesn't ensure a "matching" between the columns of $U$ and $V$. Computing the SVD starting from the eigendecompositions is still a great learning example. $\endgroup$ – Federico Poloni Feb 8 '17 at 17:37
5
$\begingroup$

As I outlined in a comment to @whuber's answer, this method to compute the SVD doesn't work for every matrix. The issue is not limited to signs.

The problem is that there may be repeated eigenvalues, and in this case the eigendecomposition of $A'A$ and $AA'$ is not unique and not all choices of $U$ and $V$ can be used to retrieve the diagonal factor of the SVD. For instance, if you take any non-diagonal orthogonal matrix (say, $A=\begin{bmatrix}3/5&4/5\\-4/5&3/5\end{bmatrix}$), then $AA'=A'A=I$. Among all possible choices for the eigenvector matrix of $I$, eigen will return $U=V=I$, thus in this case $U'AV=A$ is not diagonal.

Intuitively, this is another manifestation of the same problem that @whuber outlines, that there has to be a "matching" between the columns of $U$ and $V$, and computing two eigendecompositions separately does not ensure it.

If all the singular values of $A$ are distinct, then the eigendecomposition is unique (up to scaling/signs) and the method works. Remark: it is still not a good idea to use it in production code on a computer with floating point arithmetic, because when you form the products $A'A$ and $AA'$ the computed result may be perturbed by a quantity of the order of $\|A\|^2u$, where $u \approx 2\times 10^{-16}$ is the machine precision. If the magnitudes of the singular values differ greatly (of more than $10^{-8}$, roughly), this is detrimental to the numerical accuracy of the smallest ones.

Computing the SVD from the two eigendecompositions is a great learning example, but in real life applications always use R's svd function to compute the singular value decomposition.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ This comment is good advice. Please note, though, that this thread is not concerned about the proper way to compute SVD (and I believe nobody would argue against your recommendation). The OP implicitly accepts that svd works. Indeed, they use it as a standard against which to compare a hand calculation, whose purpose is to check the understanding, not to replace svd in any way. $\endgroup$ – whuber Feb 8 '17 at 16:43
  • $\begingroup$ @whuber Correct observation; I reworded last paragraph. $\endgroup$ – Federico Poloni Feb 8 '17 at 17:43

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.