library(pracma)
library(tidyverse) Problem Set 1:
Given:
\(A= \begin{bmatrix} 1 & 2 & 3 \\ -1 & 0 & 4 \\ \end{bmatrix}\)
- Compute \(X=AA^T\) and \(Y=A^TA\).
- Compute the eigenvalues and eigenvenctors of \(X\) and \(Y\) using the built-in R commands (e.g.
eigen())
- Compute the left-singular, singular values, and right-singular vectors of \(A\) using the svd command and show that they are indeed eigenvectors of \(X\) and \(Y\).
JB:
Create A:
A1 <- matrix(c(1, 2, 3, -1, 0, 4),
nrow=2, ncol=3, byrow=T)
# sfd output:
# d = a vector containing the singular values of x
# u = a matrix whose columns contain the left singular vectors
# v = a matrix whose columns contain the right singular vectors
mySVD_A<-svd(A1)Use R to calc \(X\) and \(Y\) as stated in the question:
X <- A1 %*% t(A1)
Y <- t(A1) %*% A1Use R to calc the eigenvalues of \(X\) and \(Y\)…
e_val_X <- eigen(X)$values
e_val_Y <- eigen(Y)$values
round(e_val_X,5)## [1] 26.6018 4.3982round(e_val_Y,5)## [1] 26.6018 4.3982 0.0000…and then the eigenvectors of \(X\) and \(Y\)…
e_vec_X <- eigen(X)$vectors
e_vec_Y <- eigen(Y)$vectors
round(e_vec_X,5)## [,1] [,2]
## [1,] 0.65760 -0.75336
## [2,] 0.75336 0.65760round(e_vec_Y,5)## [,1] [,2] [,3]
## [1,] -0.01857 -0.67279 0.73960
## [2,] 0.25500 -0.71845 -0.64715
## [3,] 0.96676 0.17658 0.18490Next, we’ll use the results of using the svd function on \(A\), namely mySVD_A, to show that the two sets of singular vectors produced are eigenvectors of \(X\) and \(Y\):
mySVD_A$u # compare eigenvec## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043e_vec_X # of X## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043mySVD_A$v # compare eigenvec## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824e_vec_Y # of Y## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001rref(mySVD_A$u)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1rref(e_vec_X)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1rref(mySVD_A$v)## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
## [3,] 0 0rref(e_vec_Y)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1Problem Set 2:
Write a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors.
JB:
Below, I’ve outlined the rough steps needed to obtain the inverted matrix:
- create the matrix of minors
- get determinant of each minor-matrix
- repair signs to obtain the co-factor matrix \(C\), transpose it
- get determinant of original matrix
- compute \(\frac{1}{det(A)}adj(A)=\frac{1}{det(A)}C^{T}\) to get the inverse of \(A\), namely \(A^{-1}\) where \(adj(A)\) is the adjugate matrix of A.
Note that myinverse uses code like A[-i,-j] which is the same thing as saying “all rows of \(A\) except the \(i^{th}\) row and all columns of \(A\) except the \(j^{th}\) column” - this allows me to obtain the matrix of minors:
myinverse <- function(an_A){ # Matrix must be square
n <- nrow(an_A)
C <- diag(n) # Matrix to hold inv
for(i in 1:n){ # Get determinent of
for(j in 1:n){ # each sub-matrix
C[i,j] <- det(an_A[-i,-j])*(-1)^(i+j)# change sign as needed
} # an_A[-i,-j] = minor matrix =
} # "not row i & not col j"
C <- C %>% t() # transpose C
B <- 1/det(an_A)*C # calculate A inverse
return(B)
}To test my function, I’ll use the following matrix 3x3 matrix \(A\) and 4x4 matrix, \(X\):
\(A= \begin{bmatrix} 1 & -2 & 2\\ 2 & 1 & 1\\ 3 & 4 & 5\\ \end{bmatrix} \), \(X= \begin{bmatrix} -1 & -2 & 3 & 4\\ 1 & 2 & 7 & 5 \\ -5 & -3 & 1 & 7\\ 1 & 4 & 5 & -6 \\ \end{bmatrix}\)
To show that \(AA^{-1}=I\) and \(XX^{-1}=I\) for both my created function, myinverse, and the base R function solve:
A <- matrix(c(-1, -2, 2, 2, 1, 1, 3, 4, 5),
nrow=3, ncol=3, byrow=T)
X1 <- matrix(c(-1, -2, 3, 4, 1, 2, 7, 5, -5, -3, 1, 7, 1, 4, 5, -6),
nrow=4, ncol=4, byrow=T)
round(A%*%myinverse(A),10) # Test my function ## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1round(A%*%solve(A), 10) # compare to R output## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1round(X1%*%myinverse(X1),10) # Test my function ## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1round(X1%*%solve(X1),10) # compare to R output## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1As shown above, my function myinverse provides the same output as the base R solve function. Note that I have used round with 10 digits of accuracy to clean up the vanishingly small numbers in the resulting matrix and to illustrate how close to zero each of these numbers was - please scroll to the bottom to review the raw output of A %*% myinverse(A):
Notes:
- An important portion of my algorithm comes from the weekly module that provided the logic required to build the co-factor matrix: \((-1)^{i+j}\). This allowed me to create the “checkerboard” plus/minus 1’s within the
forloop to ensure each element of the adjugate matrix had the correct sign.
- This Khan Acadamy video was very helpful in showing the basic steps needed for calculating the inverse of a 3x3 matrix using co-factors.
- Below, after the notes, I’ve shown a couple of non-rounded examples to show the raw results.
A%*%myinverse(A)## [,1] [,2] [,3]
## [1,] 1.000000e+00 3.053113e-16 -1.110223e-16
## [2,] 5.551115e-17 1.000000e+00 0.000000e+00
## [3,] 2.220446e-16 -9.992007e-16 1.000000e+00A%*%solve(A)## [,1] [,2] [,3]
## [1,] 1.000000e+00 -5.551115e-17 1.110223e-16
## [2,] -5.551115e-17 1.000000e+00 -2.775558e-17
## [3,] 0.000000e+00 0.000000e+00 1.000000e+00