Problem Set 1

Verify that SVD and Eigenvalues are related. Given the 2x2 matrix A \[\mathbf{A} = \left[\begin{array} {rrr} 1 & 2 & 3 \\ -1 & 0 & 4 \end{array}\right] \] Compute \(X=AA^T\) and \(Y=A^TA\). Compute eigenvalues and eigenvectors of X and Y (using built-in R commands). Compute the left-singular, singular values and right-singular vectors of A (svd command). Show the two sets of singular vectors are eigenvectors of X and Y. Store scalars and vectors in variables.

  1. Find \(X=AA^T\).
A = matrix(c(1, 2, 3, -1, 0, 4), nrow = 2, byrow = T)
X = A %*% t(A)
X
##      [,1] [,2]
## [1,]   14   11
## [2,]   11   17
  1. Find \(Y=A^TA\).
Y = t(A) %*% A
Y
##      [,1] [,2] [,3]
## [1,]    2    2   -1
## [2,]    2    4    6
## [3,]   -1    6   25
  1. Eigenvalues/Eigenvectors for X and Y.
eigens = function(A){
  Val = eigen(A)$values #eigenvalues
  Vec = eigen(A)$vectors #eigenvectors
  return(list(round(Val,4),round(Vec,4)))
}

Xeigens = eigens(X)
Yeigens = eigens(Y)
Xeigens
## [[1]]
## [1] 26.6018  4.3982
## 
## [[2]]
##        [,1]    [,2]
## [1,] 0.6576 -0.7534
## [2,] 0.7534  0.6576
Yeigens
## [[1]]
## [1] 26.6018  4.3982  0.0000
## 
## [[2]]
##         [,1]    [,2]    [,3]
## [1,] -0.0186 -0.6728  0.7396
## [2,]  0.2550 -0.7185 -0.6472
## [3,]  0.9668  0.1766  0.1849
  1. Computer the left-singular, singular values and right-singular vectors of A using the svd command.
svd.left = round(svd(A)$u,4)
svd.left
##         [,1]    [,2]
## [1,] -0.6576 -0.7534
## [2,] -0.7534  0.6576
svd.sing =  round(svd(A)$d,4)
svd.sing
## [1] 5.1577 2.0972
svd.right =  round(svd(A)$v,4)
svd.right
##         [,1]    [,2]
## [1,]  0.0186 -0.6728
## [2,] -0.2550 -0.7185
## [3,] -0.9668  0.1766
  1. Show the two sets of singular vectors are eigenvectors of X and Y.
# We can see that by comparing the X eigenvectors and the left-singular values, the are multiples of each other
list(Xeigens[2], svd.left)
## [[1]]
## [[1]][[1]]
##        [,1]    [,2]
## [1,] 0.6576 -0.7534
## [2,] 0.7534  0.6576
## 
## 
## [[2]]
##         [,1]    [,2]
## [1,] -0.6576 -0.7534
## [2,] -0.7534  0.6576
# The same can be said for the Y eigenvectors and the right-singular values.
# The only thing with Y is that the svd.right function only provides the first two eigenvectors.
list(Yeigens[2], svd.right)
## [[1]]
## [[1]][[1]]
##         [,1]    [,2]    [,3]
## [1,] -0.0186 -0.6728  0.7396
## [2,]  0.2550 -0.7185 -0.6472
## [3,]  0.9668  0.1766  0.1849
## 
## 
## [[2]]
##         [,1]    [,2]
## [1,]  0.0186 -0.6728
## [2,] -0.2550 -0.7185
## [3,] -0.9668  0.1766
  1. Show that two non-zero eigenvalues of both X and Y are the same and are squares of the non-zero singular values of A.
list(Xeigens[1], Yeigens[1], svd.sing^2)
## [[1]]
## [[1]][[1]]
## [1] 26.6018  4.3982
## 
## 
## [[2]]
## [[2]][[1]]
## [1] 26.6018  4.3982  0.0000
## 
## 
## [[3]]
## [1] 26.601869  4.398248

Problem Set 2

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

myinverse = function(A) {
  # Check that the entry is a square matrix
  rows = nrow(A)
  cols = ncol(A)
  if(rows!=cols | qr(A)$rank!= dim(A)[1]){
    return("Nope. Must be square matrix and full-rank.")
  }else{
      # create an empty matrix for cofactors
    cofactors = matrix(rep(0, length(A)), ncol=cols, nrow = rows)
    # fill in every entry with a determinant 
    # go first by row
    for (i in 1:cols) {
      # then by column
      for (j in 1:rows) {
        # exclude the entry to calculate determinate
        dets = det(A[-i,-j])
        # alternate signs
        cofactors[i,j] = dets*((-1)^(i+j))
      }
      }
    # Transpose matrix and divide by determinant
    inverse = t(cofactors)/det(A)
    return(inverse)
  }
}
  
A = matrix(c(1,2,3,-1,0,4,5,7,9), nrow = 3, byrow = T)
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4
## [3,]    5    7    9
B = myinverse(A)
B
##            [,1]       [,2]       [,3]
## [1,] -3.1111111  0.3333333  0.8888889
## [2,]  3.2222222 -0.6666667 -0.7777778
## [3,] -0.7777778  0.3333333  0.2222222

\(AB=I\). Meaning that, to check whether our function works properly, the matrix multiplication of A and B should equal an identity matrix. To ensure that no numerical precision errors occur (due to some values being extremely small), we round all values withing the matrix.

round(B%*%A,4)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

The inverse function is successful.