LU Decomposition

LU decomposition is another approach to Gauss-Jordan elimination. It is generalized as A = LU.

Example

LU decomposition is a strategy for simplifying math with A.

julia> using LinearAlgebra

julia> A = [1 2 1; 3 8 1; 0 4 1]
3×3 Matrix{Int64}:
 1  2  1
 3  8  1
 0  4  1

julia> lu(A, NoPivot())
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 3.0  1.0  0.0
 0.0  2.0  1.0
U factor:
3×3 Matrix{Float64}:
 1.0  2.0   1.0
 0.0  2.0  -2.0
 0.0  0.0   5.0

Note that NoPivot() must be specified because Julia wants to do a more efficient factorization.

Computation of Elimination Matrices

Gauss-Jordan elimination begins with identifying the pivot and transforming this:

┌        ┐
│ [1] 2 1│
│  3  8 1│
│  0  4 1│
└        ┘

...into this:

┌         ┐
│ [1] 2  1│
│  0  2 -2│
│  0  4  1│
└         ┘

This transformation was ultimately a linear combination of rows: subtracting three of row 1 from row 2. This can be reformulated with matrices.

julia> A = [1 2 1; 3 8 1; 0 4 1]
3×3 Matrix{Int64}:
 1  2  1
 3  8  1
 0  4  1

julia> E2_1 = [1 0 0; -3 1 0; 0 0 1]
3×3 Matrix{Int64}:
  1  0  0
 -3  1  0
  0  0  1

julia> E2_1 * A
3×3 Matrix{Int64}:
 1  2   1
 0  2  -2
 0  4   1

This elimination matrix is called E_{2 1} because it eliminated cell (2,1), the elimination cell. An elimination matrix is always the identity matrix with the negated multiplier in the elimination cell.

The Gauss-Jordan approach continues with subtracting two of row 2 from row 3. Formulated as matrices instead:

julia> E3_2 = [1 0 0; 0 1 0; 0 -2 1]
3×3 Matrix{Int64}:
 1   0  0
 0   1  0
 0  -2  1

julia> E3_2 * E2_1 * A
3×3 Matrix{Int64}:
 1  2   1
 0  2  -2
 0  0   5

Decomposition as LU

In this specific example, elimination can be written out as (E_{3 2}E_{2 1})A = U.

A preferable form is A = LU, where L takes on the role of all elimination matrices. L may be a lower triangular matrix. In this specific example, L = E_{2 1}^-1E_{3 2}^-1.

See below how E_{3 2}E_{2 1} is messier than E_{2 1}^-1E_{3 2}^-1 despite needing to compute inverses:

julia> E3_2 * E2_1
3×3 Matrix{Int64}:
  1   0  0
 -3   1  0
  6  -2  1

julia> convert(Matrix{Int64}, inv(E2_1) * inv(E3_2))
3×3 Matrix{Int64}:
 1  0  0
 3  1  0
 0  2  1

Furthermore, in this expression, elimination matrices are iteratively appended instead of prepended.

CategoryRicottone