Data 605 HW4: Linear Transformations, Representations
Please refer to the Assignment 4 Document.
1 Problem Set 1
In this problem, we’ll verify using R that SVD and Eigenvalues are related as worked out in the weekly module. Given a \(2 \times 3\) matrix \(A\) \[A = \begin{bmatrix} 1 & 2 & 3 \\ -1 & 0 & 4 \end{bmatrix} \]
Write code in R to compute \(X = AA^{T}\) and \(Y = A^{T}A\). Then, compute the eigenvalues and eigenvectors of X and Y using the built-in commans in R.
Then, compute the left-singular, singular values, and right-singular vectors of \(A\) using the svd command. Examine the two sets of singular vectors and show that they are indeed eigenvectors of \(X\) and \(Y\). In 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.
1.1 X and Y
Given the \(2\times3\) matrix \(A\), we have \(X\) as a \(2\times2\) matrix and \(Y\) as a \(3\times3\) matrix.
## [,1] [,2]
## [1,] 14 11
## [2,] 11 17## eigen() decomposition
## $values
## [1] 26.601802 4.398198
##
## $vectors
## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043X_evalues <- eigen(X)$values #store the X eigenvalues in X_evalues
X_evectors <- eigen(X)$vectors #store the X eigenvectors in X_evectors
Y <- t(A) %*% A
Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25## eigen() decomposition
## $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.1849001Y_evalues <- eigen(Y)$values #store the Y eigenvalues in Y_evalues
Y_evectors <- eigen(Y)$vectors #store the Y eigenvectors in Y_evectorsBy the above calculations, we can see that:
\(X\) is a \(2\times2\) matrix \(\begin{bmatrix} 14 & 11\\ 11 & 17\end{bmatrix}\) with 2 eigenvalues \(26.602\) and \(4.398\) and their corresponding eigenvectors \(\begin{bmatrix} 0.6576043\\0.7533635\end{bmatrix}\) and \(\begin{bmatrix} -0.7533635\\0.6576043\end{bmatrix}\).
\(Y\) is a \(3\times3\) matrix \(\begin{bmatrix} 2 & 2 & -1 \\ 2 & 4 & 6\\ -1 & 6 & 25\end{bmatrix}\) with 3 eigenvalues \(26.602\), \(4.398\), and \(0\), and their corresponding eigenvectors \(\begin{bmatrix} -0.01856629\\0.25499937\\0.96676296\end{bmatrix}\), \(\begin{bmatrix} -0.6727903\\-0.7184510\\0.1765824\end{bmatrix}\), and \(\begin{bmatrix} 0.7396003\\-0.6471502\\0.1849001\end{bmatrix}\).
1.2 SVD of A
Let’s compute the left-singular, singular values, and right-singular vectors of \(A\) using the svd command.
## $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.1849001As we can see, the matrix \(D\) is a \(3\times3\) matrix. To get the correct dimensions for multiplication, we need to add a blank column to D.
## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0To check if our matrices are correct, compare the multiplication result with our matrix \(A\) by using the equation \(A=U\times \Sigma \times V^{T}\).
## [1] TRUETherefore, we have the left-singular, singular values, and right-singular vectors of \(A\) be: \[U = \begin{bmatrix} -0.6576043 & -0.7533635 \\ -0.7533635 & 0.6576043 \end{bmatrix} \] \[D = \begin{bmatrix} 5.157693 & 0 & 0\\ 0 & 2.097188 & 0\end{bmatrix} \] \[V = \begin{bmatrix} 0.01856629 & -0.6727903 & -0.7396003 \\ -0.25499937 & -0.7184510 & 0.6471502\\ -0.96676296 & 0.1765824 & -0.1849001\end{bmatrix} \]
1.3 Examine X & Y with A
I will examine the two sets of singular vectors and show that they are indeed eigenvectors of \(X\) and \(Y\). I will also show that the two non-zero eigenvalues of both \(X\) and \(Y\) are the same and are squares of the non-zero singular values of \(A\).
1.3.1 U vs X
First, I will compare the left-singular vectors of \(A\) in \(U\) with the eigenvectors of \(X\).
By observation, the first eigenvector of \(X\) has opposite sign of the first singular vector in \(U\), thus I will multiply the first eigenvector of \(X\) by \(-1\) and compare the resulting matrix X_evectors with matrix \(U\).
Note that, multiplying the first eigenvector of \(X\) by a scalar will not affect the comparison with the singular vectors in \(U\) as the modified eigenvector will still be in its original span.
## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043# Compare the eigenvectors of X with the left-singular vectors of A in U
all(round(X_evectors,4) == round(u,4))## [1] TRUEBased on the calculation above, the left-singular vectors of \(A\) in matrix \(U\) are indeed the eigenvectors of \(X\).
1.3.2 V vs Y
Next, I will compare the right-singular vectors of \(A\) in \(V\) with the eigenvectors of \(Y\).
By observation, the first and third eigenvectors of \(Y\) are in the opposite sign of the first and third singular vectors in \(V\), thus I will multiply the first and third eigenvectors of \(Y\) by \(-1\) and compare the resulting matrix Y_evectors with matrix \(V\).
Note that, multiplying the eigenvectors by a scalar will not affect the comparison with the singular vectors in \(V\) as the modified eigenvectors will still be in their original span.
## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001# Multiply the first and third eigenvectors of Y by -1
Y_evectors[,1] <- -1*Y_evectors[,1]
Y_evectors[,3] <- -1*Y_evectors[,3]
Y_evectors## [,1] [,2] [,3]
## [1,] 0.01856629 -0.6727903 -0.7396003
## [2,] -0.25499937 -0.7184510 0.6471502
## [3,] -0.96676296 0.1765824 -0.1849001# Compare the eigenvectors of Y with the right-singular vectors of A in V
all(round(Y_evectors,4) == round(v,4))## [1] TRUEBased on the calculation above, the right-singular vectors of \(A\) in matrix \(V\) are indeed the eigenvectors of \(Y\).
1.3.3 Singular Values vs Eigenvalues
Last, I will compare the non-zero eigenvalues of \(X\) and that of \(Y\), show that they are the same and are the squares of the non-zero singular values of \(A\).
As the third eigenvalue of \(Y\) is nearly zero, I will take it as zero and omit it from the below calculations. Also, by comparing the singular values and the square root of the two non-zero eigenvalues, I will get rid of the third column (all zeros) of matrix \(D\).
## [1] "Eigenvalues of X"## [1] 26.601802 4.398198## [1] "Non-zero eigenvalues of Y"## [1] 26.601802 4.398198# Compare the non-zero eigenvalues of X and Y to see if they are the same
print('Same eigenvalues?')## [1] "Same eigenvalues?"## [1] TRUE## [1] "Square-root of the two non-zero eigenvalues:"## [,1] [,2]
## [1,] 5.157693 0.000000
## [2,] 0.000000 2.097188## [1] "Singular values in matrix D:"## [,1] [,2] [,3]
## [1,] 5.157693 0.000000 0
## [2,] 0.000000 2.097188 0# Compare the square root of the two non-zero eigenvalues with the two singular values
print('Compare the non-zero eigenvalues with singular values, are they equal?')## [1] "Compare the non-zero eigenvalues with singular values, are they equal?"## [1] TRUEThus, the non-zero eigenvalues of \(X\) and \(Y\) are the same and they are equal to the squares of the singular values of A.
2 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. 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 \times 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.
Please submit PS1 and PS2 in an R-markdown document with your first initial and last name.
2.1 General Ideas
According to the Week 4 Lecture, a co-factor \(C_{ij}\) is defined as the determinant of the sub-matrix \(M_{ij}\) together with the appropriate sign \((-1)^{i+j}\), or \((-1)^{i+j}det(M_{ij})\), which simplifies the formula for the determinant of the matrix \(A\) to: \[det(A) = a_{i1}C_{i1} + a_{i2}C_{i2} + a_{i3}C_{i3} + \cdots + a_{in}C_{in}\] where \(C_{ij}\) are the co-factors, for \(n\) be the dimension of the square matrix \(A\) and \(i\) be any row of the square matrix \(A\) where \(i \leq n\).
Also, the inverse of a full-rank square matrix \(A\) can be found by the formula \[A^{-1} = C^{T}/det(A) = \begin{bmatrix} C_{11} & \cdots & C_{n1}\\ \vdots & \ddots & \vdots \\ C_{1n} & \cdots & C_{nn}\end{bmatrix} / \begin{bmatrix} det(A) & \cdots & det(A)\\ \vdots & \ddots & \vdots \\ det(A) & \cdots & det(A)\end{bmatrix}\]
We can compute the determinant of the matrix \(A\) and all its co-factors, take the transpose of the co-factor matrix and divide it by the determinant of the matrix \(A\) to get the inverse.
2.2 Create the function
Create the function myinverse to find the inverse of full-rank square matrix \(A\).
Steps:
Check if the input \(A\) is a matrix
Check if the input matrix \(A\) is square and full-rank
Find all submatrices \(M_{ij}\)
Calculate all co-factors \(C_{ij}\) by using \(C_{ij} = (-1)^{i+j}det(M_{ij})\)
Get the inverse of \(A\) by using \(A^{-1} = C^{T}/det(A)\)
myinverse <- function(A){
#check if A is (1) a matrix, (2) square, and (3) full-rank
if (class(A) != 'matrix')
return('Please input a full-rank square matrix.')
if ((nrow(A)!=ncol(A)) | (nrow(A)!=rankMatrix(A)[1])) #can't use det(A), error if A is not square
return('Please input a full-rank square matrix.')
#initialize the dimension of the full-rank square matrix A
n <- nrow(A)
#create an identity matrix with dimension n to store the cofactors below
cofactor <- diag(n)
#generate all co-factors
for (r in 1:n){
for (c in 1:n){
#find all submatrices by removing one row and one column correspondingly
#in order to prevent the error of removing 1 row and 1 col from a 2x2 matrix
# resulting with a single number, we have to re-define the submatrix as 'matrix'
#please note that, 'byrow' must be FALSE to have the correct submatrix
submatrix <- matrix(A[-r,-c], nrow=n-1, ncol=n-1, byrow=FALSE)
#find all co-factors by multiplying appropriate sign to the det(submatrices)
cofactor[r,c] <- (-1)^(r+c) * det(submatrix)
}
}
#return the inverse of A
return(t(cofactor)/det(A))
}2.3 Test the function
Test the myinverse() function with matrices with different dimensions and verify the results with the R command solve() and matlib function inv().
2.3.1 A Number
Test with a number \(10\).
## [1] "My input:"## [1] 10## [1] "myinverse(A) result:"## [1] "Please input a full-rank square matrix."2.3.2 1x1 matrix
Test with a \(1\times1\) matrix \([10]\).
Note that, myinverse() and solve() can find the inverse of a \(1\times1\) matrix but inv() cannot.
## [1] "My input:"## [,1]
## [1,] 10## [1] "myinverse(A) result:"## [,1]
## [1,] 0.1## [1] "Verify the function using the R command solve()"## [,1]
## [1,] 0.1## Error in apply(abs(X[, 1:n]) <= sqrt(.Machine$double.eps), 1, all) :
## dim(X) must have a positive length## [1] "Are they equal?"## [1] TRUE## [1] "Check A*myinverse(A):"## [,1]
## [1,] 12.3.3 2x2 matrix
Test with a \(2\times2\) matrix \(A = \begin{bmatrix} 1 & 2 \\ 2 & 1\end{bmatrix}\)
## [1] "My input:"## [,1] [,2]
## [1,] 1 2
## [2,] 2 1## [1] "myinverse(A) result:"## [,1] [,2]
## [1,] -0.3333333 0.6666667
## [2,] 0.6666667 -0.3333333## [1] "Verify the function using the R command solve()"## [,1] [,2]
## [1,] -0.3333 0.6667
## [2,] 0.6667 -0.3333## [1] "Verify the function using the matlib function inv()"## [,1] [,2]
## [1,] -0.3333 0.6667
## [2,] 0.6667 -0.3333## [1] "Are they equal?"## [1] TRUE## [1] "Check A*myinverse(A):"## [,1] [,2]
## [1,] 1 0
## [2,] 0 12.3.4 2x3 matrix
Test with a \(2\times3\) matrix \(A = \begin{bmatrix} 1 & 2 & 3 \\ 0 & 1 & -2\end{bmatrix}\).
## [1] "Input matrix A:"## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 0 1 -2## [1] "myinverse(A) result:"## [1] "Please input a full-rank square matrix."2.3.5 3x3 singular matrix
Test with a singular square matrix \(A = \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 6\\ 0 & 1 & -2\end{bmatrix}\).
## [1] "Input matrix A:"## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 2 4 6
## [3,] 0 1 -2## [1] "myinverse(A) result:"## [1] "Please input a full-rank square matrix."2.3.6 4x4 matrix
Test with a full-rank \(4 \times 4\) matrix \(A = \begin{bmatrix} 1 & 2 & 3 & 4\\ -1 & 0 & 1 & 3\\ 0 & 1 & -2 & 1\\ 5 & 4 & -2 & -3\end{bmatrix}\).
A <- matrix(c(1,2,3,4,-1,0,1,3,0,1,-2,1,5,4,-2,-3), nrow=4, ncol=4, byrow=TRUE)
print('Input matrix A:')## [1] "Input matrix A:"## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] -1 0 1 3
## [3,] 0 1 -2 1
## [4,] 5 4 -2 -3## [1] "myinverse(A) result:"## [,1] [,2] [,3] [,4]
## [1,] -2.7778 6.7778 -2.8889 2.1111
## [2,] 3.0000 -7.0000 3.0000 -2.0000
## [3,] 0.8889 -1.8889 0.4444 -0.5556
## [4,] -1.2222 3.2222 -1.1111 0.8889## [1] "Verify the function using the R command solve()"## [,1] [,2] [,3] [,4]
## [1,] -2.7778 6.7778 -2.8889 2.1111
## [2,] 3.0000 -7.0000 3.0000 -2.0000
## [3,] 0.8889 -1.8889 0.4444 -0.5556
## [4,] -1.2222 3.2222 -1.1111 0.8889## [1] "Verify the function using the matlib function inv()"## [,1] [,2] [,3] [,4]
## [1,] -2.7778 6.7778 -2.8889 2.1111
## [2,] 3.0000 -7.0000 3.0000 -2.0000
## [3,] 0.8889 -1.8889 0.4444 -0.5556
## [4,] -1.2222 3.2222 -1.1111 0.8889## [1] "Are they equal?"## [1] TRUE## [1] "Check A*myinverse(A):"## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 12.3.7 6x6 matrix
Test with a \(6\times6\) matrix \(A = \begin{bmatrix}3 & -9 & -13 & 22 & 18 & -24\\ -25 & 21 & -12 & -7 & 23 & -27\\ 12 & -24 & 27 & -8 & 24 & 16\\ 9 & -4 & -9 & 17 & 26 & -22\\ 21 & -29 & 11 & -30 & 10 & -31\\ -21 & 3 & 26 & -14 & 37 & 19\end{bmatrix}\)
A <- matrix(c(3,-9,-13,22,18,-24,-25,21,-12,-7,23,-27,12,-24,27,-8,24,16,9,-4,-9,17,26,-22,21,-29,11,-30,10,-31,-21,3,26,-14,37,19), nrow=6, ncol=6, byrow=TRUE)
print('Input matrix A:')## [1] "Input matrix A:"## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 3 -9 -13 22 18 -24
## [2,] -25 21 -12 -7 23 -27
## [3,] 12 -24 27 -8 24 16
## [4,] 9 -4 -9 17 26 -22
## [5,] 21 -29 11 -30 10 -31
## [6,] -21 3 26 -14 37 19## [1] "myinverse(A) result:"## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.0656 0.0006 0.0154 0.0663 -0.0036 -0.0241
## [2,] -0.0501 -0.1886 -0.2779 0.1494 0.0742 0.1968
## [3,] 0.0458 -0.5575 -0.7600 0.2480 0.2367 0.5790
## [4,] 0.0463 -0.2481 -0.3284 0.1036 0.0891 0.2478
## [5,] -0.0296 0.1725 0.2381 -0.0468 -0.0738 -0.1674
## [6,] -0.0357 0.2747 0.3953 -0.1222 -0.1302 -0.2888## [1] "Verify the function using the R command solve()"## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.0656 0.0006 0.0154 0.0663 -0.0036 -0.0241
## [2,] -0.0501 -0.1886 -0.2779 0.1494 0.0742 0.1968
## [3,] 0.0458 -0.5575 -0.7600 0.2480 0.2367 0.5790
## [4,] 0.0463 -0.2481 -0.3284 0.1036 0.0891 0.2478
## [5,] -0.0296 0.1725 0.2381 -0.0468 -0.0738 -0.1674
## [6,] -0.0357 0.2747 0.3953 -0.1222 -0.1302 -0.2888## [1] "Verify the function using the matlib function inv()"## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.0656 0.0006 0.0154 0.0663 -0.0036 -0.0241
## [2,] -0.0501 -0.1886 -0.2779 0.1494 0.0742 0.1968
## [3,] 0.0458 -0.5575 -0.7600 0.2480 0.2367 0.5790
## [4,] 0.0463 -0.2481 -0.3284 0.1036 0.0891 0.2478
## [5,] -0.0296 0.1725 0.2381 -0.0468 -0.0738 -0.1674
## [6,] -0.0357 0.2747 0.3953 -0.1222 -0.1302 -0.2888## [1] "Are they equal?"## [1] TRUE## [1] "Check A*myinverse(A):"## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 1 0 0 0 0
## [3,] 0 0 1 0 0 0
## [4,] 0 0 0 1 0 0
## [5,] 0 0 0 0 1 0
## [6,] 0 0 0 0 0 1