##Packages
library(pracma) #for rref
library(matlib) #for cofactor, going to use this as a check in part 2
#inputting the given matrix
A=matrix(c(1,-1,2,0,3,4),nrow = 2)
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
Now computing the eigenvalues and eigenvectors of X and Y. The first two eigen values of Y are equal to those of X.
eX<-eigen(X)
eX
## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
eY<-eigen(Y)
eY
## 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
Now showing the SVD of A
sA<-svd(A)
sA
## $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
dsA<-sA$`d` #singular values
usA<-sA$u #left singular vectors
vsA<-sA$v #right singular vectors
dsA^2
## [1] 26.601802 4.398198
usA and eX$vectors should be equal, but there is a factor of -1 on the first vector. Similary there is a -1 factor on the first vsA and eY$vectors. I have not ben able to determine why this is, although given that the vectors make the same space, this should not be an issue. The square of the singular values and the eigen values are equal.
Firts thing is creating a sample matrix to try this out with when it works.
c1<-c(1:4,7:10,1,3,5,11,-1,-2,0,1)
A<-matrix(c1,nrow=4)#our sample matrix
rref(A) #showing that it has rank 4
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
I wrote two functions to create a cofactor matrix. The big one is just two nested loops. It’s not pretty, but it works. The other one gives the determinant of a submatrix of A.
createCmatrix<- function(A){
CA<-matrix(c(1:nrow(A)^2),nrow=nrow(A))
for(i in 1:nrow(A)){
for(j in 1:ncol(A)){
CA[i,j]<-minidet(A,i,j)*(-1)^(i+j)
}
}
return(CA)
}
minidet<-function(A,i,j){
B<-A[-i,-j]
db<-det(B)
return(db)
}
Now, creating the cofactor matrix, comparing it to a couple of test cases and then proving that we have an inverse.
CF<-createCmatrix(A)
CF
## [,1] [,2] [,3] [,4]
## [1,] -85 25 6 2.400000e+01
## [2,] 23 -11 6 -4.800000e+01
## [3,] 89 -5 -30 2.400000e+01
## [4,] -39 3 18 2.664535e-15
cofactor(A,1,2)
## [1] 25
cofactor(A,3,2) #comparing my cofactor matrix to the one created by the function
## [1] -5
invA<-t(CF)/det(A) #creating the invserse matrix
resul<-invA%*%A #multiplying the inverse by A. This results in something that is just about the identity matrix, within the bounds of
#floating point math
resul
## [,1] [,2] [,3] [,4]
## [1,] 1.000000e+00 -7.105427e-15 0.000000e+00 8.881784e-16
## [2,] -5.551115e-17 1.000000e+00 -5.551115e-17 -1.387779e-17
## [3,] 2.220446e-16 8.881784e-16 1.000000e+00 -5.551115e-17
## [4,] 2.590520e-16 8.141636e-16 6.291264e-16 1.000000e+00
ident<-diag(4) #creating a 4x4 identity matrix
ident
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
all.equal(ident,resul) #finally, verifying that resul and ident are just about the same, despite tiny differences
## [1] TRUE