X <-matrix(c(1,1,3,4,1,1,3,2,5),nrow=3)
X## [,1] [,2] [,3]
## [1,] 1 4 3
## [2,] 1 1 2
## [3,] 3 1 5t(X)## [,1] [,2] [,3]
## [1,] 1 1 3
## [2,] 4 1 1
## [3,] 3 2 5X %*% t(X)## [,1] [,2] [,3]
## [1,] 26 11 22
## [2,] 11 6 14
## [3,] 22 14 35t(X) %*% X## [,1] [,2] [,3]
## [1,] 11 8 20
## [2,] 8 18 19
## [3,] 20 19 38\[If A = \left[ \begin{array}{ccc} 1 & 4 & 3 \\ 1 & 1 & 2\\3 & 1 & 5\end{array}\right], then~ A^T = \left[ \begin{array}{ccc} 1 & 1 & 3 \\ 4 & 1 & 1\\3 & 2 & 5\end{array}\right]\]
\[AA^T = \left[ \begin{array}{ccc} 1 & 4 & 3 \\ 1 & 1 & 2\\3 & 1 & 5\end{array}\right] \left[ \begin{array}{ccc} 1 & 1 & 3 \\ 4 & 1 & 1\\3 & 2 & 5\end{array}\right] = \left[ \begin{array}{ccc} 26 & 11 & 22 \\ 11 & 6 & 14 \\ 22 & 14 & 35 \end{array}\right]\] \[A^TA = \left[ \begin{array}{ccc} 1& 1 & 3 \\ 4 & 1 & 1\\3 & 2 & 5\end{array}\right] \left[ \begin{array}{ccc} 1 & 4 & 3 \\ 1 & 1 & 2\\3 & 1 & 5\end{array}\right] = \left[ \begin{array}{ccc} 11 & 8 & 20 \\ 8 & 18 & 19 \\ 20 & 19 & 38 \end{array}\right]\]
\[\left[ \begin{array}{ccc} 26 & 11 & 22 \\ 11 & 6 & 14 \\ 22 & 14 & 35 \end{array}\right] \ne \left[ \begin{array}{ccc} 11 & 8 & 20 \\ 8 & 18 & 19 \\ 20 & 19 & 38 \end{array}\right], therefore ~A^TA \ne AA^T\]
When a matrix can be inverted, it holds true. For example:
A <- matrix(c(1, 2, 2, 4), nrow = 2)
identical(A %*% t(A), t(A) %*% A)## [1] TRUEMatrix factorization is a very important problem. There are supercomputers built just to do matrix factorization. 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.
A <- matrix(rexp(25), ncol = 5)
LU.decomp <- function(A){
# The creation of upper matrix U starts with the original matrix A
U = A
n = nrow(U)
L = diag(1, ncol(A), nrow(A))
# https://stackoverflow.com/questions/7102117/how-to-create-a-diagonal-matrix-in-r
# In order to loop through two rows at a time, the first loop must complete at n-1
# https://stackoverflow.com/questions/27802246/how-to-loop-through-a-matrix-rows-from-in-and-columns
for (i in 1:(n - 1)){
for (j in (i + 1):n){
if (U[i, i] != 0){
L[j, i] <- U[j, i] / U[i, i]
U[j, ] <- U[j, ] - L[j, i] * U[i, ]
}
}
}
print("Lower Matrix L")
print(L)
print("Upper Matrix U")
print(U)
print("L * U")
print(L%*%U)
}Finally, let’s print the original matrix followed by the function’s output of the Lower, Upper, and L*U
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.2418593 1.4270450 1.1373453 1.51981880 2.66164119
## [2,] 0.0636447 1.3399831 1.6590887 0.04192281 0.20056258
## [3,] 2.0657826 2.5842405 2.3293750 1.23329804 1.42024700
## [4,] 0.2525338 0.4407095 0.4308589 0.43832121 0.30043340
## [5,] 0.9288036 2.8310146 0.4310677 0.05077835 0.04253641LU.decomp(A)## [1] "Lower Matrix L"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00000000 0.0000000 0.00000000 0.000 0
## [2,] 0.05124953 1.0000000 0.00000000 0.000 0
## [3,] 1.66345952 0.1660886 1.00000000 0.000 0
## [4,] 0.20335140 0.1188130 0.05468411 1.000 0
## [5,] 0.74791371 1.3922021 -15.43491800 -102.587 1
## [1] "Upper Matrix U"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.241859 1.427045 1.1373453 1.51981880 2.66164119
## [2,] 0.000000 1.266848 1.6008003 -0.03596718 0.06415472
## [3,] 0.000000 0.000000 0.1715724 -1.28888527 -3.01794075
## [4,] 0.000000 0.000000 0.0000000 0.20401885 -0.08340407
## [5,] 0.000000 0.000000 0.0000000 0.00000000 -57.17529826
## [1] "L * U"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.2418593 1.4270450 1.1373453 1.51981880 2.66164119
## [2,] 0.0636447 1.3399831 1.6590887 0.04192281 0.20056258
## [3,] 2.0657826 2.5842405 2.3293750 1.23329804 1.42024700
## [4,] 0.2525338 0.4407095 0.4308589 0.43832121 0.30043340
## [5,] 0.9288036 2.8310146 0.4310677 0.05077835 0.04253641The output confirms A=LU