Solved – Singular value decomposition procedure in R

rsvd

I'm trying to follow Prof Strang exercise
but I have a problem with the signs of the resultant matrix.

He also had a problem with the signs during the exercise, so I have no way to find what is the error.

Can somebody point me what I'm doing wrong?

Thanks in advance,

Diego

These are my results: 

> A
     [,1] [,2]
[1,]    4    4
[2,]   -3    3
> SVD_of_A
     [,1] [,2]
[1,]   -4    4
[2,]   -3   -3
> U
     [,1] [,2]
[1,]   -1    0
[2,]    0   -1
> sigma
     [,1]     [,2]
[1,] 5.656854 0.000000
[2,] 0.000000 4.242641
> V
      [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

These are the results of the svd function

> svd(A)
$d
[1] 5.656854 4.242641

$u
     [,1] [,2]
[1,]   -1    0
[2,]    0    1

$v
       [,1]       [,2]
[1,] -0.7071068 -0.7071068
[2,] -0.7071068  0.7071068

This is the code I'm using for that exercise:

rm(list=ls())

#------------------------------------------------------------------
# Find the Singular Value Decomposition of the matrix A
#------------------------------------------------------------------

A = matrix(c(4,-3,4,3),2)

#------------------------------------------------------------------
# We want to obtain:
# SVD = U*sigma*t(V)
#------------------------------------------------------------------

#------------------------------------------------------------------
# First we need to obtain the transpose of A
# t = Given a matrix or data.frame x, t returns the transpose of x.
#------------------------------------------------------------------

At= t(A)

#------------------------------------------------------------------
# Multiply the 2 matrices, At and A to obtain AtA
# %*% = Matrix Multiplication
#------------------------------------------------------------------

AtA=At%*%A

#------------------------------------------------------------------
# Now we need to obtain the spectral decomposition of the matrix AtA
# eigen = Computes eigenvalues and eigenvectors of real 
#         (double, integer, logical) or complex matrices.
#------------------------------------------------------------------

eig1 = eigen(AtA)

#------------------------------------------------------------------
# With the eigenvalues of AtA we can build sigma.
# I did this manually here for the purpose of the procedure.
# However note that sigma is equal to sqrt(AAt)
#------------------------------------------------------------------

sigma=matrix(c(sqrt(eig1$values[1]),0,0,sqrt(eig1$values[2])),2)

#------------------------------------------------------------------
# The eigenvectors of AtA are the values of V
#------------------------------------------------------------------

V=as.matrix(eig1$vectors)

#------------------------------------------------------------------
# Next we need to compute A*At
#------------------------------------------------------------------

AAt=A%*%At

#------------------------------------------------------------------
# Now we need to obtain the spectral decomposition of the matrix AAt
#------------------------------------------------------------------

eig2=eigen(AAt)

#------------------------------------------------------------------
# U is equal to the eigen vectors of AAt
#------------------------------------------------------------------

U=as.matrix(eig2$vectors)

#------------------------------------------------------------------
# Finally the singular value decomposition of A is equal to U*sigma*V
#------------------------------------------------------------------

SVD_of_A=U%*%sigma%*%V

#------------------------------------------------------------------
# Results
#------------------------------------------------------------------

A
SVD_of_A
U
sigma
V

PS: I've also tried V transpose, but the results are still wrong…

Best Answer

I did the exact same thing with Matlab/Octave and got the right result.

A = [4 4;-3 3]

At= transpose(A)
AtA = At * A
AAt = A * At

[vec1, val1] = eig(AtA)
V = vec1

[vec2, val2] = eig(AAt)
U = vec2

sigma=sqrt(val1)

SVD_of_A = U * sigma * transpose(V)

Result:

SVD_of_A =

   4.0000   4.0000
  -3.0000   3.0000

The error is that $eigen(AAt)$ should produce eigenvector equals to have the right result to:

      [,1] [,2]
[1,]    1    0
[2,]    0    1

See comment of Cardinal

Related Question