In this problem, we’ll verify using R that SVD and Eigenvalues are related as worked out in the weekly module. Given a 3 × 2 matrix A
A <- matrix(c(1,-1,2,0,3,4),nrow=2, ncol = 3)
(A)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
compute X = AAT and Y = ATA. Then, compute the eigenvalues and eigenvectors of X and Y using the built-in commans in R.
Then, compute the left-singular, singular values, and right-singular vectors of A using the svd command.
Create a function to transpose matrix
myTranspose <- function (A) {
#for matrix A with m x n
#create a new empty matrix with n x m for transpose
A_T <- matrix(, nrow = ncol(A), ncol = nrow(A))
# transpose of A [i,j] <==> A[j,i]
for(i in 1:nrow(A)) {
for(j in 1:ncol(A)) {
A_T[j,i] <- A[i,j]
}
}
return(A_T)
}
(myTranspose(A))
## [,1] [,2]
## [1,] 1 -1
## [2,] 2 0
## [3,] 3 4
Try X = A times transpose of A and Y = transpose of A times A
computeXY <- function (A) {
X <- A %*% myTranspose(A)
Y <- myTranspose(A) %*% A
return(list(X=X,Y=Y))
}
XY <- computeXY(A)
(XY)
## $X
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17
##
## $Y
## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25
Calculate eigenvalues and vectors of X and Y
eigen_X <- eigen(XY$X) #built-in in R for eigenvalues and eigenvectors calculation
(eigen_X)
## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
eigen_Y <- eigen(XY$Y) #built-in in R for eigenvalues and eigenvectors calculation
(eigen_Y)
## eigen() decomposition
## $values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001
Now, calculate SVD
svd_A <- (svd(A))
(svd_A)
## $d
## [1] 5.157693 2.097188
##
## $u
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043
##
## $v
## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824
square of singular value of A is equal to eigen value of X
(svd_A$d*svd_A$d)
## [1] 26.601802 4.398198
(eigen_X$values)
## [1] 26.601802 4.398198
the singular vectors are simlar but the sign of the value and number of columns are different
(svd_A$u)
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043
(eigen_X$vectors)
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
(svd_A$v)
## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824
(eigen_Y$vectors)
## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001
Write a invese function with Using co-factors method and prove AxB = I where matrix B is invese of matrix A and I is identify matrix
co-factor method
https://www.youtube.com/watch?v=g7TFJUJXErU
step 1: minor matrix step 2: get cofactor matrix step 3: get adjoint matrix of cofactor matrix setp 4: inverse matrix = adjoint matrix divided by determine of matrix
Define a 3x3 matrix B
B <- matrix(c(1,4,5,3,0,3,-2,1,2),3,3)
(B)
## [,1] [,2] [,3]
## [1,] 1 3 -2
## [2,] 4 0 1
## [3,] 5 3 2
Create inverse function
myInverse <- function(M) {
#combine step 1 and step2 to get cofactor matrix. i don't know how to store matrix of minors in matrix, so i combine both step 1 and step 2.
#remove row and colunm in R
#https://stackoverflow.com/questions/12919984/how-to-delete-specific-rows-and-columns-from-a-matrix-in-a-smarter-way
cofactor_matrix <- matrix(,nrow=nrow(M),ncol = ncol(M))
for(i in 1:nrow(M)) {
for(j in 1:ncol(M)) {
tmpM <- M
tmpM <- tmpM[-i,]
tmpM <- tmpM[,-j]
cofactor_matrix[i,j] <- (-1)^(i+j) * det(tmpM)
}
}
#step3: adjoint matrix of M, where adj(M) = transpose of cofactor matrix
adjoint_M <- myTranspose(cofactor_matrix)
#step4: invese of M = adj(M) / det(M)
return(adjoint_M / det(M))
}
Calculate inverse of B
inverse_B <- myInverse(B)
(inverse_B)
## [,1] [,2] [,3]
## [1,] 0.08333333 0.3333333 -0.08333333
## [2,] 0.08333333 -0.3333333 0.25000000
## [3,] -0.33333333 -0.3333333 0.33333333
Try B times inverse of B = I, round digit to 2 avoid precision error
(round(B %*% inverse_B,2))
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1