Show that \(A^TA \neq AA^T\) in general.
Let A be a square 2x2 matrix.
Let \(x_1, x_2, x_3, x_4\) be unique elements of A.
Let \(A = \begin{bmatrix} x_1 & x_2 \\ x_3 & x_4 \end{bmatrix}\).
Then \(A^T = \begin{bmatrix} x_1 & x_3 \\ x_2 & x_4 \end{bmatrix}\).
From \(A\) and \(A^T\), a general example can be found.
\(A^TA = \begin{bmatrix} x_1 & x_3 \\ x_2 & x_4 \end{bmatrix} \times \begin{bmatrix} x_1 & x_2 \\ x_3 & x_4 \end{bmatrix} = \begin{bmatrix} x_1 \cdot x_1 + x_3 \cdot x_3 & x_1 \cdot x_2 + x_3 \cdot x_4 \\ x_2 \cdot x_1 + x_4 \cdot x_3 & x_2 \cdot x_2 + x_4 \cdot x_4 \end{bmatrix}\)
\(AA^T = \begin{bmatrix} x_1 & x_2 \\ x_3 & x_4 \end{bmatrix} \times \begin{bmatrix} x_1 & x_3 \\ x_2 & x_4 \end{bmatrix} = \begin{bmatrix} x_1 \cdot x_1 + x_2 \cdot x_2 & x_1 \cdot x_3 + x_2 \cdot x_4 \\ x_3 \cdot x_1 + x_4 \cdot x_2 & x_3 \cdot x_3 + x_4 \cdot x_4 \end{bmatrix}\)
\(\begin{bmatrix} x_1 \cdot x_1 + x_3 \cdot x_3 & x_1 \cdot x_2 + x_3 \cdot x_4 \\ x_2 \cdot x_1 + x_4 \cdot x_3 & x_2 \cdot x_2 + x_4 \cdot x_4 \end{bmatrix}\) \(\neq\) \(\begin{bmatrix} x_1 \cdot x_1 + x_2 \cdot x_2 & x_1 \cdot x_3 + x_2 \cdot x_4 \\ x_3 \cdot x_1 + x_4 \cdot x_2 & x_3 \cdot x_3 + x_4 \cdot x_4 \end{bmatrix}\)
Now, let A be a square 2x2 matrix.
Let \(A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\).
Then \(A^T = \begin{bmatrix} 1 &3 \\ 2 & 4 \end{bmatrix}\).
It follows that \(A^TA = \begin{bmatrix} 10 & 14 \\ 14 & 20 \end{bmatrix}\)
Similarly, \(AA^T = \begin{bmatrix} 5 & 11 \\ 11 & 25 \end{bmatrix}\)
Thus, in general, \(A^TA \neq AA^T\)
If all elements of A were identical, that would be an instance where \(A^TA = AA^T\).
From the above conclusion, if the elements of A are not unique, then \(A^TA = AA^T\). If \(x_1 = x_2 = x_3 = x_4\), then \(A = A^T\) and it follows that \(A^TA = AA^T\).
In general, if \(x \in \mathbb{N}\), then \(A = x \cdot I\), where \(I\) is the Identity matrix, and \(A^TA = AA^T\)
factorizetwo <- function(x) {
if(!is.matrix(x)){
return("This is not a matrix. Try again.")
break
}
if(nrow(x) != ncol(x)){
return("This is not a square matrix, I cannot help you.")
break
}
if (nrow(x) == 2 & ncol(x) == 2){
l1 <- 1
l2 <- 0
l4 <- 1
l3 <- x[2,1] / x[1,1]
u1 <- x[1,1]
u2 <- x[1,2]
u3 <- 0
u4 <- (x[2,2] - (l3 * u2)) / l4
return(list("L: ", matrix(c(l1,l3,l2,l4), ncol = 2), "U: ", matrix(c(u1,u3,u2,u4), ncol = 2)))
}
}
test <- matrix(c(4,6,3,3), ncol = 2)
test2 <- matrix(c(5,1,2,3), ncol = 2)
factorizetwo(test)## [[1]]
## [1] "L: "
##
## [[2]]
## [,1] [,2]
## [1,] 1.0 0
## [2,] 1.5 1
##
## [[3]]
## [1] "U: "
##
## [[4]]
## [,1] [,2]
## [1,] 4 3.0
## [2,] 0 -1.5factorizetwo(test2)## [[1]]
## [1] "L: "
##
## [[2]]
## [,1] [,2]
## [1,] 1.0 0
## [2,] 0.2 1
##
## [[3]]
## [1] "U: "
##
## [[4]]
## [,1] [,2]
## [1,] 5 2.0
## [2,] 0 2.6The above function is sufficient for 2x2 matrices.
\(A = \begin{bmatrix} 4 & 3 \\ 6 & 3 \end{bmatrix} = LU = \begin{bmatrix} 1 & 0 \\ 1.5 & 1 \end{bmatrix} \times \begin{bmatrix} 4 & 3 \\ 0 & -1.5 \end{bmatrix}\)
\(B = \begin{bmatrix} 5 & 1 \\ 2 & 3 \end{bmatrix} = LU = \begin{bmatrix} 1 & 0 \\ 0.2 & 1 \end{bmatrix} \times \begin{bmatrix} 5 & 2 \\ 0 & -2.6 \end{bmatrix}\)
After considerable effort to understand the processes behind the general form of the LU decomposition, I was able to find a reasonable procedure. Because LU decomposition does not work if the input matrix is singular, I found a function that tests if a matrix is singular, and I included it in my conditions.
I have included four test cases to validate the function.
factorize <- function(x) {
if(is.singular.matrix(x)){
return("Singular Matrix, no LU decomposition")
break
}
if(!is.matrix(x)){
return("This is not a matrix. Try again.")
break
}
if(nrow(x) != ncol(x)){
return("This is not a square matrix, I cannot help you.")
break
}
L <- matrix(0, nrow = nrow(x), ncol = ncol(x))
U <- matrix(0, nrow = nrow(x), ncol = ncol(x))
diag(L) <- rep(1,nrow(x))
for (i in 1:nrow(x)){
for ( j in 1:ncol(x)) {
U[i,j] <- x[i,j]
if ( i > 1 ) {
for ( k in 1:(i - 1)) {
U[i,j] <- U[i,j] - L[i,k] * U[k,j]
}
}
}
if (i < nrow(x)){
for (j in i:nrow(x)){
L[j,i] <- x[j,i]
if (i > 1){
for (k in 1:(i - 1)){
L[j,i] <- L[j,i] - (L[j,k] * U[k,i])
}
}
L[j,i] <- L[j,i] / U[i,i]
}
}
}
return(list("L: ", L, "U: ", U))
}
a <- matrix(c(4,6,3,3), ncol = 2, byrow = TRUE)
b <- matrix(c(5,6,4,3,5,7,8,9,2), ncol = 3, byrow = TRUE)
c <- matrix(c(3,4,2,4,5,4,3,7,6,5,8,3,4,5,6,9), ncol = 4, byrow = TRUE)
#d <- matrix(c(1,1,1,1,1,1,1,1,1), ncol = 3, byrow = TRUE)
factorize(a)## [[1]]
## [1] "L: "
##
## [[2]]
## [,1] [,2]
## [1,] 1.00 0
## [2,] 0.75 1
##
## [[3]]
## [1] "U: "
##
## [[4]]
## [,1] [,2]
## [1,] 4 6.0
## [2,] 0 -1.5factorize(b)## [[1]]
## [1] "L: "
##
## [[2]]
## [,1] [,2] [,3]
## [1,] 1.0 0.0000000 0
## [2,] 0.6 1.0000000 0
## [3,] 1.6 -0.4285714 1
##
## [[3]]
## [1] "U: "
##
## [[4]]
## [,1] [,2] [,3]
## [1,] 5 6.0 4.000000
## [2,] 0 1.4 4.600000
## [3,] 0 0.0 -2.428571factorize(c)## [[1]]
## [1] "L: "
##
## [[2]]
## [,1] [,2] [,3] [,4]
## [1,] 1.000000 0.000 0.0000000 0
## [2,] 1.666667 1.000 0.0000000 0
## [3,] 2.000000 1.125 1.0000000 0
## [4,] 1.333333 0.125 0.7714286 1
##
## [[3]]
## [1] "U: "
##
## [[4]]
## [,1] [,2] [,3] [,4]
## [1,] 3 4.000000e+00 2.0000000 4.0000000
## [2,] 0 -2.666667e+00 -0.3333333 0.3333333
## [3,] 0 -4.440892e-16 4.3750000 -5.3750000
## [4,] 0 3.425831e-16 0.0000000 7.7714286\(A = \begin{bmatrix} 4 & 6 \\ 3 & 3 \end{bmatrix} = LU = \begin{bmatrix} 1 & 0 \\ 0.75 & 1 \end{bmatrix} \times \begin{bmatrix} 4 & 6 \\ 0 & -1.5 \end{bmatrix}\)
\(B = \begin{bmatrix} 5 & 6 & 4 \\ 3 & 5 & 7 \\ 8 & 9 & 2 \end{bmatrix} = LU = \begin{bmatrix} 1 & 0 & 0 \\ 0.6 & -0.4 & 0 \\ 1.6 & 0.17 & 1 \end{bmatrix} \times \begin{bmatrix} 5 & 6 & 4 \\ 0 & 1.4 & 4.6 \\ 0 & -0.84 & -5.189 \end{bmatrix}\)
\(C = \begin{bmatrix} 3 & 4 & 2 & 4 \\ 5 & 4 & 3 & 7 \\ 6 & 5 & 8 & 3 \\ 4 & 5 & 6 & 9 \end{bmatrix} = LU = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 1.67 & 1 & 0 & 0 \\ 2 & 1.125 & 1 & 0 \\ 1.33 & 0.125 & 0.7714286 & 1 \end{bmatrix} \times \begin{bmatrix} 3 & 4 & 2 & 4 \\ 0 & -2.67 & -0.33 & 0.33 \\ 0 & 0 & 4.375 & -5.375 \\ 0 & 0 & 0 & 7.771 \end{bmatrix}\)