Inverse of a Matrix / Single Value Decomposition

Problem Set #1

Given a 3 × 2 matrix A: $A = \begin{bmatrix} 1 & 2 & 3 \\ -1 & 0 & 4 \end{bmatrix}$

write code in R to compute $$\mathbf{X} = \mathbf{AA^T}$$ and $$\mathbf{Y} = \mathbf{A^TA}$$. Then, compute the eigenvalues and eigenvectors of $$\mathbf{X}$$ and $$\mathbf{Y}$$ using the built-in commands in R.

A = matrix( c( 1,-1,2,0,3,4 ), nrow = 2, ncol = 3 )
print( A )
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4
X = A %*% t( A )
Y = t( A ) %*% A
print( X )
##      [,1] [,2]
## [1,]   14   11
## [2,]   11   17
print( Y )
##      [,1] [,2] [,3]
## [1,]    2    2   -1
## [2,]    2    4    6
## [3,]   -1    6   25

Compute the eigenvalues and eigenvectors of X and Y using the built-in commans in R.

Xeig <- eigen( X )
Yeig <- eigen( Y )
print( 'eigenvalues of X:' )
## [1] "eigenvalues of X:"
print( Xeig['values'] )
## $values ## [1] 26.601802 4.398198 print( 'eigenvalues of Y:' ) ## [1] "eigenvalues of Y:" print( Yeig['values'] ) ##$values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
print( 'eigenvectors of X:' )
## [1] "eigenvectors of X:"
print( Xeig['vectors'] )
## $vectors ## [,1] [,2] ## [1,] 0.6576043 -0.7533635 ## [2,] 0.7533635 0.6576043 print( 'eigenvectors of Y:' ) ## [1] "eigenvectors of Y:" print( Yeig['vectors'] ) ##$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

Compute the left-singular, singular values, and right-singular vectors of A using the svd command:

components = svd( A )
print( 'singular values of A' )
## [1] "singular values of A"
print( components['d'] )
## $d ## [1] 5.157693 2.097188 print( 'left singular vectors of A' ) ## [1] "left singular vectors of A" print( components['u'] ) ##$u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
print( 'right singular vectors of A' )
## [1] "right singular vectors of A"
print( components['v'] )
## \$v
##             [,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:

#will use the asrenal library's comparedf
library( arsenal )
#recast as dataframes and use comparedf
#get the summary and assign to a variables
cmp_Xu <- summary( comparedf( as.data.frame( Xeig['vectors'] ), as.data.frame( components['u'] ) ) )
cmp_Yv <- summary( comparedf( as.data.frame( Yeig['vectors'] ), as.data.frame( components['v'] ) ) )
#we are only interested in the differences between the values in the dataframes
numdiffs_Xu <- n.diffs( cmp_Xu ) #determines how many differences were found between values
numdiffs_Yv <- n.diffs( cmp_Yv )
paste( 'Num differences between the eigenvectors of X and the left singular vectors of A: ', numdiffs_Xu)
## [1] "Num differences between the eigenvectors of X and the left singular vectors of A:  0"
paste( 'Num differences between the eigenvectors of Y and the right singular vectors of A: ', numdiffs_Yv)
## [1] "Num differences between the eigenvectors of Y and the right singular vectors of A:  0"
#ref 1 for more info

#comparing the eigenvalues of X, Y and the singular values of A
numdigits <- 10
test1 <- round( Xeig['values'][[1]][1], numdigits ) == round( Yeig['values'][[1]][1], numdigits )
test2 <- round( Xeig['values'][[1]][2], numdigits ) == round( Yeig['values'][[1]][2], numdigits )

#comparing X to the squares of A
test3 <- round( Xeig['values'][[1]][1], numdigits ) == round( components['d'][[1]][1]^2, numdigits )
test4 <- round( Xeig['values'][[1]][2], numdigits ) == round( components['d'][[1]][2]^2, numdigits )

if ( test1 == TRUE & test2 == TRUE )
print( 'The 2 nonzero eigenvalues of X & Y are equal' )
## [1] "The 2 nonzero eigenvalues of X & Y are equal"
if ( test3 == TRUE & test4 == TRUE )
paste( 'The 2 eigenvalues of X are equal to the square of the singular values of A (it follows that this holds for the non zero eigenvalues of Y as well)' )
## [1] "The 2 eigenvalues of X are equal to the square of the singular values of A (it follows that this holds for the non zero eigenvalues of Y as well)"

Problem Set #2

write a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors.

The Inverse Matrix will be calculated with these steps (ref#2):

1. calculate the matrix of minors
2. find the matrix of cofactors
3. find the adjugate
4. multiply by 1/Determinant
#functionalize finding the co-factor
findCofactor = function( M ){
cofactor_M <- M
stopifnot( length( unique( dim( M ) ) )==1 ) #make sure M is square
sizeM <- dim( M )[1]
for(i in 1:sizeM ){
for( j in 1:sizeM ){
cofactor_M[i,j] <- ( det(M[-i,-j]) * (-1)^(i+j) )
#STEP1 the first term (det(M[-i,-j])) finds the matrix of minors
#STEP2 the second term ((-1)^(i+j)) alternates +/- for matrix of cofactors
}
}
return( cofactor_M )
}

myinverse = function(M){
cofactor <- findCofactor( M )
T_cofactor <- t( cofactor ) #STEP3 take the transpose of the cofactor to find the adjugate (adjoint)
det_M <- det(M)
M_inverse <- T_cofactor/det_M #STEP4 multiply by 1/det to find the inverse
return( M_inverse )
}
C <- matrix(c(1,5,5,3,3,1,7,5,4), nrow=3, ncol=3, byrow=TRUE)
invC <- myinverse( C )
print( C )
##      [,1] [,2] [,3]
## [1,]    1    5    5
## [2,]    3    3    1
## [3,]    7    5    4
print( invC )
##            [,1]       [,2]       [,3]
## [1,] -0.1458333 -0.1041667  0.2083333
## [2,]  0.1041667  0.6458333 -0.2916667
## [3,]  0.1250000 -0.6250000  0.2500000
#damage control
damage_control <- C %*% invC
print( round(damage_control, digits = 14 ) )
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1