Assignment 4

Fundamentals of Computational Mathematics

CUNY MSDS DATA 605, Fall 2018

Rose Koh

09/21/2018

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

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

# Compute X, Y
X <- A %*% t(A)
X

##      [,1] [,2]
## [1,]   14   11
## [2,]   11   17

Y <- t(A) %*% A
Y

##      [,1] [,2] [,3]
## [1,]    2    2   -1
## [2,]    2    4    6
## [3,]   -1    6   25

# Compute Eigenvectors
X.eigenvalues <- eigen(X)$values
X.eigenvectors <- eigen(X)$vectors
Y.eigenvalues <- eigen(Y)$values
Y.eigenvectors <- eigen(Y)$vectors
X.eigenvalues

## [1] 26.601802  4.398198

X.eigenvectors

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

Y.eigenvalues

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

Y.eigenvectors

##             [,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(u), singular values(d), right singular(v) vectors of A
svd.A <- svd(A)
left.sing <- svd.A$u
right.sing <- svd.A$v
sing.values <- svd.A$d

# Examine 2 sets of singular vectors and show they are eigenvectors of X and Y
compare.ux <- cbind(abs(left.sing),abs(X.eigenvectors))
colnames(compare.ux) <- c("SVD.u1", "SVD.u2", "X.eigenvectors=SVD.u1", "X.eigenvectors=SVD.u2")
compare.ux

##         SVD.u1    SVD.u2 X.eigenvectors=SVD.u1 X.eigenvectors=SVD.u2
## [1,] 0.6576043 0.7533635             0.6576043             0.7533635
## [2,] 0.7533635 0.6576043             0.7533635             0.6576043

compare.vy <- cbind(abs(right.sing),abs(Y.eigenvectors[,1:2]))
colnames(compare.vy) <- c("SVD.v1", "SVD.v2", "Y.eigenvectors1=SVD.v1", "Y.eigenvectors2=SVD.v2")
compare.vy

##          SVD.v1    SVD.v2 Y.eigenvectors1=SVD.v1 Y.eigenvectors2=SVD.v2
## [1,] 0.01856629 0.6727903             0.01856629              0.6727903
## [2,] 0.25499937 0.7184510             0.25499937              0.7184510
## [3,] 0.96676296 0.1765824             0.96676296              0.1765824

# two non-zero eigenvalues of both X and Y are the same and are squares of the non-zero singular values of A
sing.values

## [1] 5.157693 2.097188

sqrt(X.eigenvalues)

## [1] 5.157693 2.097188

We can clearly see that the singular vectors calculated using SVD function returns the same value as eigenvectors of X and Y. the vectors are decared in unit vector form, thus the absolute value was used to compare them.

myinverse <- function(A) {
  # Initialize co.factor matrix  
  co.factor <- matrix(c(0:0), nrow=nrow(A), ncol=ncol(A), byrow = TRUE)
  print("====== Initial Co-factors matrix ======")
  print(co.factor)
  # Loop through every elements of the matrix
  for (i in 1:nrow(A)) {
    for (j in 1:ncol(A)){
      # Eliminate the row and the column of current element when calculating the co-factor
      co.factor[i,j] <- (-1)^(i+j)*det(A[-i,-j])
    }
  }
  print("====== Co-factors matrix ======")
  print(co.factor)
  
  # Calculate inverse of matrix A: divide A transpose by det A
  A.inv <- t(co.factor)/det(A)
  print("====== Inverse Matrix ======")
  print(A.inv)
  
  # Identity Matrix
  print("====== Identity matrix ======")
  round(A.inv %*% A, 0)
}

test <- matrix(sample.int(15, size = 9*100, replace = TRUE), nrow = 5, ncol = 5)
myinverse(test)

## [1] "====== Initial Co-factors matrix ======"
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    0    0    0    0
## [4,]    0    0    0    0    0
## [5,]    0    0    0    0    0
## [1] "====== Co-factors matrix ======"
##       [,1]  [,2]  [,3]  [,4]  [,5]
## [1,] -4920  3428  1890  5144 -2276
## [2,]  4525  5573 -2660 -5960     4
## [3,] -4707 -7855  3836  3594  6450
## [4,]  7063  4837 -8428  1120 -1624
## [5,] -1468 -9104 10150 -1808   894
## [1] "====== Inverse Matrix ======"
##             [,1]          [,2]        [,3]        [,4]        [,5]
## [1,] -0.08136804  7.483544e-02 -0.07784540  0.11680945 -0.02427811
## [2,]  0.05669302  9.216750e-02 -0.12990772  0.07999537 -0.15056395
## [3,]  0.03125724 -4.399166e-02  0.06344061 -0.13938412  0.16786293
## [4,]  0.08507260 -9.856779e-02  0.05943836  0.01852281 -0.02990110
## [5,] -0.03764099  6.615288e-05  0.10667152 -0.02685807  0.01478517
## [1] "====== Identity matrix ======"

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

# same as inverse matrix
solve(test)

##             [,1]          [,2]        [,3]        [,4]        [,5]
## [1,] -0.08136804  7.483544e-02 -0.07784540  0.11680945 -0.02427811
## [2,]  0.05669302  9.216750e-02 -0.12990772  0.07999537 -0.15056395
## [3,]  0.03125724 -4.399166e-02  0.06344061 -0.13938412  0.16786293
## [4,]  0.08507260 -9.856779e-02  0.05943836  0.01852281 -0.02990110
## [5,] -0.03764099  6.615288e-05  0.10667152 -0.02685807  0.01478517