Recall the rule, if A is an n × m matrix and B is an m × p matrix, then AB is an n x p matrix by definition.
Assume A is an m x n matrix.
\[ \begin{equation*} \mathbf{A} = \begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{m1}&a_{m2}&\cdots &a_{mn}\\\end{bmatrix} \end{equation*} \]
Therefore AT is an n x m matrix.
\[ \begin{equation*} \mathbf{A^T} =\begin{bmatrix}a_{11}&a_{21}&\cdots &a_{m1}\\a_{12}&a_{22}&\cdots &a_{m2}\\\vdots &\vdots &\ddots &\vdots \\a_{1n}&a_{2n}&\cdots &a_{nm}\\\end{bmatrix} \end{equation*} \]
A(AT) is an m x m matrix.
\[ \begin{equation*} \mathbf{A_{m x n}} * \mathbf{A^T_{n x m}} = \mathbf{B_{m x m}} \end{equation*} \]
(AT)A is an n x n matrix.
\[ \begin{equation*} \mathbf{A^T_{n x m}} * \mathbf{A_{m x n}} = \mathbf{C_{n x n}} \end{equation*} \]
Since m ≠ n, the two matrices are not of the same dimensions, and thus they cannot be equal.
Illustrate with a 2 x 2 matrix first.
\[ \begin{equation*} \mathbf{A} = \begin{bmatrix}a&b\\c&d\end{bmatrix} \end{equation*} \]
\[ \begin{equation*} \mathbf{A^T} = \begin{bmatrix}a&c\\b&d\end{bmatrix} \end{equation*} \]
\[ \begin{equation*} \mathbf{A(A^T)} = \begin{bmatrix}a&b\\c&d\end{bmatrix} * \begin{bmatrix}a&c\\b&d\end{bmatrix} = \begin{bmatrix}a^2 + b^2&ac + bd\\ac+bd&c^2+d^2\end{bmatrix} \end{equation*} \]
\[ \begin{equation*} \mathbf{A^T(A)} = \begin{bmatrix}a&c\\b&d\end{bmatrix} * \begin{bmatrix}a&b\\c&d\end{bmatrix} = \begin{bmatrix}a^2 + c^2&ab + cd\\ab+cd&b^2+d^2\end{bmatrix} \end{equation*} \]
The matrices are equal if and only if the following expression is true.
\[ \begin{equation*} \begin{bmatrix}a^2 + b^2&ac + bd\\ac+bd&c^2+d^2\end{bmatrix} = \begin{bmatrix}a^2 + c^2&ab + cd\\ab+cd&b^2+d^2\end{bmatrix} ? \end{equation*} \]
The system of necessary conditions is rewritten here:
\[ \begin{equation*} a^2 + b^2 = a^2+c^2\\ c^2 + d^2 = b^2+d^2\\ ac + bd = ab+cd\\ \end{equation*} \]
First, we note that the following must be true.
\[ \begin{equation*} b^2 = c^2 \end{equation*} \]
\[ \begin{equation*} b = c \end{equation*} \]
In this case, \[ \begin{equation*} \mathbf{A} = \mathbf{A^T} \end{equation*} \]
\[ \begin{equation*} If (c = b), \mathbf{A(A^T)} = \mathbf{A^T(A)} = \begin{bmatrix}a&b\\b&d\end{bmatrix} * \begin{bmatrix}a&b\\b&d\end{bmatrix} = \begin{bmatrix}a^2 + b^2&ab + bd\\ab+bd&b^2+d^2\end{bmatrix} \end{equation*} \]
\[ \begin{equation*} c = -b \end{equation*} \]
\[ \begin{equation*} a^2 + b^2 = a^2+(-b)^2\\ (-b)^2 + d^2 = b^2+d^2\\ a(-b) + bd = ab+(-b)d\\ \end{equation*} \]
The first two equations hold and the third can be simplified as follows.
\[ \begin{equation*} -ab + bd = ab-bd\\ 2ab = 2bd\\ ab = bd\\ a = d\\ \end{equation*} \]
Therefore, if c = -b, a must equal d for A(AT) = AT(A).
\[ \begin{equation*} If (c = -b ∧ a = d), \mathbf{A(A^T)} = \mathbf{A^T(A)} = \\\begin{bmatrix}a&b\\-b&a\end{bmatrix} * \begin{bmatrix}a&-b\\b&a\end{bmatrix} =\\ \begin{bmatrix}a&-b\\b&a\end{bmatrix} * \begin{bmatrix}a&b\\-b&a\end{bmatrix} =\\\begin{bmatrix}a^2 + b^2&0\\0&a^2+b^2\end{bmatrix} =\\(a^2 + b^2) \begin{bmatrix}1&0\\0&1\end{bmatrix} \end{equation*} \]
Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
Since I already did a LU Decomposition for a 3x3 matrix in Week 1’s assignment, here I generalize the function to any square matrix, and additionally decompose the diagonal matrix.
## Warning: package 'MASS' was built under R version 3.5.2decompLDU <- function(mat){
m <- sqrt(length(mat))
A <- matrix(mat, nrow=m, byrow=T)
U <- A
L <- matrix(0, nrow = m, ncol = m, byrow = T) # initialize Lower Matrix
for (k in 1:m) {
L[k,k] <- 1
}
D <- matrix(0, nrow = m, ncol = m, byrow = T) # initialize Diagonal Matrix
#LU Composition of m x m matrix
i <- 1
a <- NA
while (i < m) {
j <- i
while (j < m) {
a[j] <- U[j+1,i]/U[i,i]
U[j+1,] <- U[i,] * -a[j] + U[j + 1,]
L[j+1,i] <- a[j]
j <- j + 1
}
i <- i + 1
}
## Selecting the diagonal values for the Diagonal Matrix
for (k in 1:m) {
D[k,k] <- U[k,k]
}
## Converting the Upper Diagonal Matrix pivots to 1
for (k in 1:m) {
U[k,] <- U[k,] / U[k,k]
}
L <- fractions(L)
D <- fractions(D)
U <- fractions(U)
print("Lower Triangular Matrix")
print(L)
print("")
print("Diagonal Matrix")
print(D)
print("")
print("Upper Triangular Matrix")
print(U)
print("")
recomposed <- round(L %*% D %*% U)
print("Original Matrix")
print(A)
print("")
print("Re-composed Matrix for Verification")
print(recomposed)
print("")
print(paste0("The two matrices are identical: ", identical(recomposed, A)))
}Verifying the decomposition using Archetype K from the text.
ArchK <- c(10,18,24,24, -12, 12,-2,-6,0,-18,-30,-21,-23,-30,39,27, 30,36,37,-30,18,24,30,30,-20)
decompLDU(ArchK)## [1] "Lower Triangular Matrix"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 6/5 1 0 0 0
## [3,] -3 -165/118 1 0 0
## [4,] 27/10 93/118 -81/20 1 0
## [5,] 9/5 21/59 -12/5 12/19 1
## [1] ""
## [1] "Diagonal Matrix"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 10 0 0 0 0
## [2,] 0 -118/5 0 0 0
## [3,] 0 0 20/59 0 0
## [4,] 0 0 0 19/10 0
## [5,] 0 0 0 0 -2/19
## [1] ""
## [1] "Upper Triangular Matrix"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 9/5 12/5 12/5 -6/5
## [2,] 0 1 87/59 72/59 9/59
## [3,] 0 0 1 51/10 -6
## [4,] 0 0 0 1 -30/19
## [5,] 0 0 0 0 1
## [1] ""
## [1] "Original Matrix"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 10 18 24 24 -12
## [2,] 12 -2 -6 0 -18
## [3,] -30 -21 -23 -30 39
## [4,] 27 30 36 37 -30
## [5,] 18 24 30 30 -20
## [1] ""
## [1] "Re-composed Matrix for Verification"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 10 18 24 24 -12
## [2,] 12 -2 -6 0 -18
## [3,] -30 -21 -23 -30 39
## [4,] 27 30 36 37 -30
## [5,] 18 24 30 30 -20
## [1] ""
## [1] "The two matrices are identical: TRUE"Note: My function currently breaks down when zeroes appear in the lower portion of the Original Matrix and/or when they appear after row operations. For the a-factorizations that occurs in the function, I would have to build a test to conditionally enable row swaps if zeroes are discovered.