Problem set 1
A <- matrix(data = 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
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
Compute Eigenvalues and Eigenvectors of X and Y
eigen_X <- eigen(X)
eigen_values_X <- round(eigen_X$values, 3)
eigen_values_X
## [1] 26.602 4.398
eigen_Y <- eigen(Y)
eigen_values_Y <- round(eigen_Y$values, 3)
eigen_values_Y
## [1] 26.602 4.398 0.000
eigen_vectors_X <- eigen_X$vectors
eigen_vectors_X
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043
eigen_vectors_Y <- eigen_Y$vectors
eigen_vectors_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
Compute left-singular, singular values, and right-singular vectors of A
svd_A <- svd(A)
left_singular_A <- svd_A$u
left_singular_A
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043
singular_values_A <- round(svd_A$d, 3)
singular_values_A
## [1] 5.158 2.097
right_singular_A <- svd_A$v
right_singular_A
## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824
Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y. In addition, the 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 are squares of the non-zero singular values of A.
round(svd(X)$d, 3)
## [1] 26.602 4.398
round(eigen_values_X, 3)
## [1] 26.602 4.398
round(svd(Y)$d, 3)
## [1] 26.602 4.398 0.000
round(eigen_values_Y, 3)
## [1] 26.602 4.398 0.000
squares <- round(singular_values_A ^ 2, 3)
squares
## [1] 26.605 4.397
Problem set 2
Using the procedure outlined in section 1 of the weekly handout, 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)
myinverse <- function(A){
dimension_A <- nrow(A)
matrix <- matrix(NA, nrow = dimension_A, ncol = dimension_A)
for (i in 1:dimension_A) {
for (j in 1:dimension_A) {
matrix[i,j] <- det(A[-i,-j]) * (-1) ^ (i + j)
}
}
tranpose_matrix <- t(matrix)
inverse_matrix <- tranpose_matrix / det(A)
return(inverse_matrix)
}
A <- matrix(c(1,2,3,6,5,6,7,8,9), nrow = 3)
B <- round(myinverse(A), 5)
B
## [,1] [,2] [,3]
## [1,] -0.25 -1 1.08333
## [2,] 0.50 -1 0.50000
## [3,] -0.25 1 -0.58333
round(A%*%B, 3)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1