In this problem, we’ll verify using R that SVD and Eigenvalues are related as worked out in the weekly module. Given a 3×2 matrix \(A\)
\(A = \left( \begin{array}{ccc} 1 & 2 & 3\\ -1 & 0 & 4\end{array} \right)\)
write code in R to compute \(X = AA^T\) and \(Y = A^TA\).
# Create matrix A
A <- matrix(c(1,-1,2,0,3,4),nrow=2)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4# Calculate X and Y
X <- A %*% t(A)
Y <- t(A) %*% A
X## [,1] [,2]
## [1,] 14 11
## [2,] 11 17Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25Then, compute the eigenvalues and eigenvectors of \(X\) and \(Y\) using the built-in commands in R.
# Calculate and display eigenvalues and eigenvectors for X and Y
ex <- eigen(X)
ex_vals <- ex$values
ex_vals## [1] 26.601802 4.398198ex_vect <- ex$vectors
ex_vect## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043ey <- eigen(Y)
ey_vals <- ey$values
ey_vals## [1] 2.660180e+01 4.398198e+00 1.058982e-16ey_vect <- ey$vectors
ey_vect## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001Then, compute the left-singular, singular values, and right-singular vectors of \(A\) using the svd command.
# Use svd function to compute the singular-value decomposition of X
svd_x <- svd(X)
# The left-singular vectors of X
svd_x$u## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043# The singular values of X
svd_x$d## [1] 26.601802 4.398198# The right-singular values of X
svd_x$v## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043# Use svd function to compute the singular-value decomposition of Y
svd_y <- svd(Y)
# The left-singular vectors of X
svd_y$u## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001# The singular values of X
svd_y$d## [1] 2.660180e+01 4.398198e+00 2.836098e-15# The right-singular values of X
svd_y$v## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 -0.7396003
## [2,] 0.25499937 -0.7184510 0.6471502
## [3,] 0.96676296 0.1765824 -0.1849001Examine the two sets of singular vectors and show that they are indeed eigenvectors of \(X\) and \(Y\).
# First for X - create diagonal matrix from singular values
x_sigma <- diag(svd_x$d)
x_sigma## [,1] [,2]
## [1,] 26.6018 0.000000
## [2,] 0.0000 4.398198# Mutiply u * sigma * the conjugate transpose of v
# to verify factorization of X
svd_x$u %*% x_sigma %*% Conj(t(svd_x$v))## [,1] [,2]
## [1,] 14 11
## [2,] 11 17# Compare to X
X## [,1] [,2]
## [1,] 14 11
## [2,] 11 17# Next for Y - create diagonal matrix from singular values
y_sigma <- diag(svd_y$d)
y_sigma## [,1] [,2] [,3]
## [1,] 26.6018 0.000000 0.000000e+00
## [2,] 0.0000 4.398198 0.000000e+00
## [3,] 0.0000 0.000000 2.836098e-15# Mutiply u * sigma * the conjugate transpose of v
# to verify factorization of Y
svd_y$u %*% y_sigma %*% Conj(t(svd_y$v))## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25# Compare to Y
Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25In addition, 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\).
Your code should compute all these vectors and scalars and store them in variables. Please add enough comments in your code to show me how to interpret your steps.
# You can see by inspection that the non-zero singular values
# of X and Y are the same, not withstanding R rounding quirks
svd_x$d## [1] 26.601802 4.398198svd_y$d## [1] 2.660180e+01 4.398198e+00 2.836098e-15# We can run the svd function on A to get the non-singular values
svd_a <- svd(A)
svd_a$d## [1] 5.157693 2.097188# Squares of non-singular values of A
svd_a$d^2## [1] 26.601802 4.398198# We can then compare the squares of the non-singular values of
# A to the non-zero singular values of X and Y (using X, because X
# is a cleaner matrix) and rounding to 6 digits (to eliminate
# differences that might exists 14 or 15 digits out)
round(svd_a$d^2, 6) == round(svd_x$d, 6)## [1] TRUE TRUEUsing 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)
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.
Solution: Compute cofactor for each sub-matrix of \(A\) and insert into a cofactor matrix using formula below where \(S\) is the sub-matrix
\(C_{ij} = (-1)^{i+j}det(S)\)
Then, compute the inverse of \(A\) by taking the cofactor matrix \((C^T)\) and dividing it by the determinant of \(A\)
\(A^-1 = C^T/det(A)\)
# Function to compute the inverse of A
myinverse <- function(A) {
# check to make sure matrix is square
if (nrow(A) != ncol(A)) {
print("Your matrix is not square")
break}
# initialize cofactor matrix of NAs
cfmatrix <- matrix(nrow=nrow(A),ncol=ncol(A))
# create vectors of rows and column to loop through
seqr <- seq(nrow(A))
seqc <- seq(ncol(A))
for (i in seqr) {
for (j in seqc) {
submatrix <- A[seqr[-i],seqc[-j]]
# apply formula to find cofactor for sub-matrix
cfmatrix[i,j] = ((-1)^(i+j))*det(submatrix)
}}
inv <- t(cfmatrix)/det(A)
}
# Use function to compute inverse of sample A matrix from handout
A <- matrix(c(1,2,3,2,1,3,1,1,3),nrow=3)
A## [,1] [,2] [,3]
## [1,] 1 2 1
## [2,] 2 1 1
## [3,] 3 3 3# Use the myinverse function to compute the inverse of A
invA <- myinverse(A)
invA## [,1] [,2] [,3]
## [1,] 0 1 -0.3333333
## [2,] 1 0 -0.3333333
## [3,] -1 -1 1.0000000# Check that A times the inverse of A is the identity matrix
A %*% invA## [,1] [,2] [,3]
## [1,] 1 0 -1.110223e-16
## [2,] 0 1 -1.110223e-16
## [3,] 0 0 1.000000e+00