PROBLEM SET 1
Verify that SVD and Eigenvalues of matrix A are related.
# Set output values to round off to 2 digits
options(digits=2)
# Matrix A
A <- matrix(c(1,-1,2,0,3,4), ncol=3)
A
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4
# Compute X = A x transpose of A
X <- A %*% t(A)
# Compute Y = transpose of A x A
Y <- t(A) %*% A
# Eigenvalues / eigenvectors of X
Xegnvls <- round(eigen(X)$values, digits=3)
Xegnvts <- eigen(X)$vectors
# Eigenvalues / eigenvectors of Y
Yegnvls <- round(eigen(Y)$values, digits=3)
Yegnvts <- eigen(Y)$vectors
# Rank of matrix A
rank <- qr(A)$rank
# SVD of A
Asvd <- svd(A)
Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y.
Left-singular = U = columns of X eigenvectors
# Columns of X eigenvectors
Xegnvts
## [,1] [,2]
## [1,] 0.66 -0.75
## [2,] 0.75 0.66
# Columns of left-singular matrix
Asvd$u
## [,1] [,2]
## [1,] -0.66 -0.75
## [2,] -0.75 0.66
The first column vector of U is multiple of -1 of first eigenvector of X. The second columns are are exactly the same in X and U.
Left-singular = V = columns of X eigenvectors
# Columns of Y eigenvectors
Yegnvts
## [,1] [,2] [,3]
## [1,] -0.019 0.67 0.74
## [2,] 0.255 0.72 -0.65
## [3,] 0.967 -0.18 0.18
# Columns of left-singular matrix
Asvd$v
## [,1] [,2]
## [1,] 0.019 -0.67
## [2,] -0.255 -0.72
## [3,] -0.967 0.18
The first column vector of V is multiple of -1 of first eigenvector of X. The second columns are exactly the same in Y and V.
Show that two non-zero eigenvalues (the 3rd value will be very close to zero, if not zero) of both X and Y are the same and squares of non-zero singular values of A
# Eigenvalues of X
Xegnvls
## [1] 26.6 4.4
# Eigenvalues of Y
Yegnvls
## [1] 26.6 4.4 0.0
# Square roots of singular values of A
(Asvd$d)^2
## [1] 26.6 4.4
The Y eigenvalues contain the same eigenvalues as X and a third eigenvalue that is zero. The squre of singular values are equal to eigenvalues of X and Y.
PROBLEM SET 2
Create function to compute inverse of well-conditioned full-rank square marix using co-factors.
### Function: myinverse -------->
myinverse <- function(m) {
# check if matrix is square
if(nrow(m) != ncol(m)) { return('provide square matrix only')}
# intialize co-factor matrix with zeros
cofactor <- em <- mat.or.vec(nrow(m), ncol(m))
# Find co-factor matrix
for ( i in 1:nrow(m)){
for(j in 1:nrow(m)) {
c <- (-1)^(i+j) * det(m[-i,-j])
cofactor[i,j] <- c
}
}
cat('Co-factor matrix:\n')
print(cofactor)
cat('\n')
# determinant of A:
# can be calculated by dot product of row i of A with column of i of transpose(cofactor) matrix
detm <- (m[1,] %*% t(cofactor)[,1])[1,1]
cat('det(A) = ', detm, '\n')
inv <- t(cofactor) / detm
return(inv)
}
### <-------- end function: myinverse
Test with 3X3 matrix
A <- matrix(c(1,0,5,2,1,6,3,4,0), ncol=3)
A
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 0 1 4
## [3,] 5 6 0
B <- myinverse(A) # my inverse function
## Co-factor matrix:
## [,1] [,2] [,3]
## [1,] -24 20 -5
## [2,] 18 -15 4
## [3,] 5 -4 1
##
## det(A) = 1
B
## [,1] [,2] [,3]
## [1,] -24 18 5
## [2,] 20 -15 -4
## [3,] -5 4 1
# check for A * B = I
round(A %*% B, digits=2)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
Test with 4X4 matrix
A <- matrix(c(1,3,-1,-3,-2,-9,2,-6,-2,0,4,26,-3,-9,7,2), ncol=4)
B <- myinverse(A) # my inverse function
## Co-factor matrix:
## [,1] [,2] [,3] [,4]
## [1,] 882 396 201 -102
## [2,] -212 -94 -48 24
## [3,] 90 42 21 -12
## [4,] 54 24 12 -6
##
## det(A) = -6
B
## [,1] [,2] [,3] [,4]
## [1,] -147 35 -15.0 -9
## [2,] -66 16 -7.0 -4
## [3,] -34 8 -3.5 -2
## [4,] 17 -4 2.0 1
# check for A * B = I
round(A %*% B, digits=2)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1