library(pracma)In this problem, we’ll verify using R that SVD and Eigenvalues are related as worked out in the weekly module.
Given Matrix A:
A = matrix(c(1,2,3,-1,0,4), nrow=2, ncol=3, byrow=TRUE)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4Write code in R to compute \(X = AA^T\) and \(Y = A^TA\). Then, compute the eigenvalues and eigenvectors of X and Y using the built-in commands in R.
Then, compute the left-singular, singular values, and right-singular vectors of A using the svd command.
Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y.
In addition, the two non-zero eigenvalues (the 3rd value will be very close to zero, if not zero) of both X and Y are the same and are squares of the non-zero singular values of A.
Compute the left-singular, singular values, and right-singular vectors of A using the svd command.
d contains the singular valuesu contains the left singular vectors (rows).v contains the right singular vectors (columns).A_svd <- svd(A, nu=nrow(A), nv=ncol(A))
singular_values <- as.matrix(diag(A_svd$d, nrow=nrow(A), ncol=ncol(A)))
left_singluar_vectors <- as.matrix(A_svd$u)
right_singluar_vectors <- as.matrix(A_svd$v)The singular values from SVD function:
singular_values## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0The left singular vectors from SVD function:
left_singluar_vectors## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043The right singular vectors from SVD function:
right_singluar_vectors## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001Compare left singular vector by computing \(AA^T\) and taking its eigenvector:
X <- A %*% t(A)
X_eigen <- eigen(X) #AAT
X_left_singluar_vectors <- as.matrix(X_eigen$vectors)
X_eigenvalues <- X_eigen$values The SVD left singular vectors and the eigenvectors of \(AA^T\) match but there are some sign differences. According to the sources linked below, these sign differences are to be expected.
http://rpubs.com/aaronsc32/singular-value-decomposition-r
http://rpubs.com/aaronsc32/singular-value-decomposition-r
left_singluar_vectors## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043X_left_singluar_vectors## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043The square root of the eigenvalues of X match the singular values of the SVD.
singular_values## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0diag(sqrt(X_eigenvalues), nrow=nrow(A), ncol=ncol(A))## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0Compare right singular vector by computing \(A^TA\) and taking its eigenvector:
Y <- t(A)%*% A #ATA
Y_eigen <- eigen(Y)
Y_right_singluar_vectors <- as.matrix(Y_eigen$vectors)
Y_eigenvalues <- Y_eigen$valuesThe SVD right singular vectors and \(A^TA\) eigenvectors match with sign differences. As mentioned earlier, sign differences are to be expected.
right_singluar_vectors## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001Y_right_singluar_vectors## [,1] [,2] [,3]
## [1,] -0.01856629 0.6727903 0.7396003
## [2,] 0.25499937 0.7184510 -0.6471502
## [3,] 0.96676296 -0.1765824 0.1849001The eigenvalues of X [\(AA^T\)] and Y [\(A^TA\)] match for nonzero values. Zero values are to be ignored.
round(X_eigenvalues,4)## [1] 26.6018 4.3982round(Y_eigenvalues,4)## [1] 26.6018 4.3982 0.0000The code below assumes that the matrix is square and that it is invertible.
#The function to get cofactors
get_cofactors <- function(M) {
cofactors <- M
for(i in 1:nrow(M)){
for(j in 1:ncol(M)){
cofactors[i,j] <- (det(M[-i,-j])*(-1)^(i+j))
}
}
return(cofactors) # output of cofactors matrix
}
#The function to get the inverse of a square matrix that is invertible
myinverse <- function(M){
M_cofactors <- get_cofactors(M)
M_cofactors_transpose <- t(M_cofactors)
M_det <- det(M)
return(M_cofactors_transpose/M_det)
}
#Example of an invertible matrix
C <- matrix(c(5,6,6,8,2,2,2,8,6,6,2,8,2,3,6,7), nrow=4, ncol=4, byrow=TRUE)
C## [,1] [,2] [,3] [,4]
## [1,] 5 6 6 8
## [2,] 2 2 2 8
## [3,] 6 6 2 8
## [4,] 2 3 6 7#Find the inverse of the example
C_inverse <- myinverse(C)
#Confirm the inverse of matrix * matrix results in identity
C_inverse %*% C## [,1] [,2] [,3] [,4]
## [1,] 1.000000e+00 4.973799e-14 4.263256e-14 4.973799e-14
## [2,] -6.750156e-14 1.000000e+00 -1.136868e-13 -1.350031e-13
## [3,] -4.440892e-15 -8.881784e-16 1.000000e+00 -8.881784e-16
## [4,] -8.881784e-16 -1.776357e-15 -2.664535e-15 1.000000e+00