We begin by demonstrating that \[A^{T}A \neq AA^T\] by considering a small 2x2 counterexample.
Let \(A\) be the matrix \[A = \left[ \begin{array}{rr} 0 & 0 \\ 1 & 2 \end{array} \right] \text{ and } A^T = \left[ \begin{array}{rr} 0 & 1 \\ 0 & 2 \\ \end{array} \right]\]
We now compute the product of the matrix with its transpose in both orders.
\[ A \cdot A^T = \left[ \begin{array}{rr} 0 & 0 \\ 0 & 5 \end{array} \right] \text{ while the other product gives } A^T \cdot A = \left[ \begin{array}{rr} 1 & 2 \\ 2 & 4 \end{array} \right] \]
This is sufficient to prove that a matrix and its transpose don’t commute by finding one counterexample.
The following numerical experiment shows that a randomly chosen matrix does not commute with its transpose. I conduct an experiment by measuring the set of \(n=1000\) randomly generated 2 by 2 matrices where entries are chosen independently from a uniform \([0,1]\) interval distribution.
I then measure the norm of its commuter matrix defined as: \[ comm(A) = A \cdot A^T - A^T \cdot A\] Clearly A commutes with its transpose iff \[||comm(A)|| = 0\]. As a control, I check if the matrix A is also invertible.
#
# This function expects a square non-empty matrix
# It returns the zero square matrix iff the matrix and its transpose commute
commuter <- function(A)
{
return( A %*% t(A) - t(A) %*% A)
}
is_commuter <- function(A)
{
return( norm( commuter(A), type="m") == 0 )
}
is_invertible <- function(A)
{
return( abs( det(A)) > 1.0e-10 )
}
random_test <- function(n=100)
{
list_matrices = lapply(1:n, function(x) matrix( runif(4), nrow=2, ncol =2 ))
print("Which random 2x2 matrices are commute with their transpose? ")
print( table( unlist( lapply( list_matrices, is_commuter ) ) ) )
print("Which random 2x2 matrices are invertible? ")
print( table( unlist( lapply( list_matrices, is_invertible ) ) ) )
}
random_test(1000)## [1] "Which random 2x2 matrices are commute with their transpose? "
##
## FALSE
## 1000
## [1] "Which random 2x2 matrices are invertible? "
##
## TRUE
## 1000The above findings show that a random matrix does not commute with its transpose but is invertible.
A few standard categories of matrices do commute with their transposes:
According to wikipedia, the spectral theorem solves the problem of matrices which commute with their transpose. It states that a matrix is normal iff it is unitarily similar to a diagonal matrix. However, that theorem is stated in terms of its conjugate transpose \(A^*\). [https://en.wikipedia.org/wiki/Normal_matrix] However we have not studied the spectral theorem yet.
In this exercise, I present an R function that calculates an LU decomposition of a real-valued matrix. The algorithm chosen was called Doolittle’s algorithm and found online in multiple locations. After implementing the algorithm as an R function, I write a test harness to allow multiple test cases to be validated. Afterwards, I show that multiple test cases can be correctly replicated within the test harness. Finally, I conclude by showing a corner case of a 2x2 matrix where no LU decomposition exists.
The best exposition of the algorithm that was implementable in R was: [https://www3.nd.edu/~zxu2/acms40390F11/Alg-LU-Crout.pdf]
This exposition used indices from 1 to \(n\) in contrast to expositions in Wikipedia which used \(0\)-based index notation for the algorithm. I found some typos in exposition which are corrected below.
The R code directly reflects the notation in this section which serves as documentation for the R.
\[ A = \left[ \begin{array}{rrrr} a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{array} \right] = LU = \left[ \begin{array}{rrrr} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{array} \right] \left[ \begin{array}{rrrr} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{array} \right] \] Each entry \(a_{ij}\) is associated with the product of two vectors \(L[i,]\) and \(U[,j]\) which defines an equation in terms of the desired variables which need to be solved. The key idea is that by choosing the correct order in which \(a_{ij}\) are examined, the equations are solvable in terms of previously solved quantities.
The strategy of Doolittle’s algorithm is to start by solving for the first row of the \(U\) matrix by using the associated equations from the first row of \(A\). Then solving for the first column of the \(L\) matrix which is only dependent on the entry \(u_{11}\). The iterative step is to solve the variables in the \(i\)-th row of U and then the \(i\)-th column of L. The final step is to solve for \(u_{nn}\) in terms of previously derived variables.
Thus, the key formulas for the first step are:
\[ \begin{align} u_{1j} = a_{1j} && j = 1,2, \cdots, n &&\text{the first row of U} \\ && && \\ l_{j1} = \frac{a_{j1}}{u_{11}} && j = 2, \cdots ,n && \text{the first column of L} \\ \end{align} \] In the above, note that the \(u_{11}\) must be solved before solving \(l_{j1}\).
The next steps are solving the interior variables along the \(i\)-th row of \(U\) and the \(i\)-th column of \(L\).
\[ \begin{align} \text{for i = 2, 3, } \cdots n-1 && && && \text{ calculate the following terms} \\ u_{ii} && = && a_{ii} - \sum_{t=1}^{i-1} l_{it}u_{ti} && \text{solve the diagonal term} \\ u_{ij} && = && a_{ij} - \sum_{t=1}^{i-1} l_{it}u_{tj} && \text{for } j = i+1, \cdots , n \text{ i-th row of U} \\ l_{ji} && = && \frac{a_{ji} - \sum_{t=1}^{i-1} l_{jt}u_{ti} }{u_{ii}} && \text{for } j = i+1, \cdots , n \text{ i-th column of L} \\ \end{align} \]
Lastly, we have to solve for the right bottom corner variable of the \(U\) matrix:
\[ u_{nn} = a_{nn} - \sum_{t=1}^{n-1} l_{nt}u_{tn} \]
In the next section, I code this logic in R.
The code implementing Doolittle’s algorithm is shown below.
# FUNCTION: DoolittleLUAlgo
#
# PURPOSE: Calculate the Lower and Upper Triangular decomposition of a square
# matrix in R using the Dolittle Algorithm described above.
#
# INPUT: A matrix object with equal rows and columns. Matrix must be non-empty.
# Values must be real numbers.
#
# OUTPUTS: A list of 3 objects:
#
# Element 1: the number of rows in the the square matrix.
# Element 2: The lower triangular matrix L for input A
# Element 3: The upper triangular matrix U for input A
#
# Note that L * U = (input matrix) A
# ----------------------------------------------------------------------------------
DoolittleLUAlgo <- function(A )
{
n = nrow(A) # the input argument must
L = matrix(rep(0,n*n), nrow = n)
U = matrix(rep(0,n*n), nrow = n)
# Set up first row of Upper triangular matrix
# Set up the first column of the
for(j in 1:n)
{
U[1,j] = A[1,j]
L[j,1] = A[j,1]/U[1,1]
L[j,j] = 1
}
i = 2
while( i < n )
{
sum_i = 0
t = 1
while( t < i )
{
sum_i = sum_i + L[i,t]*U[t,i]
t = t + 1
}
# Define the main diagonal entry of U in row i
# ----------------------------------------------
U[i,i] = A[i,i] - sum_i
sum_i = 0
j = i + 1
# Define the U entries of row i
# (to the right side of the main diagonal)
# -----------------------------------------------
while( j <= n)
{
sum_i = 0
for(t in 1:(i-1) )
{
sum_i = sum_i + L[i,t] * U[t,j]
}
U[i,j] = A[i,j] - sum_i
j = j + 1
}
# Define the L entries of column i
# (below the main diagonal of row)
# -------------------------------------------
j = i + 1
while( j <= n)
{
sum_i = 0
for(t in 1:(i-1))
{
sum_i = sum_i + L[j,t] * U[t, i]
}
L[j, i ] = (A[j, i ] - sum_i ) / U[i, i]
j = j + 1
}
i = i + 1
}
# Final corner U[n,n]
# ---------------------------------------------------
sum_i = 0
t = 1
while(t < n)
{
sum_i = sum_i + L[n, t] * U[ t, n]
t = t + 1
}
U[n,n] = A[n,n] - sum_i
# Return the row count, lower matrix L, upper matrix M
# -------------------------------------------------------
retList = list(n, L, U )
}A few helper functions allow us to succinctly specify test cases. Matrix equality is the first key requirement as the object comparison functions in R are not really suitable. In addition, we need a way to specify known LU decompositions for small examples and verify our R function generates LU matrices consistent with expected outputs.
#
# matequal returns TRUE if the two objects are matrices with equal dimension and equal
# numerical value within a small tolerance otherwise it returns FALSE
# ------------------------------------------------------------------------------------
matequal <- function( A, B)
{
return( is.matrix(A) && is.matrix(B) && dim(A) == dim(B) && all( abs(A - B) < 1.0e-10 ) )
}
solveAndValidateTestCase <- function( tag, A, solL, solU )
{
print(paste0("solveAndValidateTestCase for ", tag ) )
ret = DoolittleLUAlgo(A)
L = ret[[2]]
U = ret[[3]]
n = ret[[1]]
print(paste0( " Dimensions: ", n ) )
print(paste0( " Lower triangular matrix is: "))
print( L)
print(paste0( " Does L matches initial guess? ", matequal(L, solL)))
print(paste0(" Upper triangular matrix is: "))
print( U)
print(paste0(" Does U matches initial guess? ", matequal(U, solU)))
print(paste0(" Does the solution work? Does A equal L %*% U ? " , matequal( A, L %*% U) ) )
print( L %*% U)
print("------------------------------------------------------")
}We develop 5 test cases ordered by dimension of the matrix \(A\) from \(dim(A) = 1, 2, 3, 4\). All test cases have known solutions and the examples below when run using my R function all pass the tests.
S1 = matrix(c( 9 ), nrow = 1, byrow = TRUE)
LS1 = matrix(c( 1 ) , nrow = 1, byrow = TRUE)
US1 = matrix(c( 9 ) , nrow = 1, byrow = TRUE)
solveAndValidateTestCase("Test Case 1: 1x1 (Trivial) ", S1, LS1, US1 )## [1] "solveAndValidateTestCase for Test Case 1: 1x1 (Trivial) "
## [1] " Dimensions: 1"
## [1] " Lower triangular matrix is: "
## [,1]
## [1,] 1
## [1] " Does L matches initial guess? TRUE"
## [1] " Upper triangular matrix is: "
## [,1]
## [1,] 9
## [1] " Does U matches initial guess? TRUE"
## [1] " Does the solution work? Does A equal L %*% U ? TRUE"
## [,1]
## [1,] 9
## [1] "------------------------------------------------------"Ex2 = matrix(c( 3, 1, 4, 2), nrow = 2, byrow = TRUE)
ExL2 = matrix(c( 1, 0, 4/3, 1 ), nrow = 2, byrow = TRUE)
ExU2 = matrix(c( 3, 1, 0, 2/3) , nrow = 2, byrow = TRUE)
solveAndValidateTestCase(" Test Case 2 (2x2) ", Ex2, ExL2, ExU2)## [1] "solveAndValidateTestCase for Test Case 2 (2x2) "
## [1] " Dimensions: 2"
## [1] " Lower triangular matrix is: "
## [,1] [,2]
## [1,] 1.000000 0
## [2,] 1.333333 1
## [1] " Does L matches initial guess? TRUE"
## [1] " Upper triangular matrix is: "
## [,1] [,2]
## [1,] 3 1.0000000
## [2,] 0 0.6666667
## [1] " Does U matches initial guess? TRUE"
## [1] " Does the solution work? Does A equal L %*% U ? TRUE"
## [,1] [,2]
## [1,] 3 1
## [2,] 4 2
## [1] "------------------------------------------------------"A3 = matrix(c( 1, 0, 1, 4, 4, 4, 5, 5, 4), nrow = 3, byrow= TRUE )
solL3 = matrix(c( 1, 0, 0, 4, 1, 0, 5, 5/4,1 ), nrow = 3, byrow = TRUE)
solU3 = matrix(c( 1, 0, 1, 0, 4, 0, 0, 0, -1 ), nrow = 3, byrow = TRUE)
solveAndValidateTestCase( "Test Case 3 (3x3) ", A3, solL3, solU3 )## [1] "solveAndValidateTestCase for Test Case 3 (3x3) "
## [1] " Dimensions: 3"
## [1] " Lower triangular matrix is: "
## [,1] [,2] [,3]
## [1,] 1 0.00 0
## [2,] 4 1.00 0
## [3,] 5 1.25 1
## [1] " Does L matches initial guess? TRUE"
## [1] " Upper triangular matrix is: "
## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 0 4 0
## [3,] 0 0 -1
## [1] " Does U matches initial guess? TRUE"
## [1] " Does the solution work? Does A equal L %*% U ? TRUE"
## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 4 4 4
## [3,] 5 5 4
## [1] "------------------------------------------------------"A4 = matrix( c( 1, 2, 4, 3, 8, 14, 2, 6, 13 ), nrow = 3, byrow = TRUE)
solL4 = matrix( c( 1, 0, 0, 3, 1, 0, 2, 1, 1), nrow = 3, byrow = TRUE)
solU4 = matrix( c( 1, 2, 4, 0, 2, 2, 0, 0, 3), nrow = 3, byrow = TRUE)
solveAndValidateTestCase("Test Case 4 (3x3)", A3, solL3, solU3)## [1] "solveAndValidateTestCase for Test Case 4 (3x3)"
## [1] " Dimensions: 3"
## [1] " Lower triangular matrix is: "
## [,1] [,2] [,3]
## [1,] 1 0.00 0
## [2,] 4 1.00 0
## [3,] 5 1.25 1
## [1] " Does L matches initial guess? TRUE"
## [1] " Upper triangular matrix is: "
## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 0 4 0
## [3,] 0 0 -1
## [1] " Does U matches initial guess? TRUE"
## [1] " Does the solution work? Does A equal L %*% U ? TRUE"
## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 4 4 4
## [3,] 5 5 4
## [1] "------------------------------------------------------"A5 = matrix( c( 2, 1, 1, 0, 4, 3,3,1, 8, 7, 9, 5, 6,7,9,8 ), nrow = 4, byrow = TRUE)
solL5 = matrix( c( 1, 0, 0, 0, 2, 1, 0, 0, 4, 3, 1, 0, 3, 4, 1, 1 ), nrow = 4, byrow = TRUE)
solU5 = matrix( c( 2, 1, 1, 0, 0, 1, 1, 1, 0, 0, 2, 2, 0, 0, 0, 2), nrow = 4, byrow = TRUE)
solveAndValidateTestCase("Test Case 5 (4x4)", A5, solL5, solU5)## [1] "solveAndValidateTestCase for Test Case 5 (4x4)"
## [1] " Dimensions: 4"
## [1] " Lower triangular matrix is: "
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 2 1 0 0
## [3,] 4 3 1 0
## [4,] 3 4 1 1
## [1] " Does L matches initial guess? TRUE"
## [1] " Upper triangular matrix is: "
## [,1] [,2] [,3] [,4]
## [1,] 2 1 1 0
## [2,] 0 1 1 1
## [3,] 0 0 2 2
## [4,] 0 0 0 2
## [1] " Does U matches initial guess? TRUE"
## [1] " Does the solution work? Does A equal L %*% U ? TRUE"
## [,1] [,2] [,3] [,4]
## [1,] 2 1 1 0
## [2,] 4 3 3 1
## [3,] 8 7 9 5
## [4,] 6 7 9 8
## [1] "------------------------------------------------------"Necessary and sufficient conditions for the existence of an LU factorization are now known. It was solved in a 1997 paper by Johnson and Okunev. [https://arxiv.org/pdf/math/0506382v1.pdf] The paper gives a well known example of a matrix with no LU factorization. The simplest example is:
\[ A = \left[ \begin{array}{rr} 0 & 1 \\ 1 & 0 \end{array} \right] = LU = \left[ \begin{array}{rr} 1 & 0 \\ l_{21} & 1 \end{array} \right] \left[ \begin{array}{rr} u_{11} & u_{12} \\ 0 & u_{22} \end{array} \right] \] The reason is simple. Obviously \(0 = a_{11} = l_{11}u_{11} = 1 \cdot u_{11}\). This implies \(u_{11} = 0\). Which implies \(U\) is singular since it has a zero column. However, \(A\) is invertible which is a contradiction.
The following test cases has no solution. We set \(U\) equal to the identify matrix. We set \(L\) to some plausible 2x2 matrix.
CCA = matrix(c( 0, 1, 1, 0), nrow = 2, byrow = TRUE)
CCL = matrix(c( 1, 0, 1, 1 ), nrow = 2, byrow = TRUE)
CCU = matrix(c( 1, 0, 0, 1 ) , nrow = 2, byrow = TRUE)
solveAndValidateTestCase(" Corner Case with no solution (2x2) ", CCA, CCL, CCU)## [1] "solveAndValidateTestCase for Corner Case with no solution (2x2) "
## [1] " Dimensions: 2"
## [1] " Lower triangular matrix is: "
## [,1] [,2]
## [1,] 1 0
## [2,] Inf 1
## [1] " Does L matches initial guess? FALSE"
## [1] " Upper triangular matrix is: "
## [,1] [,2]
## [1,] 0 1
## [2,] 0 -Inf
## [1] " Does U matches initial guess? FALSE"
## [1] " Does the solution work? Does A equal L %*% U ? NA"
## [,1] [,2]
## [1,] 0 NaN
## [2,] NaN NaN
## [1] "------------------------------------------------------"Clearly, the above example will cause the code to produce a nonsense answer as the matrix entries are NaN and Inf. This is easily solved by permuting the rows of \(A\) to obtain the identity matrix.
Clearly, LU decomposition may require permutation of rows or columns to allow a solution.