library(tinytex)
library(matrixcalc)\[A^T A \neq AA^T\]
Let’s create a 2x 2 square matrix for \(A\): \[ \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \]
Its transpose, \(A^T\): \[ \begin{bmatrix} a & c \\ b & d \\ \end{bmatrix} \]
Verifying that the initial equation is not true: \[ \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \begin{bmatrix} a & c \\ b & d \\ \end{bmatrix} \neq \begin{bmatrix} a & c \\ b & d \\ \end{bmatrix} \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \]
Multiplying the matrices: \[ \begin{bmatrix} a^2 + b^2 & ac + bd \\ ca + db & c^2 + d^2 \\ \end{bmatrix} \neq \begin{bmatrix} a^2 + c^2 & ab + cd \\ ba + dc & b^2 + d^2 \\ \end{bmatrix} \]
Clearly not equivalent for many values.
\[A^T A = AA^T\]
This is true when the matrix is symmetric - meaning a square matrix where the below is true: \[A=A^T\]
The identity matrix is symmetric, but note that being able to be reduced to the identity matrix does not make a matrix symmetric.
Showing with a quick example. The below symmetric matrix is its own transpose, so the stated conditions would be true: \[ \begin{bmatrix} x & y \\ y & x \\ \end{bmatrix} \]
Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer. You don’t have to worry about permuting rows of A and you can assume that A is less than 5x5, if you need to hard-code any variables in your code.
Rough plan: 1. First, verify it’s a matrix. 2. Next, verify it’s a square matrix. 2. Finally, start with \(A\) and go toward figuring out \(LU\) using Gaussian elimination and a shortcut from the YouTube video so that: \[A = LU\]
LU_Function <- function(A) {
#verify matrix
if (is.matrix(A)==FALSE) {
return('Please provide a matrix')
}
#verify matrix is square
#could verify if rows = columns using dim, but we'll be lazy and use package
if (is.square.matrix(A)==FALSE) {
return('Please provide a square matrix (same numbers of rows and columns)')
}
#derive the L and U
m_count <- nrow(A)
#build lower triangular matrix starting with diagonal matrix
#diagonal of 1s with only 0s above
L <- diag(m_count)
#start with U as A for upper triangular matrix
#everything below diag has to be 0 but diag does not have to be 1
U <- A
#start the work
#shortcut method - use opposites of multiples
for (j in (1:m_count)){
for (i in (2:(m_count))){
if(i > j){
working_row <- U[j,]
U_mult <- U[i,j]/working_row[j]
U[i,] <- U[i,] - (U_mult * working_row)
L[i,j] <- U_mult
}
}
}
return(list(A=A,L=L,U=U))
}Test case:
A = matrix(c(2, -3, 1, 7), 2, 2, byrow = TRUE)
LU_Function(A)## $A
## [,1] [,2]
## [1,] 2 -3
## [2,] 1 7
##
## $L
## [,1] [,2]
## [1,] 1.0 0
## [2,] 0.5 1
##
## $U
## [,1] [,2]
## [1,] 2 -3.0
## [2,] 0 8.5Using a function:
lu.decomposition(A)## $L
## [,1] [,2]
## [1,] 1.0 0
## [2,] 0.5 1
##
## $U
## [,1] [,2]
## [1,] 2 -3.0
## [2,] 0 8.5Testing with a bigger matrix:
B = matrix(c(2, -3, 1, 7, 5, 2, 4, 9, 23, 34, 1, 5, 7, 9, 10, 13), 4, 4, byrow = TRUE)
LU_Function(B)## $A
## [,1] [,2] [,3] [,4]
## [1,] 2 -3 1 7
## [2,] 5 2 4 9
## [3,] 23 34 1 5
## [4,] 7 9 10 13
##
## $L
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0.000000 0.0000000 0
## [2,] 2.5 1.000000 0.0000000 0
## [3,] 11.5 7.210526 1.0000000 0
## [4,] 3.5 2.052632 -0.1604938 1
##
## $U
## [,1] [,2] [,3] [,4]
## [1,] 2 -3.0 1.00000 7.000000
## [2,] 0 9.5 1.50000 -8.500000
## [3,] 0 0.0 -21.31579 -14.210526
## [4,] 0 0.0 0.00000 3.666667lu.decomposition(B)## $L
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0.000000 0.0000000 0
## [2,] 2.5 1.000000 0.0000000 0
## [3,] 11.5 7.210526 1.0000000 0
## [4,] 3.5 2.052632 -0.1604938 1
##
## $U
## [,1] [,2] [,3] [,4]
## [1,] 2 -3.0 1.00000 7.000000
## [2,] 0 9.5 1.50000 -8.500000
## [3,] 0 0.0 -21.31579 -14.210526
## [4,] 0 0.0 0.00000 3.666667References:
https://stat.ethz.ch/R-manual/R-devel/library/base/html/dim.html
https://www.youtube.com/watch?v=UlWcofkUDDU
https://www.quantstart.com/articles/LU-Decomposition-in-Python-and-NumPy/