Problem Set 1

Verify using R that SVD and Eigenvalues are related.

A. Given a 3 X 2 matrix A, write code in R to compute \(X = AA^T\) and \(Y = A^TA\). \[A = \begin{bmatrix}1 & 2 & 3 \\ -1 & 0 & 4\end{bmatrix}\]

# Define matrix A.
A <- matrix(c(1, 2, 3, -1, 0, 4), 2, 3)

# Calculate X and Y.
X <- A %*% t(A)
Y <- t(A) %*% A

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

# Compute the eigenvalues and eigenvectors of X and Y.
eigens_x <- eigen(X)
eigens_y <- eigen(Y)

Eigenvalues of X.

eigens_x$values
## [1] 21.09017  9.90983

Eigenvectors of X.

eigens_x$vectors
##            [,1]       [,2]
## [1,] -0.0898056 -0.9959593
## [2,]  0.9959593 -0.0898056

Eigenvalues of Y.

eigens_y$values
## [1] 21.09017  9.90983  0.00000

Eigenvectors of Y.

eigens_y$vectors
##            [,1]       [,2]       [,3]
## [1,] -0.4141868 -0.3734355  0.8300574
## [2,]  0.2755368 -0.9206109 -0.2766858
## [3,] -0.8674842 -0.1141117 -0.4842001

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

# Compute the left-singular vector of A.
left_singular_vectors <- svd(A)$u

# Compute the singular values of A.
singular_values <- svd(A)$d

# Compute the right-singular vector of A.
right_singular_vectors <- svd(A)$v

Left singular vectors of A.

left_singular_vectors
##            [,1]      [,2]
## [1,] -0.0898056 0.9959593
## [2,]  0.9959593 0.0898056

Singular values of A.

singular_values
## [1] 4.592404 3.147988

Right singular vectors of A.

right_singular_vectors
##            [,1]      [,2]
## [1,]  0.4141868 0.3734355
## [2,] -0.2755368 0.9206109
## [3,]  0.8674842 0.1141117

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

# Compare the left singular vectors of A with the eigenvectors of X.
left_singular_vectors
##            [,1]      [,2]
## [1,] -0.0898056 0.9959593
## [2,]  0.9959593 0.0898056
eigens_x$vectors
##            [,1]       [,2]
## [1,] -0.0898056 -0.9959593
## [2,]  0.9959593 -0.0898056
# Confirm that they are indeed equal.
round(abs(left_singular_vectors)) == round(abs(eigens_x$vectors))
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
# Compare the right singular vectors of A with the eigenvectors of Y.
right_singular_vectors
##            [,1]      [,2]
## [1,]  0.4141868 0.3734355
## [2,] -0.2755368 0.9206109
## [3,]  0.8674842 0.1141117
eigens_y$vectors
##            [,1]       [,2]       [,3]
## [1,] -0.4141868 -0.3734355  0.8300574
## [2,]  0.2755368 -0.9206109 -0.2766858
## [3,] -0.8674842 -0.1141117 -0.4842001

E. Show that the two non-zero eigenvalues of both X and Y are the same.

# Compare the two non-zero eigenvalues of both X and Y to confirm that they the same.
eigens_y$values
## [1] 21.09017  9.90983  0.00000
eigens_x$values
## [1] 21.09017  9.90983

F. Show that the two non-zero eigenvalues of both X and Y are squares of the non-zero singular values of A.

# Compare the two non-zero eigenvalues of both X and Y to confirm that they the same.
eigens_x$values
## [1] 21.09017  9.90983
# Eigenvalues of Y.
eigens_y$values
## [1] 21.09017  9.90983  0.00000
# Singular values of A.
singular_values^2
## [1] 21.09017  9.90983
# Confirm that the eigenvalues of X are squares of the non-zero singular values of A.
round(abs(eigens_x$values[1:2])) == round(abs(singular_values^2))
## [1] TRUE TRUE
# Confirm that the eigenvalues of Y are squares of the non-zero singular values of A.
round(abs(eigens_y$values[1:2])) == round(abs(singular_values^2))
## [1] TRUE TRUE

Problem Set 2

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

Your function should have the following signature: B = myinverse(A) where A is a matrix and B is its inverse and AxB = 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.

# Define the function myinverse().
myinverse <- function(A) {
  # Check if the passed matrix is a square matrix.
  if (nrow(A) != ncol(A)) {
    return('The provided matrix is not a square matrix.')
  }

  # Check to ensure that the passed matrix is invertable.
  if (det(A) == 0) {
    return('The provided matrix is not invertable.')
  }

  # Assign the matrix rows and columns to variables (rows and columns).
  rows <- nrow(A)
  columns <- ncol(A)

  # Define the co-factor matrix.
  cofactor_matrix <- diag(0, rows, columns)
  for(i in 1:rows) {
    for(j in 1:columns) {
      # Populate the matrix with rows and columns.
      cofactor_matrix[i,j] <- (-1) ^ (i + j) * det(A[-i, -j])
    }
  }

  # Define the inverted matrix and return it.
  inverted_matrix <- (t(cofactor_matrix) / det(A))

  return(inverted_matrix)
}

Test the myinverse() function with a 3x3 matrix.

# Define the 3x3 matrix A.
A <- matrix(c(1, 0, 5, 3, 1, 5, 3, 4, 0), 3, 3)
A
##      [,1] [,2] [,3]
## [1,]    1    3    3
## [2,]    0    1    4
## [3,]    5    5    0
# Pass the matrix to the myinverse() function and print the inverted matrix.
B <- myinverse(A)
B
##      [,1] [,2]  [,3]
## [1,] -0.8  0.6  0.36
## [2,]  0.8 -0.6 -0.16
## [3,] -0.2  0.4  0.04
# Check that the result is correct.
accuracy_check <- solve(A)
round(accuracy_check, 2) == round(B, 2)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
# Check that the off-diagonal elements of I are zero, or close to zero,
# and that the diagonal elements are 1, or close to 1.
round(A %*% B, 2)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Test the myinverse() function with a 4x4 matrix.

# Define the 4x4 matrix A.
A <- matrix(c(1, 6,-2, -3, 2, -8, 3, -6, -1, 0, 3, 16, -5, -2, 7, 9), 4, 4)
A
##      [,1] [,2] [,3] [,4]
## [1,]    1    2   -1   -5
## [2,]    6   -8    0   -2
## [3,]   -2    3    3    7
## [4,]   -3   -6   16    9
# Pass the matrix to the myinverse() function and print the inverted matrix.
B <- myinverse(A)
B
##             [,1]       [,2]       [,3]        [,4]
## [1,]  0.25607477 0.20140187 0.28971963 -0.03831776
## [2,]  0.21682243 0.01542056 0.18691589 -0.02149533
## [3,]  0.18504673 0.01962617 0.05607477  0.06355140
## [4,] -0.09906542 0.04252336 0.12149533 -0.02897196
# Check that the result is correct.
accuracy_check <- solve(A)
round(accuracy_check, 2) == round(B, 2)
##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE
# Check that the off-diagonal elements of I are zero, or close to zero,
# and that the diagonal elements are 1, or close to 1.
round(A %*% B, 2)
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1