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

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


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

eigen(Y)
## $values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
## 
## $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

X =A x t(A) = U x D1 X t(U)

We know U is the eignevector of A x t(A) which is 2 x 2 matrix Lets calculate eigen of X

eigen(X)
## $values
## [1] 26.601802  4.398198
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043

We see that eigen vector of X is symmetric orthonormal matrix.

U <- matrix(c(0.6576043,0.7533635,-0.7533635, 0.6576043), nrow=2)

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

U is the left singular vector

We should also note that eigen values of X, or A x t(A), this eigen values make up a 2 x 2 diagonal matrix D1, A x t(A) = U x D1 X t(U)

D1 <- matrix(c(26.601802,0,0,4.398198), nrow =2)
D1
##         [,1]     [,2]
## [1,] 26.6018 0.000000
## [2,]  0.0000 4.398198

For V,

Y =t(A) x A = V x D X t(V)

We know V is the eignevector of t(A) x A which is 3 x 3 matrix Lets calculate eigen of Y

eigen(Y)
## $values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
## 
## $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

We see that eigen vector of Y is symmetric orthonormal matrix, V.

Lets look at the inverse of Eigen (Y)

V <- eigen(Y)$vectors

solve(V)
##             [,1]       [,2]      [,3]
## [1,] -0.01856629  0.2549994 0.9667630
## [2,] -0.67279027 -0.7184510 0.1765824
## [3,]  0.73960026 -0.6471502 0.1849001

We also notice it just the transpose of V

We see that t(V) = solve(V) #inverse of V

t(V) is the right singular vector

We should also note that eigen values of Y, or t(A) x A, this eigen values make up a 3 x 3 diagonal matrix D, A x t(A) = V x D X t(V)

D <- matrix(c(2.660180e+01,0,0,0, 4.398198e+00,0,0,0,1.058982e-16), byrow = TRUE, nrow =3)
D
##         [,1]     [,2]         [,3]
## [1,] 26.6018 0.000000 0.000000e+00
## [2,]  0.0000 4.398198 0.000000e+00
## [3,]  0.0000 0.000000 1.058982e-16

We notice that D consists of eigen values that are similar to D1 expect that one value that is close to zero.

The singular values of S is 2x3 matrix where lamda(ii) is sqrt(D)

S <- matrix(c(sqrt(26.6018),0,0,0,sqrt(4.398198),0), nrow=2, byrow = TRUE)
S
##          [,1]     [,2] [,3]
## [1,] 5.157693 0.000000    0
## [2,] 0.000000 2.097188    0

Lets confirm this SVD decomposition by multiplying the factors A = U x S x t(V)

A_svd <- U %*% S %*% t(V)
A_svd
##            [,1]          [,2] [,3]
## [1,]  0.9999999  2.000000e+00    3
## [2,] -1.0000000 -4.661433e-08    4
#compare to
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4

notice that the product of U x S x t(V) gives a value approximately A

Using the svd function we get similar values execpt that U and V values are switched on mine.

svd(A,2,3)
## $d
## [1] 5.157693 2.097188
## 
## $u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
## 
## $v
##             [,1]       [,2]       [,3]
## [1,]  0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510  0.6471502
## [3,] -0.96676296  0.1765824 -0.1849001
B <-svd(A,2,3)

Computing the product of factors from the svd command gives A= U x S x t(V)

d <- S
A_compute_Svd <- B$u %*% d %*% t(B$v)
A_compute_Svd
##      [,1]         [,2] [,3]
## [1,]    1 2.000000e+00    3
## [2,]   -1 7.962203e-09    4
#compare both
A_svd
##            [,1]          [,2] [,3]
## [1,]  0.9999999  2.000000e+00    3
## [2,] -1.0000000 -4.661433e-08    4

We can compare this values gotten from svd command (A_compute_Svd) factorization to the one A_svd computed earlier and we notice, it is very similar and all approximately equal to A

PRoblem 2

Using the procedure outlined in section 1 of the weekly handout, 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)

Solution

We know inv(A) = t(C)/det(A) where C is the cofactor of the matrix

Cij= (-1)^(i+j)det(Mij)

where M is the matrix remaining after removing the ith row and jth column of Matrix A

myinverse <- function(A)
{
  Rownum <- nrow(A)
  Colnum <- ncol(A)
  
  #initialize cofactor matrix C
  
  C <- diag(Rownum)
  
  #calculate Cofactor matrix C
  
  for(i in 1:Rownum )
  {
    for(j in 1:Colnum)
    {
      p <- i+j
      C[i,j] <- ((-1)^p)*det(as.matrix(A[-i,-j])) 
    }
  }

  A_inv <- t(C)/det(A)
  return(A_inv)
}

#as.matrix(A[-1,-1]) %*% diag(2)

#Test function 
A<- matrix(c(1,2,4,-1,-3,4,1,8,1), nrow=3, byrow=TRUE)
#A<- matrix(c(1,2,3,4), nrow=2, byrow=TRUE)
B= myinverse(A)
B
##            [,1]        [,2]        [,3]
## [1,]  0.7777778 -0.66666667 -0.44444444
## [2,] -0.1111111  0.06666667  0.17777778
## [3,]  0.1111111  0.13333333  0.02222222
#compare with solve function
solve(A)
##            [,1]        [,2]        [,3]
## [1,]  0.7777778 -0.66666667 -0.44444444
## [2,] -0.1111111  0.06666667  0.17777778
## [3,]  0.1111111  0.13333333  0.02222222
all.equal(B,solve(A))
## [1] TRUE