Each element in a product of two matrices, AB, is a dot product of row and column of A and B respectively.
The product of \({A}^{T}A\) is not always equal to \(A{A}^{T}\) because dot product contains different rows and columns.
\[\mathbf{A} = \left[\begin{array} {rrr} {x}_{1} & {x}_{2} \\ {y}_{1} & {y}_{2} \end{array}\right] \] \[\mathbf{A^T} = \left[\begin{array} {rrr} {x}_{1} & {y}_{1} \\ {x}_{2} & {y}_{2} \end{array}\right] \] \[\mathbf{{A}^{T}A}= \left[\begin{array} {rrr} {x}_{1}{x}_{1}+{y}_{1}{y}_{1} & {x}_{1}{x}_{2}+{y}_{1}{y}_{2} \\ {x}_{2}{x}_{1}+{y}_{2}{y}_{1} & {x}_{2}{x}_{2}+{y}_{2}{y}_{2} \end{array}\right] \] \[\mathbf{A{A}^{T}}= \left[\begin{array} {rrr} {x}_{1}{x}_{1}+{x}_{2}{x}_{2} & {x}_{1}{y}_{1}+{x}_{2}{y}_{2} \\ {y}_{1}{x}_{1}+{y}_{2}{x}_{2} & {y}_{1}{y}_{1}+{y}_{2}{y}_{2} \end{array}\right] \]
As you can see the dot product of first element in \({A}^{T}A\) and \(A{A}^{T}\) are not equal. IE: \({x}_{1}{x}_{1}+{y}_{1}{y}_{1} \ne {x}_{1}{x}_{1}+{x}_{2}{x}_{2}\).
Moreover, If A is m x n matrix then \(A^T\) is n x m, given \(m \ne n\). Product of \(A{A}^{T}\) would be m x m and product of \({A}^{T}A\) would be n x n. These dimensions are not equal and thus their products are not either.
# for examle
A <- matrix(c(1,2,5,
4,3,2,
2,1,9),nrow=3)
# matrix A
A## [,1] [,2] [,3]
## [1,] 1 4 2
## [2,] 2 3 1
## [3,] 5 2 9# transpose of matrix A
t(A)## [,1] [,2] [,3]
## [1,] 1 2 5
## [2,] 4 3 2
## [3,] 2 1 9# product of A with transpose of A
A %*% t(A)## [,1] [,2] [,3]
## [1,] 21 16 31
## [2,] 16 14 25
## [3,] 31 25 110# product of transpose of A with A
t(A) %*% A## [,1] [,2] [,3]
## [1,] 30 20 49
## [2,] 20 29 29
## [3,] 49 29 86This occurs when \({A}^{T}\) is inverse of A and their product result is an identity matrix. When a matrix is multiplied by its inverse, its result is always an identity matrix.
Let’s try to prove it.
\[[\begin{array} {rrr} {A}^{T}A = A{A}^{T} = I \\ {A}^{T}A = I \\ A{A}^{T}A = AI \\ A{A}^{T}A{A}^{-1} = AI{A}^{-1}\\ A{A}^{T}I = A{A}^{-1} \\ A{A}^{T} = I \end{array}] \]
# Factorize matrix into L and U matrices
factorize <- function(A){
U <- A
L_list <- list()
U_list <- list(A)
# loop through each column of elimination matrix
for (i in 1:(ncol(U)-1)) {
#loop through each row to reduce element to zero
for(j in (i+1):nrow(U)){
em <- mat.or.vec(nrow(U), ncol(U)) # create matrix of zeros
for(e in 1:nrow(U)){em[e,e] <-1} # add ones to create identity matrix
em[j,i] <- -1 * U[j,i] / U[i,i] # start factoring
#add elimination matrix to create upper matrix
U_list <- append(list(em),U_list)
#add inverse of elimination matrix in reverse order to create lower matrix
L_list <- append(L_list, list(solve(em)))
}
# multiply matrices in order to get upper/lower matrices
U <- Reduce("%*%", U_list)
L <- Reduce("%*%", L_list)
}
return(list("L"=L,"U"=U))
}
# Run factorize function on matrix given in week 2 module
M <- t(matrix(c(1,1,2,
2,1,0,
3,1,1), ncol=3))
M## [,1] [,2] [,3]
## [1,] 1 1 2
## [2,] 2 1 0
## [3,] 3 1 1factorize(M)## $L
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 2 1 0
## [3,] 3 2 1
##
## $U
## [,1] [,2] [,3]
## [1,] 1 1 2
## [2,] 0 -1 -4
## [3,] 0 0 3factorize(M)$L %*% factorize(M)$U #check if LU = M## [,1] [,2] [,3]
## [1,] 1 1 2
## [2,] 2 1 0
## [3,] 3 1 1# test with 4 by 4 matrix example obtained from
# http://math.stackexchange.com/posts/735420/revisions
M <- t(matrix(c(1,-2,-2,-3,
3,-9,0,-9,
-1,2,4,7,
-3,-6,26,2), ncol=4))
M## [,1] [,2] [,3] [,4]
## [1,] 1 -2 -2 -3
## [2,] 3 -9 0 -9
## [3,] -1 2 4 7
## [4,] -3 -6 26 2factorize(M)## $L
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 3 1 0 0
## [3,] -1 0 1 0
## [4,] -3 4 -2 1
##
## $U
## [,1] [,2] [,3] [,4]
## [1,] 1 -2 -2 -3
## [2,] 0 -3 6 0
## [3,] 0 0 2 4
## [4,] 0 0 0 1factorize(M)$L %*% factorize(M)$U ## [,1] [,2] [,3] [,4]
## [1,] 1 -2 -2 -3
## [2,] 3 -9 0 -9
## [3,] -1 2 4 7
## [4,] -3 -6 26 2# test with 4 by 4 matrix example obtained from
# http://mathfaculty.fullerton.edu/mathews/n2003/CholeskyMod.html
M <- t(matrix(c(2,1,1,3,2,
1,2,2,1,1,
1,2,9,1,5,
3,1,1,7,1,
2,1,5,1,8), ncol=5))
M## [,1] [,2] [,3] [,4] [,5]
## [1,] 2 1 1 3 2
## [2,] 1 2 2 1 1
## [3,] 1 2 9 1 5
## [4,] 3 1 1 7 1
## [5,] 2 1 5 1 8factorize(M)## $L
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0 0.0000000 0.0000000 0.0000000 0
## [2,] 0.5 1.0000000 0.0000000 0.0000000 0
## [3,] 0.5 1.0000000 1.0000000 0.0000000 0
## [4,] 1.5 -0.3333333 0.0000000 1.0000000 0
## [5,] 1.0 0.0000000 0.5714286 -0.8571429 1
##
## $U
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.000000e+00 1.000000e+00 1.0 3.000000e+00 2
## [2,] 0.000000e+00 1.500000e+00 1.5 -5.000000e-01 0
## [3,] 0.000000e+00 0.000000e+00 7.0 0.000000e+00 4
## [4,] 0.000000e+00 0.000000e+00 0.0 2.333333e+00 -2
## [5,] 4.440892e-16 2.220446e-16 0.0 1.332268e-15 2factorize(M)$L %*% factorize(M)$U ## [,1] [,2] [,3] [,4] [,5]
## [1,] 2 1 1 3 2
## [2,] 1 2 2 1 1
## [3,] 1 2 9 1 5
## [4,] 3 1 1 7 1
## [5,] 2 1 5 1 8