IS 605 FUNDAMENTALS OF COMPUTATIONAL MATHEMATICS - 2019
write code in R to 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. 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.
library(pracma)
## Warning: package 'pracma' was built under R version 3.5.2
a <- matrix(c(1,2,3,-1,0,4), byrow =T, ncol =3 )
a
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
Compute X=AAT and Y=ATA
x <- a %*% t(a)
y <- t(a) %*% a
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
Compute eigenvalues and eigenvectors of X and Y
eigenvalues_x <- eigen(x)
eigenvalues_y <- eigen(y)
eigenvalues_x
## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
eigenvalues_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
Compute the left-singular(u), singular values(d), and right-singular(v) using the SVD command.
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
Examine the two sets of singular vectors and show they are eigenvectors of x and y.
eigenvalues_x$values
## [1] 26.601802 4.398198
eigenvalues_x$vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
eigenvalues_y$values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
eigenvalues_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
myinv <- function(a){
dim_a <- dim(a)[1]
det_a <- det(a)
c <- a
for(i in 1: dim_a){
for(x in 1:dim_a){
d <- a[-i,-x]
det_d <- det(d)
c[i,x] <- (-1)^(i+x) *det_d
}
}
b <- t(c)
b <- b/det_a
return(b)
}
aa <- matrix(c(1,4,0,8,15,3,1,9,2), byrow =T, ncol =3 )
myinv(aa)
## [,1] [,2] [,3]
## [1,] -0.06122449 0.16326531 -0.24489796
## [2,] 0.26530612 -0.04081633 0.06122449
## [3,] -1.16326531 0.10204082 0.34693878