DATA 605 - Assignment 4

In this problem, we’ll verify using R that SVD and Eigenvalues are related as worked out in the weekly module. Given a 3 × 2 matrix A:

$$\mathbf{A} = \left[\begin{array} {rrr} 1 & 2 & 3 \\ -1 & 0 & 4 \end{array}\right]$$

A <- matrix(c(1,-1,2,0,3,4), nrow=2, byrow=F)
A

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4

Write code in R to compute $X = AA^{T} ; Y = A^{T}A$.

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 the eigenvalues and eigenvectors of X and Y using the built-in commands in R.

eigenvalueX<-eigen(X)$values
eigenvectorX<-eigen(X)$vectors
eigenvalueY<-eigen(Y)$values
eigenvectorY<-eigen(Y)$vectors
eigenvalueX

## [1] 26.601802  4.398198

eigenvectorX

##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043

eigenvalueY

## [1] 2.660180e+01 4.398198e+00 1.058982e-16

eigenvectorY

##             [,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.

The columns of U are the left singular vectors and the columns of V are the right singular vectors. The matrix A is a diagonal matrix, whose diagonal values are the singular values of the matrix X.

singularValues<-svd(A)$d
leftSingularVectors<-svd(A)$u
rightSingularVectors<-svd(A)$v
singularValues

## [1] 5.157693 2.097188

leftSingularVectors

##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043

rightSingularVectors

##             [,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.

The columns of V / left singular vectors are the eigenvectors of X. The negative value represents a scalar.

leftSingularVectors

##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043

eigenvectorX

##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043

The columns of U /right singular vectors are the eigenvectors of Y. The negative value represents a scalar. The third near-zero number is negligible.

rightSingularVectors

##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

eigenvectorY[,1:2]

##             [,1]       [,2]
## [1,] -0.01856629 -0.6727903
## [2,]  0.25499937 -0.7184510
## [3,]  0.96676296  0.1765824

In addition, the two non-zero eigenvalues of both X and Y are the same and are squares of the non-zero singular values of A.

all.equal(eigenvalueX,(singularValues^2))

## [1] TRUE

all.equal(eigenvalueY[1:2],(singularValues^2))

## [1] TRUE

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)$

where A is a matrix and B is its inverse and A×B = I. The off-diagonal elements of I should be close to zero, if not zero. Likewise, the diagonal elements should be close to 1, if not 1. Small numerical precision errors are acceptable but the function myinverse should be correct and must use co-factors and determinant of A to compute the inverse.

B <- matrix(c(1,-1,2,0), nrow=2)#2x2 square
C <- matrix(c(1,-1,2,0,3,4,6,7,8), nrow=3)#3x3 square
D <- matrix(c(1,-1,2,0,3,4,6,7,8,5,6,7), nrow=3)#3x3 nonsquare test

#function for determining the inverse of a matrix
myinverse<-function(A){
  n<-ncol(A)
  if(n==nrow(A)){
      #square test
    if(det(A)!=0){
      #determinant is not 0
      cofactors<- matrix(0:0,n,n)
      for(i in 1:n){
        for(j in 1:n){
          detVal<-det(A[-i,-j])
          cofactors[i,j]<- detVal*(-1^(i+j))
        }
      }
      #inverse calc cofactor matrix transpose divided by the determinant
      return(t(cofactors)/det(A))
    }else{
      stop('Determinant is 0. Break.')
    }
  }else{
    stop('Non-square matrix. Break.')
  }
}

#result
myinverse(C)

##          [,1]    [,2]      [,3]
## [1,] -0.06250 -0.3750 -0.281250
## [2,] -0.34375 -0.0625  0.203125
## [3,] -0.15625  0.0625  0.046875

#sanity test
solve(C)

##          [,1]    [,2]      [,3]
## [1,]  0.06250 -0.3750  0.281250
## [2,] -0.34375  0.0625  0.203125
## [3,]  0.15625  0.0625 -0.046875