\((1)\)
Let
\(\mathbf A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}\)
Let
\(\mathbf A^t = \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}\)
We know
\(\mathbf A^t \mathbf A = \begin{bmatrix} b_{11} a_{11} + b_{12} a_{21} & b_{11} a_{12} + b_{12} a_{22} \\ b_{21} a_{11} + b_{22} a_{21} & b_{21} a_{12} + b_{22} a_{22} \end{bmatrix}\)
We also know
\(\mathbf A \mathbf A^t = \begin{bmatrix} a_{11} b_{11} + a_{12} b_{21} & a_{11} b_{12} + a_{12} b_{22} \\ a_{21} b_{11} + a_{22} b_{21} & a_{21} b_{12} + a_{22} b_{22} \end{bmatrix}\)
Therefore, we know that \(\mathbf A^t \mathbf A\) \(\ne\) \(\mathbf A \mathbf A^t\) in general.
A <- matrix(seq(from=1,to=6), nrow=2, byrow=T)
A_t <- t(A)
A_tA = A_t%*%A
AA_t = A%*%A_t
#We know that A_tA is not equal to AA_t
identical(A_tA, AA_t)
## [1] FALSE
\((2)\)
Any symmetric matrix with \(a_{ij} = a_{ji}\) has condition of \(A^t = A\) so these matrices will be \(A^tA = AA = AA^t\)
A <- matrix(c(1,0,0,0,1,0,0,0,1), nrow=3)
A
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
A_t <- t(A)
A_t
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
A_tA = A_t%*%A
AA_t = A%*%A_t
#We know that A_tA is equal to AA_t for Identity Matrix since it is a symmetric matrix
identical(A_tA, AA_t)
## [1] TRUE
# LU_f finds U and creates a list of elimination matrices
# It returns a list of 2 items (U and the elimination list)
LU_f <- 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
A <- matrix(c(1,6,-5,3,2,1,2,6,0),nrow=3)
lu_List <- LU_f(A)
# U is in lu_List[[1]]
U <- lu_List[[1]]
U
## [,1] [,2] [,3]
## [1,] 1 3 2
## [2,] 0 -16 -6
## [3,] 0 0 4
# Elimination matrices are in lu_List[[2]
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)
L
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 6 1 0
## [3,] -5 -1 1
# Confirm that A = LU
(A == L %*% U)
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE