Problem Set 1

library(tidyverse)
A <- matrix(c(1, -1, 2, 0, 3, 4),2,3)
X <- A %*% t(A)
Y <- t(A) %*% A

eigen_x <- eigen(X)

eigen_y <- eigen(Y)
decomp <- svd(A)

gamma <- diag(decomp[[1]], 2,2)

We can see that the non-zero values of the X and Y matrices match the square of the values from gamma:

round(eigen_x$values, 3)
## [1] 26.602  4.398
round(eigen_y$values, 3)
## [1] 26.602  4.398  0.000
round(gamma ** 2, 3)
##        [,1]  [,2]
## [1,] 26.602 0.000
## [2,]  0.000 4.398

Next, we’ll compare the eigenvectors of x and y to the matching elements of the decomposition

eigen_x$vectors
##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043
decomp$u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
eigen_y$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
decomp$v
##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

SVD is not unique, so the signs of U and V don’t match the signs of X’s and Y’s eigenvectors, respectively. However, the absolute values of each entry do match.

near(abs(eigen_x$vectors), abs(decomp$u))
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
near(abs(eigen_y$vectors)[,1:2], abs(decomp$v))
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
## [3,] TRUE TRUE

We can see that both sets of eigenvectors will work:

A_prime <- decomp[[2]] %*% gamma %*% t(decomp[[3]])
A_prime2 <- eigen_x$vectors %*% gamma %*% t(eigen_y$vectors[,1:2])


near(A, A_prime)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
near(A, A_prime2)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE

Problem Set 2

Cofactor matrix inversion function

myinverse <- function(A){
  rows <- nrow(A)
  cols <- ncol(A)
  adet <- det(A)
  C = diag(0, rows, cols)
  if (rows != cols | adet == 0){
    print("This matrix is not invertible")
  }else{
    for (i in 1:rows){
      for (j in 1:cols){
        if (rows > 2){
          C[i,j] = (-1)**(i+j) * det(A[c(1:rows) != i, c(1:cols) != j])
        }else{
          C[i,j] = (-1)**(i+j) * A[c(1:rows) != i, c(1:cols) != j]
        }
      }
    }
    Ainv = t(C)/adet
    return(Ainv)
  }
}

Test

A <- matrix(c(3,4,-3,2,1,-1,-1,6,2), 3, 3)
Ainv <- myinverse(A)
Ainv
##             [,1]       [,2]       [,3]
## [1,] -0.29629630  0.1111111 -0.4814815
## [2,]  0.96296296 -0.1111111  0.8148148
## [3,]  0.03703704  0.1111111  0.1851852
round(A %*% Ainv, 3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
A <- matrix(c(1,2,3,4), 2, 2)

Ainv <- myinverse(A)
Ainv
##      [,1] [,2]
## [1,]   -2  1.5
## [2,]    1 -0.5
round(A %*% Ainv, 3)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

The identity matrix is returned, as expected.