\(A^TA \ne AA^T\) in general if and only if we cannot demonstrate a matrix \(A\) for which \(A^TA = AA^T\).
Let us define the matrix \(A\) as:
\(A = \begin{vmatrix} 1 & 3 \\ -8 & 6 \end{vmatrix}\)
A <- matrix(c(1, 3, -8, 6), nrow=2, ncol=2, byrow=TRUE)Then \(A^T = \begin{vmatrix} 1 & -8 \\ 3 & 6 \end{vmatrix}\).
Atranspose <- t(A)Does \(A^TA = AA^T\)?
We first multiply \(A^TA\) and obtain the resulting matrix \(A^TA = \begin{vmatrix} 65 & -45 \\ -45 & 45 \end{vmatrix}\).
Atranspose %*% A## [,1] [,2]
## [1,] 65 -45
## [2,] -45 45We then multiply \(AA^T\) and obtain the resulting matrix \(AA^T = \begin{vmatrix} 10 & 10 \\ 10 & 100 \end{vmatrix}\).
A %*% Atranspose## [,1] [,2]
## [1,] 10 10
## [2,] 10 100Since \(A^TA\) does not equal \(AA^T\) in this case, we can say that in general, \(A^TA \ne AA^T\).
We can also use R to check equality without showing the intermediate matrices, if we like. For Example:
A %*% Atranspose == Atranspose %*% A## [,1] [,2]
## [1,] FALSE FALSE
## [2,] FALSE FALSEFirst, let \(I\) be a \(2 \times 2\) identity matrix \(I = \begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix}\)
I <- diag(2)
I## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Then, we see that \(I^T\) is also equal to \(\begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix}\)
Itranspose <- t(I)
Itranspose## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Therefore, \(I = I^T\), which we can also check computationally with R:
I == Itranspose## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEBecause \(I\) and \(I^T\) are equal, \(I^TI\) and \(II^T\) must also be equal. The same will hold true for any \(2 \times 2\) matrix where \(a = d\) and \(b, c = 0\). For example:
M1 <- matrix(c(67, 0, 0, 67), 2, 2, TRUE)
M1 == t(M1)## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEM1 %*% t(M1) == t(M1) %*% M1## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEM2 <- matrix(c(-43, 0, 0, -43), 2, 2, TRUE)
M2 == t(M2)## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEM2 %*% t(M2) == t(M2) %*% M2## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEmatrix_decomposer_LU <- function(M) {
if (nrow(M) != ncol(M)) {
print ("Please provide a square matrix with at least 2 and at most 5 rows/ columns.")
}
else if (nrow(M) == 2) {
E21 <- matrix(c(1, 0, -(M[2, 1]/M[1, 1]), 1), 2, 2, TRUE)
U <- E21 %*% M
L <- solve(E21)
print(U)
print(L)
L %*% U == M
}
else if (nrow(M)==3) {
E21 <- matrix(c(1, 0, 0,
-(M[2, 1]/M[1, 1]), 1, 0,
0, 0, 1), 3, 3, TRUE)
M1 <- E21 %*% M
E31 <- matrix(c(1, 0, 0,
0, 1, 0,
-(M1[3, 1]/M1[1, 1]), 0, 1), 3, 3, TRUE)
M2 <- E31 %*% M1
E32 <- matrix(c(1, 0, 0,
0, 1, 0,
0, -(M2[3, 2]/M2[2, 2]), 1), 3, 3, TRUE)
U <- E32 %*% M2
L <- solve(E32) %*% solve(E31) %*% solve(E21)
print(U)
print(L)
L %*% U == M
}
else if (nrow(M) == 4) {
E21 <- matrix(c(1, 0, 0, 0,
-(M[2, 1]/M[1, 1]), 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1), 4, 4, TRUE)
M1 <- E21 %*% M
E31 <- matrix(c(1, 0, 0, 0,
0, 1, 0, 0,
-(M1[3, 1]/M1[1, 1]), 0, 1, 0,
0, 0, 0, 1), 4, 4, TRUE)
M2 <- E31 %*% M1
E32 <- matrix(c(1, 0, 0, 0,
0, 1, 0, 0,
0, -(M2[3, 2]/M2[2, 2]), 1, 0,
0, 0, 0, 1), 4, 4, TRUE)
M3 <- E32 %*% M2
E41 <- matrix(c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
-(M3[4, 1]/M3[1, 1]), 0, 0, 1), 4, 4, TRUE)
M4 <- E41 %*% M3
E42 <- matrix(c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, -(M4[4, 2]/M4[2, 2]), 0, 1), 4, 4, TRUE)
M5 <- E42 %*% M4
E43 <- matrix(c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0,-(M5[4, 3]/M5[3, 3]), 1), 4, 4, TRUE)
U <- E43 %*% M5
L <- solve(E43) %*% solve(E42) %*% solve(E41) %*% solve(E32) %*% solve(E31) %*% solve(E21)
print(U)
print(L)
L %*% U == M
}
else if (nrow(M) == 5) {
E21 <- matrix(c(1, 0, 0, 0, 0,
-(M[2, 1]/M[1, 1]), 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M1 <- E21 %*% M
E31 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
-(M1[3, 1]/M1[1, 1]), 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M2 <- E31 %*% M1
E32 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, -(M2[3, 2]/M2[2, 2]), 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M3 <- E32 %*% M2
E41 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
-(M3[4, 1]/M3[1, 1]), 0, 0, 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M4 <- E41 %*% M3
E42 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, -(M4[4, 2]/M4[2, 2]), 0, 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M5 <- E42 %*% M4
E43 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0,-(M5[4, 3]/M5[3, 3]), 1, 0,
0, 0, 0, 0, 1), 5, 5, TRUE)
M6 <- E43 %*% M5
E51 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
-(M6[5, 1]/M6[1, 1]), 0, 0, 0, 1), 5, 5, TRUE)
M7 <- E51 %*% M6
E52 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, -(M7[5, 2]/M7[2, 2]), 0, 0, 1), 5, 5, TRUE)
M8 <- E52 %*% M7
E53 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, -(M8[5, 3]/M8[3, 3]), 0, 1), 5, 5, TRUE)
M9 <- E53 %*% M8
E54 <- matrix(c(1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, -(M9[5, 4]/M9[4, 4]), 1), 5, 5, TRUE)
U <- E54 %*% M9
L <- solve(E54) %*% solve(E53) %*% solve(E52) %*% solve(E51) %*% solve(E43) %*% solve(E42) %*% solve(E41) %*% solve(E32) %*% solve(E31) %*% solve(E21)
print(U)
print(L)
L %*% U == M
}
}
B <- matrix(c(3, 0, 1, 5, 4,
2, 7, -5, 8, 1,
0, -100, 8, 3, 1,
-17, 4, 2, 2, 1,
44, 22, 1, 0, 4), 5, 5, byrow=TRUE)
matrix_decomposer_LU(B)## [,1] [,2] [,3] [,4] [,5]
## [1,] 3 0 1.000000e+00 5.000000 4.000000
## [2,] 0 7 -5.666667e+00 4.666667 -1.666667
## [3,] 0 0 -7.295238e+01 69.666667 -22.809524
## [4,] 0 0 1.776357e-15 38.080287 21.209530
## [5,] 0 0 3.920445e-15 0.000000 -3.914157
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.000000 0.0000000 0.000000 0
## [2,] 0.6666667 1.000000 0.0000000 0.000000 0
## [3,] -9.5238095 -14.285714 1.0000000 0.000000 0
## [4,] -3.8621161 2.706826 -0.1494778 1.000000 0
## [5,] 25.8264928 -2.019881 0.2731111 -2.207014 1## [,1] [,2] [,3] [,4] [,5]
## [1,] FALSE TRUE FALSE FALSE FALSE
## [2,] TRUE TRUE TRUE TRUE TRUE
## [3,] FALSE TRUE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE