SVD EigenValues and MatrixInverse Computation
Load Packages
library(knitr)
library(matlib)Problem Set 1
Answer: I’ll first code the matrix A and its transpose (will call it A_T), in the following code chunk.
A <- matrix(c(1, -1, 2, 0, 3, 4), nrow = 2)
print(A)## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4A_T <- t(A)
print(A_T)## [,1] [,2]
## [1,] 1 -1
## [2,] 2 0
## [3,] 3 4Now, in order to compute X = A.A_T and Y = A_T.A, I’ll code a general purpose matrix multiplication function my_mat_mult(), in below code chunk.
my_mat_mult <- function(m1, m2)
{
m1_rows <- c(1:nrow(m1)); m1_cols <- c(1:ncol(m1))
m2_rows <- c(1:nrow(m2)); m2_cols <- c(1:ncol(m2))
#
if(ncol(m1) != nrow(m2)) {
stop('matrices are not conformable for multiplication')
}
#
product <- matrix(0, nrow(m1), ncol(m2))
#
for(i in m1_rows) {
for(j in m2_cols) {
for(k in m1_cols) {
product[i, j] = product[i, j] + m1[i, k] * m2[k, j]
}
}
}
return(product)
}Computing X and Y, by calling my_mat_mult().
X <- my_mat_mult(A, A_T)
print(X)## [,1] [,2]
## [1,] 14 11
## [2,] 11 17Y <- my_mat_mult(A_T, A)
print(Y)## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25Computing eigenvectors and eigenvalues of X and Y, by R’s built-in function eigen().
e_X <- eigen(X)
e_X## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043e_Y <- eigen(Y)
e_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.1849001Now, I’ll do Single valued Decomposition (SVD) in the following code chunk, using R’s built-in function, svd().
SVD_A <- svd(A)
SVD_A## $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.1765824Explanation: SVD_A contains the decomposed matrices of A, namely U, Sigma and V_transposed, in that order. Sigma contains a diagonal matrix in the top left, and zeros in the rest of the locations.
The column vectors of U are the left-singular vectors of A. In the above results, observe the two columns vectors under $u.
The diagonal elements of D (top left in Sigma) are singular values of A. Note that they not vectors, because they are just elements of the diagonal matrix, D. In the above results, observe the two numbers under $d.
The column vectors of V are the right-singular vectors of A. In the above results, observe the two columns vectors $v.
Now, I’ll compare the left-singular vectors of A (i.e. U or $u), with the eigen vectors of X. In comparison_table1, we observe that the two sets of vectors (without the signs) are equal each to each.
comparison_table1 <- cbind(SVD_A$u, e_X$vectors)
colnames(comparison_table1) <- c('SVD_A_U_Col1', 'SVD_A_U_Col2', 'e_X_Col1', 'e_X_Col2')
comparison_table1## SVD_A_U_Col1 SVD_A_U_Col2 e_X_Col1 e_X_Col2
## [1,] -0.6576043 -0.7533635 0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043 0.7533635 0.6576043In the following code chunk, I’ll compare the right-singular vectors of A (i.e. V or $v), with the eigen vectors of Y. In this case, we know that while Y has 3 eigen vectors, V has 2 vectors. So, if we ignore the 3rd vector (column heading ‘IGNORE’), then the two sets of vectors (without signs) are equal, each to each.
comparison_table2 <- cbind(SVD_A$v, e_Y$vectors)
colnames(comparison_table2) <- c('SVD_A_V_Col1', 'SVD_A_V_Col2', 'e_Y_Col1', 'e_Y_Col2', 'IGNORE')
comparison_table2## SVD_A_V_Col1 SVD_A_V_Col2 e_Y_Col1 e_Y_Col2 IGNORE
## [1,] 0.01856629 -0.6727903 -0.01856629 -0.6727903 0.7396003
## [2,] -0.25499937 -0.7184510 0.25499937 -0.7184510 -0.6471502
## [3,] -0.96676296 0.1765824 0.96676296 0.1765824 0.1849001Problem Set 2
Answer: In the following chode chunk, I’ll code myinverse() function, which accept a matrix as input, and then will first check if the input matrix is square and non-singular. If both, then it’ll return the inverse.
myinverse <- function(m)
{
m_rows <- c(1:nrow(m)); m_cols <- c(1:ncol(m))
#
if(nrow(m) != ncol(m)) {
stop('the matrix is not square')
}
#
if(det(m) == 0) {
stop('the matrix is singular')
}
#
cof_m <- matrix(0, nrow(m), ncol(m))
#
for(i in m_rows) {
for(j in m_cols) {
cof_m[i, j] <- cofactor(m, i, j)
}
}
#
trans_cof_m <- t(cof_m)
inv_m <- ( 1/det(m) ) * trans_cof_m
return(inv_m)
}Now, I’ll test myinverse() with the below matrix, which is a full-rank square matrix.
m <- matrix(c(1, 2, -1, 1, -1, -2, 3, 5, 4), nrow = 3)
#
print('Original matrix:')## [1] "Original matrix:"print(m)## [,1] [,2] [,3]
## [1,] 1 1 3
## [2,] 2 -1 5
## [3,] -1 -2 4#
inv_m <- myinverse(m)
print('inverse_m:')## [1] "inverse_m:"print(inv_m)## [,1] [,2] [,3]
## [1,] -0.2727273 0.45454545 -0.36363636
## [2,] 0.5909091 -0.31818182 -0.04545455
## [3,] 0.2272727 -0.04545455 0.13636364#
product <- m %*% inv_m
print("product of m and its inverse:")## [1] "product of m and its inverse:"product## [,1] [,2] [,3]
## [1,] 1 0.000000e+00 0
## [2,] 0 1.000000e+00 0
## [3,] 0 -2.775558e-17 1Observe that the diagonal elements of the product are all ones. But, the other elements have some non-zero elements, which are almost 0. So, the product is an identity matrix.
Marker: 605-04