Section 4.4 LU Decomposition
Any nonsingular matrix \(A\) can be factored into a lower triangular matrix \(L\) and an upper triangular matrix \(U\) such that \(A = LU\text{.}\)
We obtain the matrix \(U\) by performing Gaussian elimination operations on the matrix \(A\text{.}\) We obtain the matrix \(L\) by performing the inverse of those Gaussian elimination operations on the identity matrix.
Such a factorization exists when \(A\) can be carried to row echelon form using no row interchanges.
Subsection 4.4.1 Assisted Calculation of LU factorization
To find the LU factorization manually with Sage’s assistance, we first find the echelon form of \(A\) (without full reduction to reduced row echelon form). If there is no direct command in Sage to do this (e.g., using only basic row operations), we can perform the row operations manually. This process will output the matrix \(U\text{.}\)
We then reproduce these operations in inverse order applied to the identity matrix. This will output the matrix \(L\text{.}\)
The row operations to bring \(A\) to echelon form are:
-
\(\displaystyle R_3 - R_1 \to R_3\)
-
\(\displaystyle R_3 + 2R_2 \to R_3\)
Applying these to \(A\) gives \(U\text{.}\) The inverse operations applied to the identity matrix are:
-
\(\displaystyle R_3 - 2R_2 \to R_3\)
-
\(\displaystyle R_3 + R_1 \to R_3\)
This gives \(L\text{.}\) We can then verify that \(A = LU\text{.}\)
Let’s start by defining our matrices.
Then we will create \(U\text{.}\)
Finally we will create our \(L\text{.}\)
Then we can check that our decomposition is still equal to \(A\text{.}\)
Subsection 4.4.2 Sage Calculation of LU factorization
To calculate the matrices of the factorization directly, Sage provides the method
LU.
The output of
A.LU(pivot='nonzero') is a triple of matrices \((P, L, U)\text{.}\) If the matrix admits an LU factorization without row interchanges, the first matrix \(P\) is the identity matrix.
We can then check that these values match our calculations above:
If the matrix does not admit an LU factorization, it requires at least one row interchange when being carried to row echelon form. In this case, the matrix admits a PLU factorization, where the first matrix \(P\) records the row interchanges.
Again, the output is a triple of matrices where the first matrix now is not the identity and records the row interchanges, whereas the second matrix is \(L\) and the third matrix is \(U\) as before.
We can verify that \(B = PLU\)
Subsection 4.4.3 Solving a system of equations using LU
If \(A\) can be factored as \(A = LU\text{,}\) the solution to a system \(Ax = b\) can be solved quickly in two stages:
-
First solve \(Ly = b\) for \(y\) by forward substitution.
-
Then solve \(Ux = y\) for \(x\) by back substitution.
We demonstrate these steps in Sage for \(\textbf{Case II}\text{.}\) We define the matrix \(A\) and a vector \(b\) that lies in its column space.
Since the permutation matrix is the identity for this case, we can proceed directly to solve \(Ly = b\text{.}\)
And then we solve \(Ux = y\) for the final solution.
