Singular Value Decomposition and Pseudo-Inverse

Singular Value Decomposition and Pseudo-Inverse #

Implementing Singular Value Decomposition and Pseudo-Inverse algorithms and comparing them to R’s built in algorithms to test their accuracy.

.rmd knitted with execution results available for download here

Q1#

Compute the SVD of the following matrices#

1
A = matrix(c(c(1,0), c(1,0)), nc=2, nr=2)
2
A

1
##      [,1] [,2]
2
## [1,]    1    1
3
## [2,]    0    0

1
B = matrix(c(c(-2,0), c(0, 0)), nc=2, nr=2)
2
B

1
##      [,1] [,2]
2
## [1,]   -2    0
3
## [2,]    0    0

1
C = matrix(c(c(2,0), c(0,2), c(1,0)), nc=3, nr=2)
2
C

1
##      [,1] [,2] [,3]
2
## [1,]    2    0    1
3
## [2,]    0    2    0

The SVD formula#

The SVD of any matrix ‘M’, size ‘m’ by ‘n’ is:

M = U * E * V^T,

where U is the orthonormal eigenvectors of M * M^T,

and where V is the orthonormal eigenvectors of M^T * M,

and where E is a diagonal matrix of the singular values (square root of eigenvalues) of size ‘m’ by ‘n’

SVD of A#

1
# calculating U
2
AAt = A %*% t(A)
3
eu = eigen(AAt)
4
# for some reason the eigenvectors are negative, I assume this has something to do with the eigen() implementation
5
# it should be noted that below, when using svd(), U has no negative values.
6
# So, since U needs to be the orthonormal eigenvectors of 'AAt', we can just multiply both eigenvectors
7
# by -1, and they will remain orthonormal eigenvectors of 'AAt' associated with the same eigenvalues
8

9
Ua = matrix(c(-1 * eu$vectors), nc=2, nr=2)
10
Ua

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# calculating V
2
AtA = t(A) %*% A
3
ev = eigen(AtA)
4
Va = matrix(c(ev$vectors), nc=2, nr=2)
5
Va

1
##           [,1]       [,2]
2
## [1,] 0.7071068 -0.7071068
3
## [2,] 0.7071068  0.7071068

1
# calculate E
2
# looking at the eigenvalues
3
eu$values

1
## [1] 2 0

1
ev$values

1
## [1] 2 0

1
# the only non-zero eigenvalue is 2,
2
# so the only singular value put in the E matrix is sqrt(2)
3
# so E is just sqrt(2) on the main diagonal, and the rest is zeros
4
# E needs to be size 'm' by 'n', so 2 by 2
5
Ea = matrix(c(c(sqrt(eu$values)), c(0, 0)), nc=2, nr=2)
6
Ea

1
##          [,1] [,2]
2
## [1,] 1.414214    0
3
## [2,] 0.000000    0

1
# now the result of U %*% E %*% t(V) should be equal to A again
2
A == Ua %*% Ea %*% t(Va)

1
##      [,1] [,2]
2
## [1,] TRUE TRUE
3
## [2,] TRUE TRUE

1
# using the SVD function built in to calculate svd of A
2
svdA = svd(A)
3
# U is
4
svdA$u

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# V is
2
svdA$v

1
##           [,1]       [,2]
2
## [1,] 0.7071068 -0.7071068
3
## [2,] 0.7071068  0.7071068

1
# the singular values of A are
2
svdA$d

1
## [1] 1.414214 0.000000

1
# set up the diagonal E matrix using the non-zero singular values svdA$d
2
# E is 'm' by 'n', so 2 by 2, with the singular values on the diagonal
3
E = diag(svdA$d, nc=2, nr=2)
4
E

1
##          [,1] [,2]
2
## [1,] 1.414214    0
3
## [2,] 0.000000    0

1
res = svdA$u %*% E %*% t(svdA$v)
2
res

1
##      [,1] [,2]
2
## [1,]    1    1
3
## [2,]    0    0

1
A

1
##      [,1] [,2]
2
## [1,]    1    1
3
## [2,]    0    0

1
# on my machine, this returns FALSE when comparing the 1's in 'res' to the 1's in 'A'
2
# I assume this is some floating point arithmetic error, but the value are both 1's of 'num' type
3
# I legitimately don't know why it does this, but taking the ceiling of both matrices works,
4
# so it must be some floating point error?
5
res == A

1
##       [,1]  [,2]
2
## [1,] FALSE FALSE
3
## [2,]  TRUE  TRUE

1
# what????
2
ceiling(res) == ceiling(A)

1
##      [,1] [,2]
2
## [1,] TRUE TRUE
3
## [2,] TRUE TRUE

SVD of B#

1
# calculate U
2
BBt = B %*% t(B)
3
eu = eigen(BBt)
4
# once again, the eigenvectors are negative for some reason, in this case we want a negative to get the '-2' in B
5
Ub = matrix(c(eu$vectors), nc=2, nr=2)
6
Ub

1
##      [,1] [,2]
2
## [1,]   -1    0
3
## [2,]    0   -1

1
# calculate V
2
BtB = t(B) %*% B
3
ev = eigen(BtB)
4
# these eigenvectors are all negative too,
5
# but since they're already negative in U, this will cause the 'SVD = -B'
6
# so multiply these eigenvectors by -1 (they're still othonormal eigenvectors of BtB)
7
Vb = matrix(c(-1 * ev$vectors), nc=2, nr=2)
8
Vb

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# calculate E
2
eu$values

1
## [1] 4 0

1
ev$values

1
## [1] 4 0

1
Eb = diag(c(c(sqrt(eu$values))), nc=2, nr=2)
2
Eb

1
##      [,1] [,2]
2
## [1,]    2    0
3
## [2,]    0    0

1
res = Ub %*% Eb %*% t(Vb)
2
res

1
##      [,1] [,2]
2
## [1,]   -2    0
3
## [2,]    0    0

1
res == B

1
##      [,1] [,2]
2
## [1,] TRUE TRUE
3
## [2,] TRUE TRUE

1
# using built in SVD to calculate svd of B
2
svdB = svd(B)
3
# U is
4
svdB$u

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# V is
2
svdB$v

1
##      [,1] [,2]
2
## [1,]   -1    0
3
## [2,]    0    1

1
# singular values of B are
2
svdB$d

1
## [1] 2 0

1
# now create the diagonal E matrix using the singular values on the main diagonal
2
# E has the same dimensions as B, so 2 by 2
3
E = diag(svdB$d, nc=2, nr=2)
4
E

1
##      [,1] [,2]
2
## [1,]    2    0
3
## [2,]    0    0

1
resB = svdB$u %*% E %*% t(svdB$v)
2
resB

1
##      [,1] [,2]
2
## [1,]   -2    0
3
## [2,]    0    0

1
B

1
##      [,1] [,2]
2
## [1,]   -2    0
3
## [2,]    0    0

1
resB == B

1
##      [,1] [,2]
2
## [1,] TRUE TRUE
3
## [2,] TRUE TRUE

SVD of C#

1
# calculate U
2
CCt = C %*% t(C)
3
eu = eigen(CCt)
4
# negative again, for some reason, you know the drill by now
5
Uc = matrix(c(-1 * eu$vectors), nc=2, nr=2)
6
Uc

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# calculate V
2
CtC = t(C) %*% C
3
ev = eigen(CtC)
4
# the second eigenvector is negative for some reason
5
ev$vectors[,2]

1
## [1]  0 -1  0

1
# make it non-negative, to make the SVD work (remains orthonormal eigenvector)
2
ev$vectors[,2] = -1 * ev$vectors[,2]
3
Vc = matrix(c(ev$vectors), nc=3, nr=3)
4
Vc

1
##           [,1] [,2]       [,3]
2
## [1,] 0.8944272    0  0.4472136
3
## [2,] 0.0000000    1  0.0000000
4
## [3,] 0.4472136    0 -0.8944272

1
# calculate E
2
eu$values

1
## [1] 5 4

1
# floating point error, round eigenvalues
2
round(ev$values, 14)

1
## [1] 5 4 0

1
# non-zero eigenvalues are 5 and 4
2
# size of E is the same as C, so a 2 x 3 matrix
3
Ec = diag(sqrt(eu$values), nc=3, nr=2)
4
Ec

1
##          [,1] [,2] [,3]
2
## [1,] 2.236068    0    0
3
## [2,] 0.000000    2    0

1
resC = Uc %*% Ec %*% t(Vc)
2
resC

1
##      [,1] [,2] [,3]
2
## [1,]    2    0    1
3
## [2,]    0    2    0

1
C

1
##      [,1] [,2] [,3]
2
## [1,]    2    0    1
3
## [2,]    0    2    0

1
# weird floating point stuff??
2
round(resC) == round(C)

1
##      [,1] [,2] [,3]
2
## [1,] TRUE TRUE TRUE
3
## [2,] TRUE TRUE TRUE

1
# using built in SVD to calculate svd of C
2
svdC = svd(C)
3
# U is
4
svdC$u

1
##      [,1] [,2]
2
## [1,]    1    0
3
## [2,]    0    1

1
# V is
2
svdC$v

1
##           [,1] [,2]
2
## [1,] 0.8944272    0
3
## [2,] 0.0000000    1
4
## [3,] 0.4472136    0

1
# this only gives two eigenvectors for V
2
# this is because the eigenvalue associated with the third eigenvector is 0
3
# so the third eigenvector is technically arbitrary, but I included it for matrix size consistency
4

5

6
# the singular values of C are
7
svdC$d

1
## [1] 2.236068 2.000000

1
# now create E, diagonal matrix w/ same dimensions as C, so a 2 x 3 matrix
2
# however, since 'svdC$v' only contains 2 vectors, E will have to be a 2 x 2 matrix
3
E = diag(svdC$d, nc=2, nr=2)
4
E

1
##          [,1] [,2]
2
## [1,] 2.236068    0
3
## [2,] 0.000000    2

1
resC = svdC$u %*% E %*% t(svdC$v)
2
resC

1
##      [,1] [,2] [,3]
2
## [1,]    2    0    1
3
## [2,]    0    2    0

1
C

1
##      [,1] [,2] [,3]
2
## [1,]    2    0    1
3
## [2,]    0    2    0

1
# float stuff again
2
round(resC) == round(C)

1
##      [,1] [,2] [,3]
2
## [1,] TRUE TRUE TRUE
3
## [2,] TRUE TRUE TRUE

Q2#

Compute the Pseudo-Inverses of the matrices given in Q1#

The pseudo-inverse of a matrix M uses the SVD of that matrix.

Given the SVD of M: M = U * E * V^T,

The pseudo-inverse of M, ‘M^t’ is:

M^t = V * E^t * U^T

Where E^t is the pseudo inverse of E, specifically, a diagonal matrix of the reciprocal of the non-zero singular values of ‘M’

Pseudo-Inverse of A#

1
# have U and V, need to make E^t
2

3
# get the eigenvalues
4
sing = eigen(A %*% t(A))$values
5
# invert the non-zero entries, take the square root
6
sing = sqrt(c(sing[sing != 0]^-1, sing[sing == 0]))
7
sing

1
## [1] 0.7071068 0.0000000

1
# make diagonal matrix of non-zero singular reciprocals
2
Et = diag(sing, nc=2, nr=2)
3
Et

1
##           [,1] [,2]
2
## [1,] 0.7071068    0
3
## [2,] 0.0000000    0

1
At = Va %*% Et %*% t(Ua)
2
At

1
##      [,1] [,2]
2
## [1,]  0.5    0
3
## [2,]  0.5    0

Pseudo-Inverse of B#

1
# have U and V, need to make E^t
2

3
# get the eigenvalues
4
sing = eigen(B %*% t(B))$values
5
# invert the non-zero entries, take the square root
6
sing = sqrt(c(sing[sing != 0]^-1, sing[sing == 0]))
7
sing

1
## [1] 0.5 0.0

1
# make diagonal matrix of non-zero singular reciprocals
2
Et = diag(sing, nc=2, nr=2)
3
Et

1
##      [,1] [,2]
2
## [1,]  0.5    0
3
## [2,]  0.0    0

1
Bt = Vb %*% Et %*% t(Ub)
2
Bt

1
##      [,1] [,2]
2
## [1,] -0.5    0
3
## [2,]  0.0    0

Pseudo-Inverse of C#

1
# have U and V, need to make E^t
2

3
# get the eigenvalues
4
sing = eigen(C %*% t(C))$values
5
# invert the non-zero entries, take the square root
6
sing = sqrt(c(sing[sing != 0]^-1, sing[sing == 0]))
7
sing

1
## [1] 0.4472136 0.5000000

1
# make diagonal matrix of non-zero singular reciprocals
2
# needs to be 3 x 2 to matrix multiply w/ V
3
Et = diag(sing, nc=2, nr=3)
4
Et

1
##           [,1] [,2]
2
## [1,] 0.4472136  0.0
3
## [2,] 0.0000000  0.5
4
## [3,] 0.0000000  0.0

1
Ct = Vc %*% Et %*% t(Uc)
2
Ct

1
##      [,1] [,2]
2
## [1,]  0.4  0.0
3
## [2,]  0.0  0.5
4
## [3,]  0.2  0.0

Singular Value Decomposition and Pseudo-Inverse#

Q1#