The product of a matrix with its transpose always results in a matrix that is square and symmetric, which are prerequisites to perform diagonalization. Diagonalization decomposes a complex matrix into an easier-to-use product of scalars and vectors.
When we take the product of \(A\) with its transpose (using the built-in
t()function), we see that \(X\) and \(Y\) are both square and symmetric.
# Given a matrix A
A <- matrix(c(1, -1, 2, 0, 3, 4), nrow = 2, ncol = 3)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 4# X is equal to A * AT
X <- A %*% t(A)
X## [,1] [,2]
## [1,] 14 11
## [2,] 11 17# Y is equal to AT * A
Y <- t(A) %*% A
Y## [,1] [,2] [,3]
## [1,] 2 2 -1
## [2,] 2 4 6
## [3,] -1 6 25We can use the
eigen()function to get the normalized eigenvalues and eigenvectors of \(X\) and \(Y\).Eigenvectors and eigenvalues help us diagonalize the matrices \(X\) and \(Y\), which is necessary to find the singular value decomposition of \(A\).
# Get the eigen decomposition of X
eigen_x <- eigen(X)
# Nonzero eigenvalues for X
x_values <- eigen_x$values
x_values## [1] 26.601802 4.398198# Eigenvectors for X
x_vectors <- eigen_x$vectors
x_vectors## [,1] [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635 0.6576043# Get the eigen decomposition of Y
eigen_y <- eigen(Y)
# Nonzero eigenvalues for Y
y_values <- eigen_y$values[1:2]
y_values## [1] 26.601802 4.398198# Eigenvectors for Y
y_vectors <- eigen_y$vectors
y_vectors## [,1] [,2] [,3]
## [1,] -0.01856629 -0.6727903 0.7396003
## [2,] 0.25499937 -0.7184510 -0.6471502
## [3,] 0.96676296 0.1765824 0.1849001We can use the
svd()function to compute the singular value decomposition of \(A\).
The left-singular values are equal to \(U\), the normalized eigenvalues of \(X\)
The singular values are \(\sqrt{\lambda_i}\), or the square root of the eigenvalues of \(X\) or \(Y\)
The right-singular values are equal to \(V\), the normalized eigenvalues of \(Y\)
# Calculate the singular value decomposition of A
svd_A <- svd(A)
# Set the left-singular matrix to be U
left_singular <- svd_A$u
left_singular## [,1] [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635 0.6576043# Set the singular values to be D
singular_values <- svd_A$d
singular_values## [1] 5.157693 2.097188# Set the right-singular matrix to be V
right_singular <- svd_A$v
right_singular## [,1] [,2]
## [1,] 0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296 0.1765824When we compare the two sets of singular vectors to the two sets of eigenvectors from \(X\) and \(Y\), we see they are extremely similar, only differing by a sign.
Here is a comparison between \(X\) and the left-singular matrix of \(A\):
# Eigenvector of X
round(solve(x_vectors), 3)## [,1] [,2]
## [1,] 0.658 0.753
## [2,] -0.753 0.658# Left singular matrix of A
round(left_singular, 3)## [,1] [,2]
## [1,] -0.658 -0.753
## [2,] -0.753 0.658Here is a comparison between \(Y\) and the right-singular matrix of \(A\):
# Eigenvector of Y
round(y_vectors, 3)## [,1] [,2] [,3]
## [1,] -0.019 -0.673 0.740
## [2,] 0.255 -0.718 -0.647
## [3,] 0.967 0.177 0.185# Right singular matrix of A
round(right_singular, 3)## [,1] [,2]
## [1,] 0.019 -0.673
## [2,] -0.255 -0.718
## [3,] -0.967 0.177We see that the right-singular matrix zeroed out the last column of eigenvectors for \(Y\), in order to match the original dimensions of \(A\).
We see that nonzero eigenvalues of \(X\) are equal to those of \(Y\).
# Test whether the nonzero eigenvalues of X equal those of Y
round(x_values, 3) == round(y_values, 3)## [1] TRUE TRUEround(x_values, 3)## [1] 26.602 4.398round(y_values, 3)## [1] 26.602 4.398We also see that singular values squared are equal to the eigenvalues of \(X\) and \(Y\).
# Test whether the square of the singular values is equal to the eigenvalues of X / Y
round(singular_values ^ 2, 3) == round(y_values, 3)## [1] TRUE TRUEround(singular_values ^ 2, 3)## [1] 26.602 4.398round(y_values, 3)## [1] 26.602 4.398To find the inverse of a matrix \(A\), we can use this formula:
\[ \begin{align} A^{-1} &= \frac{1}{\text{det}(A)} \ \text{adjugate}(A) \\ &= \frac{1}{\text{det}(A)} \ \text{cofactor}(A)^T \end{align} \]
where the adjugate of \(A\) is equal to the transpose of its cofactor matrix. To get a cofactor matrix:
Create a matrix of minors by replacing each element in a matrix with the determinant of the submatrix created from removing the row and column of the element in question.
Flip the signs of the matrix of minors in a checkerboard pattern – first positive, then negagtive, then positive, etc.
myInverse <- function(A) {
# Check whether the matrix is square
if (nrow(A) == ncol(A)) {
# Create a copy of the matrix to preserve the rows and columns
cofactor_A <- A
# Populate the cofactor matrix
for (i in 1:nrow(A)) {
for (j in 1:ncol(A)) {
cofactor_A[i, j] <- (det(A[-i,-j])*(-1)^(i+j))
}
}
# Calculate the adjugate matrix as the transpose of the cofactor
adjugate_A <- t(cofactor_A)
# Calculate the inverse using the determinant of A
inverse_A <- (1 / det(A)) * adjugate_A
# Print the inverse
print(inverse_A)
}
# If the matrix is not square
else {
print("Choose a square matrix!")
}
}Let’s test whether my function is accurate by comparing its output to the built-in
solvefunction.
# Test matrix
A <- matrix(c(-1,2,3,-2,1,4,2,1,5), nrow=3)
myInverse(A)## [,1] [,2] [,3]
## [1,] 0.04347826 0.78260870 -0.1739130
## [2,] -0.30434783 -0.47826087 0.2173913
## [3,] 0.21739130 -0.08695652 0.1304348solve(A)## [,1] [,2] [,3]
## [1,] 0.04347826 0.78260870 -0.1739130
## [2,] -0.30434783 -0.47826087 0.2173913
## [3,] 0.21739130 -0.08695652 0.1304348round(myInverse(A), 3) == round(solve(A), 3)## [,1] [,2] [,3]
## [1,] 0.04347826 0.78260870 -0.1739130
## [2,] -0.30434783 -0.47826087 0.2173913
## [3,] 0.21739130 -0.08695652 0.1304348## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUEThe outputs match, so my function accurately calculates the inverse of a square matrix.