# Given a 3 × 2 matrix A
A<- matrix(c(1,2,3,-1,0,4), nrow=2, ncol=3, byrow=TRUE)
A
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
# write code in R to compute X = A%*%t(A) and Y = t(A)%*%A.
# Computing X as AA(transpose)
X<- A%*%t(A)
X
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17
# Computing Y as A(transpose)A
Y<- t(A)%*%A
Y
## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25
# Then, compute the eigenvalues and eigenvectors of X and Y using the built-in commands in R
U_EIGEN <- eigen(X)$vectors
U_EIGEN
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
# Compute the svd vectors for U
U_svd <- svd(A)$u
# Compute the eigenvalues and eigenvectors of Y
V_EIGEN <- cbind(eigen(Y)$vectors[,1:2],c(0,0,0))
V_EIGEN
## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0
## [2,] 0.25499937 -0.7184510 0
## [3,] 0.96676296 0.1765824 0
# Compute the svd vectors for V
V_svd <- cbind(svd(A)$v,c(0,0,0))
V_svd
## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 0
## [2,] -0.25499937 -0.7184510 0
## [3,] -0.96676296 0.1765824 0
# Calculating Sigma...
# W is square root diagonal eigenvalues of t(A)*A extended with zero valued rows
W <- sqrt(diag(eigen(Y)$values))
Sigma_EIGEN_V<- matrix(c(W[1],0,0, 0,W[2,2],0), nrow=2, ncol=3, byrow=TRUE)
Sigma_EIGEN_V
## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0
S<- sqrt(diag(eigen(X)$values))
Sigma_EIGEN_U<- matrix(c(S[1:2],0, S[2,],0), nrow= 2, byrow=TRUE )
Sigma_EIGEN_U
## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0
Sigma_svd <- matrix(c(svd(A)$d[1],0,0,svd(A)$d[2],0,0),nrow=2)
Sigma_svd
## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0
# Decomposition and validation
U_svd %*% Sigma_svd %*% t(V_svd)
## [,1] [,2] [,3]
## [1,] 1 2.000000e+00 3
## [2,] -1 1.110223e-16 4
U_EIGEN %*% Sigma_EIGEN_U %*% t(V_EIGEN)
## [,1] [,2] [,3]
## [1,] 1 2.000000e+00 3
## [2,] -1 -9.992007e-16 4
# or using sigma eigen value of V.
U_EIGEN %*% Sigma_EIGEN_V %*% t(V_EIGEN)
## [,1] [,2] [,3]
## [1,] 1 2.000000e+00 3
## [2,] -1 -2.220446e-15 4
#require(functional)
library(functional)
myinverse <- function(A) {
# check if the matrix is square
if(length(unique(dim(A))) != 1){
stop("Matrix is not square. Please provide a square matrix")
return("Matrix is not square. Please provide a square matrix")
}
else
# check if the matrix is not invertable
if(det(A) == 0) {stop ("the determinant is zero. the matrix is not Invertable ")
return("the determinant is zero. the matrix is not Invertable ")
}
else {
# Initialize the Coefficient
C <- diag(0,nrow(A), ncol(A))
# loop to build the cofactor matrix, C
for(i in 1:nrow(A)) {
for(j in 1:ncol(A)){
C[i, j] <- (-1)^(i+j) * det(A[-i, -j])
}
}
# Using the formula from week 4 module 1: A = t(C)/det(A)
return((t(C) / det(A)))
}
}
############
# Tesing for for zero deternimant
A<- matrix(c(1,0,-1,0,1,0,-1,0,1), nrow=3, ncol=3, byrow=TRUE)
A
## [,1] [,2] [,3]
## [1,] 1 0 -1
## [2,] 0 1 0
## [3,] -1 0 1
B <- round(myinverse(A),2)
## Error in myinverse(A): the determinant is zero. the matrix is not Invertable
B
## Error in eval(expr, envir, enclos): object 'B' not found
# Testing for non square matrix
A<- matrix(c(1,2,3,4,9,6,7,8),nrow=4, ncol=2, byrow=TRUE)
A
## [,1] [,2]
## [1,] 1 2
## [2,] 3 4
## [3,] 9 6
## [4,] 7 8
B <- round(myinverse(A),2)
## Error in myinverse(A): Matrix is not square. Please provide a square matrix
B
## Error in eval(expr, envir, enclos): object 'B' not found
# Testing valid matrices
set.seed(100)
A <- matrix(sample(0:20, 9, replace=TRUE),byrow=T, nrow =3)
A
## [,1] [,2] [,3]
## [1,] 6 5 11
## [2,] 1 9 10
## [3,] 17 7 11
B <- round(myinverse(A),2)
B
## [,1] [,2] [,3]
## [1,] -0.05 -0.03 0.08
## [2,] -0.25 0.19 0.08
## [3,] 0.23 -0.07 -0.08
# Validating if A%*%B = I.
round(A%*%B,1)
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0
## [2,] 0.0 1.0 0
## [3,] -0.1 0.1 1
##
set.seed(100)
A <- matrix(sample(0:20, 16, replace=TRUE),byrow=T, nrow =4)
A
## [,1] [,2] [,3] [,4]
## [1,] 6 5 11 1
## [2,] 9 10 17 7
## [3,] 11 3 13 18
## [4,] 5 8 16 14
B <- round(myinverse(A),2)
B
## [,1] [,2] [,3] [,4]
## [1,] -0.09 0.17 0.09 -0.19
## [2,] -0.40 0.39 -0.06 -0.09
## [3,] 0.33 -0.27 -0.02 0.14
## [4,] -0.12 0.03 0.03 0.03
# Validating if A%*%B = I.
round(A%*%B,1)
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0.0 0.1 0
## [2,] 0.0 1.1 0.1 0
## [3,] -0.1 0.1 1.1 0
## [4,] -0.1 0.1 0.1 1