LU Decomposition

A fundamental step to solving systems of equation is elimination. A mathematical approach grounded in the same methods, but further generalized for automated compuation, is LU Decomposition.

Linear Algebra

Consider the below system of equations:

x + 2y + z = 2
3x + 8y + z = 12
4y + z = 2

The normal step of elimination proceeds like:

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

Critically, note that this elimination step involved subtracting 3 of row 1 from row 2.

This can instead be formulated 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.

Note that the next elimination step will involve subtracting 2 of row 2 from row 3, to the effect of eliminating cell (3,2). Therefore:

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

Simplification with Factorization

It is generally advantageous to refactor (E_{3 2} E_{2 1}) A = U into A = L U, where L takes on the role of all elimination matrices. The notation L is a parallelism to the U notation for upper triangle matrices; L will be a lower triagle matrix.

L can be trivially solved to be E_{2 1}^-1 E_{3 2}^-1.

Furthermore, because elimination matrices are modifications of the identity matrix, their inverses are also trivial to solve.

            -1
  E       E     =     I
   2,3     2,3
                  ┌      ┐
            -1    │ 1 0 0│
  E       E     = │ 0 1 0│
   2,3     2,3    │ 0 0 1│
                  └      ┘
┌       ┐         ┌      ┐
│  1 0 0│   -1    │ 1 0 0│
│ -3 1 0│ E     = │ 0 1 0│
│  0 0 1│  2,3    │ 0 0 1│
└       ┘         └      ┘
                  ┌      ┐
            -1    │ 1 0 0│
          E     = │ 3 1 0│
           2,3    │ 0 0 1│
                  └      ┘

The fundamental reason for factorization is that the singular product of E_{3 2} E_{2 1} is messier than that of E_{2 1}^-1 E_{3 2}^-1

   E        E      =      ?
    3,2      2,1
┌       ┐┌       ┐   ┌        ┐
│ 1  0 0││  1 0 0│   │  1  0 0│
│ 0  1 0││ -3 1 0│ = │ -3  1 0│
│ 0 -2 1││  0 0 1│   │  6 -2 1│
└       ┘└       ┘   └        ┘

   -1      -1
  E       E      =     L
   2,1     3,2    
┌      ┐┌      ┐   ┌      ┐
│ 1 0 0││ 1 0 0│   │ 1 0 0│
│ 3 1 0││ 0 1 0│ = │ 3 1 0│
│ 0 0 1││ 0 2 1│   │ 0 2 1│
└      ┘└      ┘   └      ┘

In particular, cell (3,1) is a 0 in L but is 6 in E_{3 2} E_{2 1}.

Permutation

If a zero pivot is reached, rows must be exchanged. This can be expressed with matrices as P A = L U.

See Permutation Matrices for more information on these matrices and their properties.

CategoryRicottone