Any non-square matrix will always have \(A^TA = AA^T\), because \(A^TA\) and \(AA^T\) will have different dimensions. Recall that an m x n matrix multiplied by an n x p matrix results in an m x p matrix. Now, using a 3x2 matrix as an example below:
A = 3 rows and 2 columns \(A^T\) = 2 rows and 3 columns \(A^TA\) would result in a 2 by 2 matrix \(AA^T\) would result in a 3 by 3 matrix Therefore, \(A^TA \neq AA^T\). We’d see this same dynamic with any non-square matrix.
For square matrices, I ran 10,000 randomly createds square 3 by 3 matrices to see if any had \(A^TA = AA^T\):
options( warn = -1 )
z <- list()
t <- vector()
for (i in 1:10000) {
A <- matrix(c(sample(-9:9, 9, replace=TRUE)),nrow=3)
z[[i]] <- A
if (sum((t(A) %*% A) == (A %*% t(A))) == 9) {
t[i] <- TRUE} else {
t[i] <- FALSE
}
}Below I extract the matrices (if any) which had \(A^TA = AA^T\) out of my 10,000 samples:
s <- list()
w <- which(t)
for (i in 1:length(w)) {s[[i]] <- z[[w[i]]]}
s## [[1]]
## [,1] [,2] [,3]
## [1,] -2 8 -2
## [2,] 8 -1 -2
## [3,] -2 -2 2
##
## [[2]]
## [,1] [,2] [,3]
## [1,] 3 9 1
## [2,] 9 -8 -8
## [3,] 1 -8 9
##
## [[3]]
## [,1] [,2] [,3]
## [1,] -3 7 6
## [2,] 7 -1 2
## [3,] 6 2 3
##
## [[4]]
## [,1] [,2] [,3]
## [1,] -5 -7 -3
## [2,] -7 -5 0
## [3,] -3 0 2Looking at the matrices where \(A^TA = AA^T\), we can see that all of them are symmetric matrices. A symmetric matrix is a matrix where \(A_{ij} = A_{ji}\) for every i and j.
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. Please submit your response in an R Markdown document using our class naming convention, E.g. LFulton_Assignment2_PS2.png
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. If you doing the entire assignment in R, then please submit only one markdown document for both the problems.
# luFunc finds U and creates a list of elimination matrices,
# returning a list of 2 items (U and the elimination list)
luFunc <- function(A)
{ r <- length(A[1,])
E <- list()
c <- 0
for (i in 2:r) {
p <- i - 1
c <- c + 1
if(A[i, p] !=0) {
elim_tmp <- (A[i, p] / A[p, p]) * -1
}
E[[c]] <- diag(r)
E[[c]][i, p] <- elim_tmp
A <- E[[c]] %*% A
if(i + 1 > r) next
else {
if(A[i + 1, p] != 0) {
elim_tmp <- (A[i + 1, p] / A[p, p]) * -1
}
c <- c + 1
E[[c]] <- diag(r)
E[[c]][i + 1, p] <- elim_tmp
}
A <- E[[c]] %*% A
}
return(list(A, E))
}
# Test on a matrix
a1=c(-1,-2,-2)
a2=c(-5,5,5)
a3=c(0,2,-1)
A=matrix(rbind(a1,a2,a3),nrow=3)
lu_List <- luFunc(A)
# U is first item in list
U <- lu_List[[1]]
# Print U
U## [,1] [,2] [,3]
## [1,] -1 -2 -2
## [2,] 0 15 15
## [3,] 5 0 -3# Elimination matrices are 2nd item in list
elims <- lu_List[[2]]
# Function to find L
findL <- function(l) {
for (i in 1:length(l)) {
tmp <- solve(elims[[i]])
if(i == 1) inv_prod <- tmp
else {inv_prod <- inv_prod %*% tmp}
}
return(inv_prod)
}
L <- findL(elims)
# Print L
L## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] 5 1.0 0
## [3,] 5 0.8 1# Confirm that A = LU
(L %*% U == A)## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE# Testing a 2nd matric
A <- matrix(c(1,2,3,1,1,1,2,0,1),nrow=3)
lu_List <- luFunc(A)
U <- lu_List[[1]]
U## [,1] [,2] [,3]
## [1,] 1 1 2
## [2,] 0 -1 -4
## [3,] 0 0 3elims <- lu_List[[2]]
findL <- function(l) {
for (i in 1:length(l)) {
tmp <- solve(elims[[i]])
if(i == 1) inv_prod <- tmp
else {inv_prod <- inv_prod %*% tmp}
}
return(inv_prod)
}
L <- findL(elims)
L## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 2 1 0
## [3,] 3 2 1(L %*% U == A)## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE