A <- matrix(c(1,2,3,-1,0,4), byrow=T, 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
#AA^T
X <- A %*% t(A)
X
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17
eigen_x <- eigen(X)
eigen_x
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
#A^TA
Y <- t(A) %*% A
Y
## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25
eigen_y <- eigen(Y)
eigen_y
## $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, singular values, and right-singular vectors of A 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
#validate/display that these eigenvalues are equal
svd_a$d
## [1] 5.157693 2.097188
eigen_x$values
## [1] 26.601802 4.398198
eigen_y$values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
#validate/display that these eigenvectors are equal
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
#validate/display that these eigenvectors are equal
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
#verify that the non-zero singular values of A, squared, are equal to the eigen values of X and Y
5.157693^2
## [1] 26.6018
2.097188^2
## [1] 4.398198
In the output, d represents a vector containing the singular values of A. U represents a matrix whose columns contains the left singular vectors of A. V represents a matrix whose columns contain the right singular vectors of A.
We can also see here that the two sets of singular vectors are eigenvectors of X and Y.
#establish the square matrix A
A <- matrix(c(1,2,3,-1,0,4,5,2,3), byrow=T, nrow=3, ncol=3)
#find the minor and cofactor of matrix A
minor <- function(A, i, j) det( A[-i,-j] )
cofactor <- function(A, i, j) (-1)^(i+j) * minor(A,i,j)
#find the adjugate of matrix A
adjugate2 <- function(A) {
n <- nrow(A)
t(outer(1:n, 1:n, Vectorize(
function(i,j) cofactor(A,i,j)
)))
}
adj_A<- adjugate2(A)
#calculate the inverse
myinverse <- function(A){
(1/det(A))*(adj_A)
}
#compute the inverse of A using the myinverse function
B <- myinverse(A)
B
## [,1] [,2] [,3]
## [1,] -0.25000 0.000 0.25000
## [2,] 0.71875 -0.375 -0.21875
## [3,] -0.06250 0.250 0.06250
#verify that A x B = I
round(A %*% B, 2)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
As you can see, \[ A x B = I \]. We have therefore verified that we’ve designed a myinverse function that properly calculates the inverse of matrix A.