DATA 605 Assignment 4

Problem Set 1

Given a 3 × 2 matrix A, A = [1 2 3, -1 0 4], (1) write code in R to compute X = AAT and Y = ATA.

First we will initialize matrix A, to ensure it’s been created in the proper form. Then we’ll compute and display X and Y by the equations noted above.

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

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

#Compute X = AAT
#3x2 matrix multiplied by a 2x3 matrix results in 3x3 matrix
X = A %*% t(A) 
X

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

#Compute Y= ATA
#2x3 matrix multiplied by a 3x2 matrix results in 2x2 matrix
Y = t(A) %*% A 
Y

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

Compute the eigenvalues and eigenvectors of X and Y using the built-in commands in R.

The two non-zero eigenvalues (the 3rd value will be very close to zero, if not zero) of both X and Y are the same and are squares of the non-zero singular values of A (to be computed later).

Using the built in eigen() function, we compute the eigenvalues and eigenvector of X and Y, and return the corresponding values and vectors.

#Compute and display the eigenvalues and eigenvectors of X and Y.
eX <- eigen(X)
eY <- eigen(Y)

eX$values

## [1] 26.601802  4.398198

eY$values

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

eX$vectors

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

eY$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

From the output above we see that the 2 non-zero eigenvalues for X and Y are the same, while the 3rd eigenvalue for Y approaches 0.

Compute the left-singular, singular values, and right-singular vectors of A using the svd command.

SVD computes the singular-value decomposition of a rectangular matrix (in our case A). A = U D VT is the factorization of A into the product of 3 matrices described as follows:

D ($d): a vector containing singular values of A sorted decreasingly

U ($u): a matrix whose columns contain the left singular values of A

VT ($v): a matrix whose columns contain the right singular vectors of x.

#Factorize A for product matrices using svd() function
svd_A <- svd(A)
svd_A

## $d
## [1] 5.157693 2.097188
## 
## $u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
## 
## $v
##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

round(eX$values, digits = 3)

## [1] 26.602  4.398

round((svd_A$d)^2, digits = 3)

## [1] 26.602  4.398

We can confirm above that (after rounding to same decimal point), X and Y (which have the same 1st 2 values) are equivalent to the squares of the non-zero singular values of A.

Next we will examine / compare eigenvectors of X and Y with the eigenvectors of the singular vectors.

Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y.

#First we will compare the eigenvectors of X to those of U.
eX$vectors

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

svd_A$u

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

#The 1st vector (column 1) has different signs so we invert the sign of this eigenvector - the result is still an eigenvector.

eX$vectors[,1] <- eX$vectors[,1] * -1

#We compare again ...

round(eX$vectors, digits = 3)

##        [,1]   [,2]
## [1,] -0.658 -0.753
## [2,] -0.753  0.658

round(svd_A$u, digits = 3)

##        [,1]   [,2]
## [1,] -0.658 -0.753
## [2,] -0.753  0.658

Voila! We confirm above that (after rounding to same decimal point), X and U have equivalent eigenvectors.

Now we do the same for Y and V.

#First we will compare the eigenvectors of Y to those of V.
eY$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

svd_A$v

##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

#The same case as before - we invert the 1st column.

eY$vectors[,1] <- eY$vectors[,1] * -1

#We compare again, this time discounting the 3rd column (associated with a 0 eigenvalue).

round(eY$vectors[,1:2], digits = 3)

##        [,1]   [,2]
## [1,]  0.019 -0.673
## [2,] -0.255 -0.718
## [3,] -0.967  0.177

round(svd_A$v, digits = 3)

##        [,1]   [,2]
## [1,]  0.019 -0.673
## [2,] -0.255 -0.718
## [3,] -0.967  0.177

Thus confirming that (after rounding to same decimal point), Y and V have equivalent eigenvectors.

Problem Set 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.

To compute the inverse of a well-conditioned full-rank matrix, we’ll initalize the matrix of cofactors, flag for non full rank inputs, construct the matrix of cofactors (via the procedure referenced above), initialize and output matrix B as the inverse of matrix A.

myinverse = function(A){
  
  #Initialize matrix of cofactors
  coMatrix = A * 0
  
  #Flag non full rank matrices:
  if(det(A) == 0){ 
    return('Please enter a full rank, square matrix.') 
    }
  
  #Construct matrix of cofactors:
  for(i in 1:ncol(A)) {
    for(j in 1:nrow(A)) {
      coMatrix[i,j] = det(A[-i,-j]) * (-1)^(i+j) 
    }}
  
  #Initialize inverse of B as inv(A) = CT / det(A)
  B = t(coMatrix / det(A))
  
  return(B)
}

Now that we have our function built, we can test for (1) a full rank, square matrix, (2) a non full rank matrix, and (3) a non square matrix.

Your function should have the following signature: B = myinverse(A) where A is a matrix and B is its inverse and A×B = I. The off-diagonal elements of I should be close to zero, if not zero. Likewise, the diagonal elements should be close to 1, if not 1. Small numerical precision errors are acceptable but the function myinverse should be correct and must use co-factors and determinant of A to compute the inverse.

#Test 1: a full rank, square 3x3 matrix
A = matrix(c(1,0,2,2,1,0,3,2,1), nrow=  3, ncol = 3, byrow = TRUE)
A

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

B = myinverse(A)
round(B, digits = 3)

##        [,1]   [,2]   [,3]
## [1,]  0.333  1.333 -0.667
## [2,] -0.667 -1.667  1.333
## [3,]  0.333 -0.667  0.333

#confirm that the matrix product of A and B is indeed the identity matrix
round(A %*% B, digits = 3)

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

#Test 2: a non full rank matrix
A = matrix(c(0,1,2,1,2,1,2,7,8), nrow = 3, ncol = 3, byrow = TRUE)
A

##      [,1] [,2] [,3]
## [1,]    0    1    2
## [2,]    1    2    1
## [3,]    2    7    8

B=myinverse(A)
B

## [1] "Please enter a full rank, square matrix."

#Test 3: a non square matrix ---COMMENTED OUT---
#A = matrix(c(1,2,3,2,1,1), nrow = 2, ncol = 3, byrow = TRUE)
#A

#B=myinverse(A)
#B

The full rank, square matrix returned a perfect matrix inverse which was confirmed via the product of A and B as the identity matrix.

The non full rank matrix was flagged as designed and the non-square matrix is automatically flagged (as an error message) because of the determinant function.

Test matrices were pulled from: https://stattrek.com/matrix-algebra/matrix-rank.aspx