A <- matrix(c(1, 2, 3, -1, 0, 4), nrow = 2, byrow = TRUE)
print(A)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
X <- A %*% t(A)
print(X)
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17
Y <- t(A) %*% A
print(Y)
## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25
Computing the Eigen values and Eigen vectors of X and Y
lambda_x <- eigen(X)$values
print(lambda_x)
## [1] 26.601802 4.398198
lambda_y <- eigen(Y)$values
print(lambda_y)
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
s_x <- eigen(X)$vectors
print(s_x)
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
s_y <- eigen(Y)$vectors
print(s_y)
## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001
Generating SVD:
svd_a <- svd(A)
print(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
As we see above, $u and $v are same as the Eigenvectors given above as these can be easily generated by a scalar multiplication of the Eigenvectors.
library(Matrix)
myinverse = function(A)
{
### First check - is the input matrix a full-rank square matrix
if (!(nrow(A) == ncol(A) & rankMatrix(A) == nrow(A)))
{
return("Not a full-rank square matrix")
}
else
{
Asize = nrow(A)
## Create a dummy matrix to create the matrix of the co-factors
Cofactor_matrix <- matrix(nrow = Asize, ncol = Asize)
for (i in 1:Asize) {
for (j in 1:Asize) {
Cofactor_matrix[i,j] <- (-1)^(i+j) * det(A[-i,-j])
}
}
## Now we have the cofactor matrix ready, we will use the fornula to get inverse of A:
## A inverse = (Transpose of Cofactor matrix of A) / det(A)
inverse_matrix <- t(Cofactor_matrix) / det(A)
return(inverse_matrix)
}
}
Testing:
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
test_matrix <- matrix(c(2,3,2,1,2,1,1,1,2),nrow = 3)
test_matrix
## [,1] [,2] [,3]
## [1,] 2 1 1
## [2,] 3 2 1
## [3,] 2 1 2
solve(test_matrix)
## [,1] [,2] [,3]
## [1,] 3 -1 -1
## [2,] -4 2 1
## [3,] -1 0 1
myinverse(test_matrix)
## [,1] [,2] [,3]
## [1,] 3 -1 -1
## [2,] -4 2 1
## [3,] -1 0 1
identical(solve(test_matrix) %>% round(6), myinverse(test_matrix) %>% round(6))
## [1] TRUE