Problem Set 1

Singular Value Decomposition: The goal of “diagonalizing” a matrix is to find an orthonormal basis that gets transformed to another orthonormal basis (times some scalar) after applying the linear transformation (i.e., multiplying by the matrix.)

In other words, if we have a matrix A, we are looking for A V = U D where V and U are matrices of orthonormal basis vectors and D is a diagonal matrix of “scalars.”

Rather than trying to find V and U simultaneously, we can use algebra to arrive at the following formula:
AT A = V D2 VT
where V is the eigenvectors and D2 is the eigenvalues of the matrix AT A (and we can then just take the square root of D2 to get the D scalar values we were originally looking for.)
Likewise, A AT = U D2 UT

# To test, we start with a matrix A
A <- rbind(c( 1, 2, 3),
           c(-1, 0, 4))

# Then we create AA^t and A^tA
X <- A %*% t(A)

Y <- t(A) %*% A

# Using R's eigen() function we can get the eigenvectors and eigenvalues of X and Y
X_eig <- eigen(X)

Y_eig <- eigen(Y)

# And we compare the output to the svd() function applied to our original matrix A
A_svd <- svd(A)

# The square root of the eigenvalues of X (and Y too) should be the same as the D in our svd of A:
sqrt(X_eig$values)
## [1] 5.157693 2.097188
sqrt(Y_eig$values)
## [1] 5.157693e+00 2.097188e+00 1.029068e-08
A_svd$d
## [1] 5.157693 2.097188
# The eigenvectors of X should be the same (or same span) as U
X_eig$vectors
##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043
A_svd$u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
# The eigenvectors of Y should be the same (or same span) as V
Y_eig$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
A_svd$v
##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

Problem Set 2

Create a function that computes the inverse of a matrix using co-factors.

myinverse <- function(A){
  C <- array(0, dim=dim(A))
  for (i in 1:nrow(A)){
    for (j in 1:ncol(A)){
      C[i,j] <- (-1)^(i+j)*det(A[-i,-j])
    }
  }
  t(C)/det(A)
}

# Test the function
A <- rbind(c(1,2,3),
           c(7,0,4),
           c(4,9,0))

B <- myinverse(A)

# Compare to solve() function
B
##             [,1]         [,2]        [,3]
## [1,] -0.19459459  0.145945946  0.04324324
## [2,]  0.08648649 -0.064864865  0.09189189
## [3,]  0.34054054 -0.005405405 -0.07567568
solve(A)
##             [,1]         [,2]        [,3]
## [1,] -0.19459459  0.145945946  0.04324324
## [2,]  0.08648649 -0.064864865  0.09189189
## [3,]  0.34054054 -0.005405405 -0.07567568
# AB=I
round(A %*% B, 2)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1