\[ A^{T}A \neq AA^{T} \]
in general.
Since the hypothesis \[ A^{T}A = AA^{T} \]
is a generalization over all matrices, a single counterexample can disprove it. Let’s use the matrix
\[ {A1} = \left[\begin{array} {rrr} 1 & 2 \\ 3 & 4 \end{array}\right] \]
A1 <- matrix(data = c(1, 2, 3, 4), ncol = 2, nrow = 2, byrow = T)
print(t(A1) %*% A1)
## [,1] [,2]
## [1,] 10 14
## [2,] 14 20
print(A1 %*% t(A1))
## [,1] [,2]
## [1,] 5 11
## [2,] 11 25
We see that the two resultant matrices are different. Therefore the hypothesis is false.
This is true when a matrix is orthogonal, i.e. when its transpose is equal to its inverse.
Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
lu_factor <- function(A) {
r <- nrow(A)
L <- matrix(0, nrow = r, ncol = r)
U <- A
for (i in 1:r) {
L[i, i] <- 1
}
for (i in 1:(r-1)) {
for (j in (i+1):r) {
fact <- -U[j, i] / U[i, i]
L[j, i] <- -fact
U[j,] <- U[j,] + fact * U[i,]
}
}
return(list(L, U))
}
Test the function with a few examples:
Test Example 1:
m1 <- matrix(data = c(4, 3, 6, 3), nrow = 2, ncol = 2, byrow=T)
r <-nrow(m1)
res <- lu_factor(m1)
L <- matrix(unlist(res[1]), nrow = r, ncol = r)
U <- matrix(unlist(res[2]), nrow = r, ncol = r)
print(L)
## [,1] [,2]
## [1,] 1.0 0
## [2,] 1.5 1
print(U)
## [,1] [,2]
## [1,] 4 3.0
## [2,] 0 -1.5
print(L %*% U)
## [,1] [,2]
## [1,] 4 3
## [2,] 6 3
Test Example 2:
m1 <- matrix(data = c(3, 1, 4, 2), nrow = 2, ncol = 2, byrow=T)
r <-nrow(m1)
res <- lu_factor(m1)
L <- matrix(unlist(res[1]), nrow = r, ncol = r)
U <- matrix(unlist(res[2]), nrow = r, ncol = r)
print(L)
## [,1] [,2]
## [1,] 1.000000 0
## [2,] 1.333333 1
print(U)
## [,1] [,2]
## [1,] 3 1.0000000
## [2,] 0 0.6666667
print(L %*% U)
## [,1] [,2]
## [1,] 3 1
## [2,] 4 2
Test Example 3:
m1 <- matrix(data = c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, ncol = 3, byrow=T)
r <-nrow(m1)
res <- lu_factor(m1)
L <- matrix(unlist(res[1]), nrow = r, ncol = r)
U <- matrix(unlist(res[2]), nrow = r, ncol = r)
print(L)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 4 1 0
## [3,] 7 2 1
print(U)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 0 -3 -6
## [3,] 0 0 0
print(L %*% U)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9