= LU Decomposition =

'''''LU'' decomposition''' is another approach to [[LinearAlgebra/Elimination|Gauss-Jordan elimination]]. It is generalized as '''''A''' = '''LU'''''.

<<TableOfContents>>

----



== Description ==

If a matrix '''''A''''' can be decomposed into [[LinearAlgebra/SpecialMatrices#Lower_Triangular_Matrices|lower]] and [[LinearAlgebra/SpecialMatrices#Upper_Triangular_Matrices|upper triangular matrices]] '''''L''''' and '''''U''''', the system can generally be solved by:
 1. Let ''y = '''U'''x''.
 2. Reformulate the original system as '''''L'''y = b''. This system's first equation simply states the value of ''y,,1,,''. The second equation gives ''y,,2,,'' in terms of ''y,,1,,'', which can be can be solved by substitution. And so on.
 3. Solve for ''y''.
 4. Reformulate the original problem as '''''U'''x = y''. As before, this system's last equation simply states the value of ''x,,n,,''. And so on.
 5. Solve for ''x''.

Furthermore, note that the [[LinearAlgebra/Invertibility|inverse]] matrix '''''A'''^-1^'' can be calculated as '''''U'''^-1^'''L'''^-1^''.

{{{
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 ==

[[LinearAlgebra/Elimination|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 [[LinearAlgebra/SpecialMatrices#Identity_Matrix|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 [[LinearAlgebra/SpecialMatrices#Lower_Triangular_Matrices|lower triangular matrix]]. In this specific example, '''''L''' = '''E''',,2 1,,^-1^'''E''',,3 2,,^-1^''.

See below how '''''E''',,3 2,,'''E''',,2 1,,'' is messier than '''''E''',,2 1,,^-1^'''E''',,3 2,,^-1^'' despite needing to compute [[LinearAlgebra/Invertibility|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