(A = matrix( c( 1, 2, 3, -1, 0, 4), nrow=2, byrow=TRUE))
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
Define X= \[AA^T\]
(X = A %*% t(A))
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17
Define Y= \[A^TA\]
(Y = t(A) %*% A)
## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25
Compute the eigenvalues and eigenvectors of X and Y
(eX = eigen(X))
## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
(eY = eigen(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
decompA = svd(A)
#u,v,d
print(c(list(decompA$u,decompA$v,decompA$d)))
## [[1]]
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043
##
## [[2]]
## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824
##
## [[3]]
## [1] 5.157693 2.097188
(u1 = decompA$u[,1])
## [1] -0.6576043 -0.7533635
(u2 = decompA$u[,2])
## [1] -0.7533635 0.6576043
(v1 = decompA$v[,1])
## [1] 0.01856629 -0.25499937 -0.96676296
(v2 = decompA$v[,2])
## [1] -0.6727903 -0.7184510 0.1765824
(d1 = decompA$d[1])
## [1] 5.157693
(d2 = decompA$d[2])
## [1] 2.097188
(X %*% u1)
## [,1]
## [1,] -17.49346
## [2,] -20.04083
(d1 * d1 * u1)
## [1] -17.49346 -20.04083
(X %*% u2)
## [,1]
## [1,] -3.313442
## [2,] 2.892274
(d2 * d2 * u2)
## [1] -3.313442 2.892274
(A %*% v1)
## [,1]
## [1,] -3.391721
## [2,] -3.885618
(d1 * u1)
## [1] -3.391721 -3.885618
(A %*% v2)
## [,1]
## [1,] -1.579945
## [2,] 1.379120
(d2 * u2)
## [1] -1.579945 1.379120
(Y %*% v1)
## [,1]
## [1,] 0.4938968
## [2,] -6.7834426
## [3,] -25.7176364
(d1 * d1 * v1)
## [1] 0.4938968 -6.7834426 -25.7176364
(Y %*% v2)
## [,1]
## [1,] -2.9590650
## [2,] -3.1598902
## [3,] 0.7766445
(d2 * d2 * v2)
## [1] -2.9590650 -3.1598902 0.7766445
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. The function myinverse should be correct and must use co-factors and determinant of A to compute the inverse.
myInverse = function(A){
#Check if matrix is square
if(dim(A)[1] != dim(A)[2]){ return('Matrix is not square!') }
#Check if matrix determinant = 0
if(det(A) == 0){ return('Matrix is singular!') }
#Create co-factor matrix
cofactorMatrix = A*0
#Co-factoring process
for (i in 1:ncol(A)) {
for (j in 1:nrow(A)) {
cofactorMatrix[i,j] = det(A[-i,-j]) * (-1)^(i+j)
}}
inverseMatrix = t((cofactorMatrix)/det(A))
return(inverseMatrix)
}
Test the function for a well-conditioned matrix
(A = matrix(c(1,2,4,4,3,1,8,7,1),nrow=3, byrow=TRUE))
## [,1] [,2] [,3]
## [1,] 1 2 4
## [2,] 4 3 1
## [3,] 8 7 1
(B = myInverse(A))
## [,1] [,2] [,3]
## [1,] -0.2 1.30 -0.50
## [2,] 0.2 -1.55 0.75
## [3,] 0.2 0.45 -0.25
(C = solve(A))
## [,1] [,2] [,3]
## [1,] -0.2 1.30 -0.50
## [2,] 0.2 -1.55 0.75
## [3,] 0.2 0.45 -0.25
round(B,5)==round(C,5)
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE