Matrix factorization is a very important problem. There are supercomputers built just to do matrix factorizations. Every second you are on an airplane, matrices are being factorized. Radars that track flights use a technique called Kalman filtering. At the heart of Kalman Filtering is a Matrix Factorization operation. Kalman Filters are solving linear systems of equations when they track your flight using radars. Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
I am using the following pseudocode algorithm for LU decomposition: https://www.researchgate.net/figure/a-LU-decomposition-pseudocode-b-Example-on-LU-decomposition_fig4_220624910
luFunction <- function(A){
# this method only words for square matrices
if (nrow(A) != ncol(A)){
print('LU Decomposition only works with square matrices!')
break
}
n <- nrow(A)
# set up blank matrices for L and U
L <- matrix(, nrow = n, ncol = n)
U <- matrix(, nrow = n, ncol = n)
for (k in 1:n){
U[k,k] <- A[k,k]
L[k,k] <- 1
if (k+1 <= n){
# pivot column calculations
for (i in (k+1):n){
L[i,k] <- A[i,k]/U[k,k]
U[k,i] <- A[k,i]
}
# remaining matrix update
for (i in (k+1):n){
for (j in (k+1):n){
A[i,j] <- A[i,j] - (L[i,k]%*%U[k,j])
}
}
} # if statement
}
L[is.na(L)] <- 0
U[is.na(U)] <- 0
print('The L matrix is:')
print(L)
print('The U matrix is:')
print(U)
#### function end
}Now that we have our function built, letโs check it with a square matrix:
A <- matrix(
c(1, 1, 3, 2, -1, 5, -1, -2, 4),
nrow = 3,
ncol = 3,
byrow = TRUE
)
luFunction(A)## [1] "The L matrix is:"
## [,1] [,2] [,3]
## [1,] 1 0.0000000 0
## [2,] 2 1.0000000 0
## [3,] -1 0.3333333 1
## [1] "The U matrix is:"
## [,1] [,2] [,3]
## [1,] 1 1 3.000000
## [2,] 0 -3 -1.000000
## [3,] 0 0 7.333333