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 \(\begin{bmatrix}1&2&3\\-1&0&4\end{bmatrix}\) write code in R to compute X = A\({A}^T\) and Y = \({A}^T\)A. 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
Here we are calculating first transpose of matrix and then values of X and Y
A <- matrix(c(1,-1,2,0,3,4), nrow=2) # given matrix
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4A_T <- t(A) #calculating transpose using t() function
A_T## [,1] [,2]
## [1,] 1 -1
## [2,] 2 0
## [3,] 3 4X <- A %*% A_T #calculating value of X
X## [,1] [,2]
## [1,] 14 11
## [2,] 11 17Y <- A_T %*% A # calulating Y
Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25Now calculating eigenvectors and eigenvalues using R functions eigen()
library(pracma)
eigenX <- eigen(X) # eigen() - decomposition of matrix of matrix X
eigen_valueX <- eigenX$values #this will give eigen value of X
paste0("eigenvalue of X = ",eigen_valueX)## [1] "eigenvalue of X = 26.6018016555873"
## [2] "eigenvalue of X = 4.39819834441274"eigen_vectorX <- eigenX$vectors # this will give eigen vectors for X
eigen_vectorX## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043eigenY <- eigen(Y) # eigen() - decomposition of matrix of matrix Y
eigen_valueY <- eigenY$values #this will give eigen value of Y
paste0("eigenvalue of Y = ",eigen_valueY)## [1] "eigenvalue of Y = 26.6018016555872"
## [2] "eigenvalue of Y = 4.39819834441274"
## [3] "eigenvalue of Y = 1.05898195693058e-16"eigen_vectorY <- eigenY$vectors # this will give eigen vectors for Y
eigen_vectorY## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001Calculating left Singular values
singular_val <- svd(A) # Compute the singular-value decomposition of a rectangular matrix
left_sing_val <- svd(A,nu = 2,nv = 0) # setting nu= 2 gives left singular values
left_sing_val## $d
## [1] 5.157693 2.097188
##
## $u
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043Calculating right Singular values
right_sing_val <- svd(A,nu = 0,nv = 3) # # setting nv= 3 gives right singular values
right_sing_val## $d
## [1] 5.157693 2.097188
##
## $v
## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001Now comparing Results of left singular vectors with eigen vectors X
left_sing_val$u## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043eigen_vectorX## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043Values looks similar other than the negative sign
Now comparing Results of right singular vector with eigen vectors Y
right_sing_val$v## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001eigen_vectorY## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001Again both values looks similar other than 1 negative
Now comparing with eigen values X and left singular values
eigen_valueX## [1] 26.601802 4.398198(left_sing_val$d)^2## [1] 26.601802 4.398198Comparing with eigen values Y and right singular values
eigen_valueY## [1] 2.660180e+01 4.398198e+00 1.058982e-16(right_sing_val$d)^2## [1] 26.601802 4.398198the non-zero eigenvalues of X and Y are the same as the non-zero singular values.
Using the procedure outlined in section 1 of the weekly handout, write a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors. In order to compute the co-factors, you may use built-in commands to compute the determinant. Your function should have the following signature: B = myinverse(A) where A is a matrix and B is its inverse and AxB=I. The off-diagonal elements of I should be close to zero, if not zero. Likewise, the diagonal elements should be close to 1, if not 1. Small numerical precision errors are acceptable but the function myinverse should be correct and must use co-factors and determinant of A to compute the inverse
myinverse <- function(A) {
n <- dim(A)[1] #get dimesnsion of square matrix
C <- matrix(0, n, n)
for(i in 1:n) { # loop for rows
for(j in 1:n) { # loop for colummns
C[i,j] <- (-1)^(i+j)*det(matrix(A[-i,-j], n-1))
}
}
# Inverse of A is equal to transpose of C divided by determinant of A
B <- t(C)/det(A)
return(B)
}
A <- matrix(c(2,1,1,3,2,1,2,1,2), nrow=3 ,byrow=TRUE)
# Calculate inverse using 'myinverse' function
B <- myinverse(A)
B## [,1] [,2] [,3]
## [1,] 3 -1 -1
## [2,] -4 2 1
## [3,] -1 0 1library(MASS) ## Warning: package 'MASS' was built under R version 3.4.3ginv(A) # checking inverse values using ginv() function built under MASS package## [,1] [,2] [,3]
## [1,] 3 -1.00000e+00 -1
## [2,] -4 2.00000e+00 1
## [3,] -1 7.21645e-16 1Checking if A * inv(A) gives Identity matrix
I <- round(A %*% B, digits=10) # rounding didgits 10 gives readable solution else sometimes the value is too long
I## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1A <- matrix(c(1,2,3,4), nrow=2)
# Calculate inverse using 'myinverse' function
B <- myinverse(A)
B## [,1] [,2]
## [1,] -2 1.5
## [2,] 1 -0.5ginv(A)## [,1] [,2]
## [1,] -2 1.5
## [2,] 1 -0.5I <- round(A %*% B, digits=10) # rounding didgits 10 gives readable solution else sometimes the value is too long
I## [,1] [,2]
## [1,] 1 0
## [2,] 0 1