If this is true in General it should work for any size matrix. To disprove the general case we need to find an instance when it does not work. We will demostrate with a simple 2x2 matrix.
Let \[ A = \begin{bmatrix} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{bmatrix} \]
And \[ A^{T} = \begin{bmatrix} a_{1,1} & a_{2,1} \\ a_{1,2} & a_{2,2} \end{bmatrix} \]
\[ A^{T}A = \begin{bmatrix} [a_{1,1} + a_{2,1}]^{2} & [(a_{1,1} + a_{2,1})(a_{1,2}+a_{2,2})] \\ [(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] & [a_{1,2} + a_{2,2}]^{2} \end{bmatrix} \]
\[ AA^{T} = \begin{bmatrix} [a_{1,1} + a_{1,2}]^{2} & [(a_{1,1} + a_{1,2})(a_{2,1}+a_{2,2})] \\ [(a_{2,1} + a_{2,2}) (a_{1,1} + a_{2,1})] & [a_{2,1} + a_{2,2}]^{2} \end{bmatrix} \]
Note, in general: \[[a_{1,1} + a_{1,2}]^{2} \neq [a_{1,1} + a_{2,1}]^{2} \] and \[[a_{1,2} + a_{2,2}]^{2} \neq [a_{2,1} + a_{2,2}]^{2} \]
In the off diagonal terms: \[ [(a_{1,1} + a_{2,1})(a_{1,2}+a_{2,2})] \neq [(a_{1,1} + a_{1,2})(a_{2,1}+a_{2,2})]\] and \[[(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] \neq [(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] \]
Since \(A^{T}A \neq AA^{T}\) in the 2x2 case, it is not true in general.
To demonsrate:
Let \[ A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} \]
And \[ A^{T} = \begin{bmatrix} 1 & 3 \\ 2 & 4 \end{bmatrix} \]
So, \[ A^{T}A = \begin{bmatrix} (1 + 3)^{2} & (1 + 3)(2+4)] \\ (2 + 4) (1 + 3) & (2 + 4)^{2} \end{bmatrix} = \begin{bmatrix} 16 & 24 \\ 24 & 36 \end{bmatrix} \]
However,
\[ AA^{T} = \begin{bmatrix} (1 + 2)^{2} & (1 + 2)(3+4)] \\ (3 + 4) (1 + 2) & (3 + 4)^{2} \end{bmatrix} = \begin{bmatrix} 9 & 21 \\ 21 & 49 \end{bmatrix} \]
Note from above that in general:
\[[a_{1,1} + a_{1,2}]^{2} \neq [a_{1,1} + a_{2,1}]^{2} \]
\[[a_{1,2} + a_{2,2}]^{2} \neq [a_{2,1} + a_{2,2}]^{2} \]
\[ [(a_{1,1} + a_{2,1})(a_{1,2}+a_{2,2})] \neq [(a_{1,1} + a_{1,2})(a_{2,1}+a_{2,2})]\] \[[(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] \neq [(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] \] Now suppose that \(a_{i,j} = a_{j,i}\)
Then you have :
\[[a_{1,1} + a_{1,2}]^{2} = [a_{1,1} + a_{2,1}]^{2} \]
\[[a_{1,2} + a_{2,2}]^{2} = [a_{2,1} + a_{2,2}]^{2} \]
\[ [(a_{1,1} + a_{2,1})(a_{1,2}+a_{2,2})] = [(a_{1,1} + a_{1,2})(a_{2,1}+a_{2,2})]\]
\[[(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] = [(a_{1,2} + a_{2,2}) (a_{1,1} + a_{2,1})] \]
All because \(a_{1,2} = a_{2,1}\)
These type of Matrices are called Symmetric and meet the condition \(A^{T} = A\) such that \(AA^{T} = A^{T}A = A^2\)
This only work for square, n x n, matrices, because the dimensional of multiplying an (m x n) x (n x m ) matrix is an m x m matrix and conversly, (n x m) x (m x n) is an n x n matrix. So for non-sqaure matricies the final dimensions of the product are not the same for \(AA^{T}\) and \(A^{T}A\).
There is also the case of Anti-symmetric matrices, which are square matrices meeting the condition:
\(a_{i,j} = 0\) for \(i = j\).
and
\(a_{i,j} = -a_{j,i}\) for \(i \neq j\)
When this is the case \(A^{T} = -A\) and \(A^{T}A = AA^{T} = -A^2\).
I’m going to use the Matrix solving code I wrote last week as a start point:
matrixSolve <- function(A, b){
A <- cbind(A,b)
c <- A[2,1]/A[1,1]
A[2,] <- A[2,]-c*A[1,]
d <- A[3,1]/A[1,1]
A[3,] <- A[3,]-d*A[1,]
e <- A[3,2]/A[2,2]
A[3,] <- A[3,]-e*A[2,]
f <- A[2,3]/A[3,3]
A[2,] <- A[2,] - f*A[3,]
g <- A[1,2]/A[2,2]
A[1,] <- A[1,] - g*A[2,]
h <- A[1,3]/A[3,3]
A[1,] <- A[1,] - h*A[3,]
x <- matrix(c(A[1,4]/A[1,1], A[2,4]/A[2,2], A[3,4]/A[3,3]), nrow = 3, ncol = 1)
return(x)
}
What I am going to do is use a for loop to get R to do the factorization automatically, instead of hard-coded.
I used https://stackoverflow.com/questions/27802246/how-to-loop-through-a-matrix-rows-from-in-and-columns to help iterate through the matrix.
matrixLU <- function(A){
# I am building L from an identity matrix
L <- diag(1, nrow = nrow(A), ncol = ncol(A))
# U is built from A using Guassian elimination
U <- A
# You need to stop iteration on the second to last row
for (i in 1:(nrow(A)-1)){
# The second loop needs to be one ahead of the first and stop on the last row
for (j in (i+1):nrow(A)){
if(U[i,i] != 0){
L[j,i] <- U[j,i]/U[i,i] #The multiple for Guassian elimination becomes an element in L.
U[j,] <- U[j,] - L[j,i]*U[i,]
}
}
}
# In Python I could return L and U as a tuple, I haven't figured out how to do similar in R yet
print(L)
print(U)
return(L%*%U) #This should return the original matrix
}
I’ll test using a matrix that was used to demonstrate the technique here: https://www.youtube.com/watch?v=UlWcofkUDDU
A <- matrix(c(1,4,-3,-2,8,5,3,4,7), nrow = 3, ncol = 3, byrow = TRUE)
A
## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7
matrixLU(A)
## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] -2 1.0 0
## [3,] 3 -0.5 1
## [,1] [,2] [,3]
## [1,] 1 4 -3.0
## [2,] 0 16 -1.0
## [3,] 0 0 15.5
## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7