Problem Set 1
1a) Given a 3 × 2 matrix …
… write code in R to compute \(X=A{ A }^{ 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.
A. Because this is a non-symmterical matrix, the transpose will have a different shape. So will our results, yet both matrices will be symmetrical.
Matrix A:
A <- matrix(c(1,2,3, -1,0,4), nrow=2, ncol=3, byrow=T)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4Note that \(X=A{ A }^{ T }\) results in a symmetrical matrix:
X <- A%*%t(A)
X## [,1] [,2]
## [1,] 14 11
## [2,] 11 17Note that \(Y={ A }^{ T }A\) results in a symmetrical matrix:
Y <- t(A)%*%A
Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25Eigenvalues and eigenvectors. Note that these eigenvectors are normalized, i.e., they are scaled to length 1. (Normalization is done by multiplying the vector by the inverse of its length.).
Eigenvalues of Y:
round(eigen(Y)$values,6)## [1] 26.601802 4.398198 0.000000Eigenvectors of Y:
round(eigen(Y)$vectors, 6)## [,1] [,2] [,3]
## [1,] -0.018566 -0.672790 0.73960
## [2,] 0.254999 -0.718451 -0.64715
## [3,] 0.966763 0.176582 0.18490Eigenvalues of X:
round(eigen(X)$values,6)## [1] 26.601802 4.398198Eigenvectors of X:
round(eigen(X)$vectors,6)## [,1] [,2]
## [1,] 0.657604 -0.753364
## [2,] 0.753364 0.657604Singular values. Given m by n matrix A, and eigenvalues \(λ1, λ2 ... λn\) that are the eigenvalues of \({ A }^{ T }A\), then the singular values of A are \({ \sigma }_{ 1 }=\sqrt { { λ }_{ 1 } } ,\quad { \sigma }_{ 2 }=\sqrt { { λ }_{ 2 } } ,\quad ...\quad { \sigma }_{ n }=\sqrt { { λ }_{ n } }\)
Also keep in mind that \(A=UDV^{ T }\), where A is an m by n matrix; V is an m by m orthogonal matrix; V is an n by n orthogonal matrix; and D is an m by n matrix with singular values of A in its main diagonal.
Also note that \(Y={ A }^{ T }A\) will always produce a symmetrical matrix with eigenvalues that are positive or zero, as demonstrated in the previous step.
Singular values using svd():
svd(A)$d## [1] 5.157693 2.097188Left singular vectors of A. Note that this is a 2 by 2 matrix according to our formula above:
svd(A)$u## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043Right singular vectors of A. Note that this is a 3 by 3 matrix according to our formula above:
svd(A)$v## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.17658241b) 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.
A. By inspection of the above and equality tests below, the absolute non-zero eigenvalues of X and Y are the same and are the squares of the singular values of A. Similarly, the absolute eigenvectors of X and Y are the absolute left and right singular vectors of A.
Eigenvalues of X, Y and squares of singular values of A:
# equivalence test
a <- round(svd(A)$d^2,5)
b <- round(eigen(X)$values,5)
c <- round(eigen(Y)$values[-3],5)
a==b## [1] TRUE TRUEa==c## [1] TRUE TRUEb==c## [1] TRUE TRUEProblem Set 2
[W]rite a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors. … [Y]ou 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.
A. After watching a demo on Khan Academy then experimenting with loops and different algorithms on Github and elsewhere, I settled on the following code, which uses a vectorized function to speed computation.
myinverse <- function(matrix) {
# get the input matrix dimensions aand quit if not square
n <- nrow(matrix)
stopifnot(n == ncol(matrix))
# cofactor function computes the determinant for each cell in the matrix of minors
cofactor <- function(matrix, i, j) (-1)^(i+j) * det(matrix[-i, -j, drop = FALSE])
# the adjoint function creates a cofactor matrix from the confacto values
# cofactor is vectorized to iterate over all cells in the input matrix
adjoint <- function(matrix) {
cofactor_matrix <- outer(seq_len(n), seq_len(n),
Vectorize(function(i, j) cofactor(matrix, i, j)))
t(cofactor_matrix)
}
# finally computes the inverse based on the adjoint matrix
1 / det(matrix) * adjoint(matrix)
}Solution using my function:
X <- matrix(c(-1, -2, 2, 2, 1, 1, 3, 4, 5), 3, 3, byrow=T)
myinverse(X)## [,1] [,2] [,3]
## [1,] 0.04347826 0.78260870 -0.1739130
## [2,] -0.30434783 -0.47826087 0.2173913
## [3,] 0.21739130 -0.08695652 0.1304348Identity check:
round(X%*%myinverse(X), 4)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1Solution using pracma package:
pracma::inv(X)## [,1] [,2] [,3]
## [1,] 0.04347826 0.78260870 -0.1739130
## [2,] -0.30434783 -0.47826087 0.2173913
## [3,] 0.21739130 -0.08695652 0.1304348Test a 4x4 matrix from C21 on p. 205 of our text:
Z <- matrix(c(1,1,-1,2,-2,-1,2,-3,1,1,0,2,-1,2,0,2), 4,4, byrow=T)
myinverse((Z))## [,1] [,2] [,3] [,4]
## [1,] 4 2 0 -1
## [2,] 8 4 -1 -1
## [3,] -1 0 1 0
## [4,] -6 -3 1 1Identity check:
round(Z%*%myinverse(Z),4)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1