Definition 1.6.1.
Let \(A\) be a square matrix. An \(LU\)-decomposition is a decomposition of the matrix \(A\) of form \(A = LU\) where \(L\) and \(U\) are lower and upper triangular matrices (of the same size), respectively.
Section 1.6 \(LU\)-Facotorization
In this section we introduce matrix factorization called LU-factorization and its use in solving a system of linear equations. We shall also explain with examples, how Sage can be use to find LU-factorization of matrices. This section also includes Sage subroutines for LU-factorization.
Subsection 1.6.1 Dooliitle and Crout’s Factorization
Let \(A = LU\) where \(L=[\ell_{ij}]\) with \(\ell_{ij}=0\)if \(i\lt j\) is lower triangular and \(U=[u_{ij}]\) with \(u_{ij}=0\) if \(i\gt j\) is an upper triangular matrix. Equating the entries of the matrices \(A\) and \(LU\) we get
\begin{equation}
\ell_{i1}u_{1j}+\ell_{i2}u_{2j}+\cdots +\ell_{in}u_{nj} = a_{ij} \tag{1.6.1}
\end{equation}
where \(\ell_{ij}=0\) if \(i\lt j\) and \(u_{ij}=0\) if \(i\gt j\text{.}\)
Note that Eqn.(1.6.1) gives us \(n^2\) equations in \(n(n+1)=n^2+n\) unknowns hence it can be solved by taking arbitrary values for any \(n\) unknowns. One of the simplest choice are
- \(\ell_{ii}=1\) for \(1\leq i\leq n\) this method is called Doolittle method
- \(u_{ii}=1\) for \(1\leq i\leq n\text{,}\) this method is called Crout’s method
The key to \(LU\)-decomposition is being able to transform the given matrix into upper triangular using elementary row operation, that involves adding multiples of rows to rows. Also every row which is replaced using the row operation in obtaining row-echelon form may be modified by using row which is above this row. We do not use row interchanges. Let us assume that we need \(k\) elementary such row operations to transform \(A\) to an upper triangular matrix \(U\text{.}\) That is,
\begin{equation*}
E_kE_{k-1}\ldots E_1 A=U
\end{equation*}
Since elementary matrices are invertible, we can multiply both sides by \((E_kE_{k-1}\ldots E_1)^{-1}\) to get the required \(LU\) -decomposition of \(A\text{.}\) Thus
\begin{equation*}
A=E^{-1}_1 E^{-1}_2\cdots E^{-1}_kU=LU
\end{equation*}
It is very easy to see that \(E^{-1}_1 E^{-1}_2\cdots E^{-1}_k\) is an upper triangular matrix with diagonal entries 1.
We list the following results without proof.
- An invertible matrix admits an \(LU\)-factorization if and only if all its principal minors are non-zero. The factorization is unique if we require that the diagonal of \(L\) (or \(U\)) consist of ones.
- If the matrix is singular, then an \(LU\)factorization may still exist. In fact, a square matrix of rank \(k\) has an \(LU\)-factorization if the first \(k\) principal minors are non-zero.
Example 1.6.2. Non existence of DoLittle \(LU-\)factorization.
Let \(A=\begin{pmatrix} 0 \amp 1 \\1 \amp 0\end{pmatrix}\text{.}\) It is easy to check that \(A\) does not have Doolittle \(LU\)-decomposition.
Example 1.6.3.
Let \(A=\begin{pmatrix} 1\amp -1\amp 5\\2\amp -3\amp 1\\1\amp 3\amp 7 \end{pmatrix}\text{.}\) Find the \(LU-\) decomposition of \(A\) using Doolittle’s method. Let
\begin{equation*}
A=\begin{pmatrix}
1\amp -1\amp 5\\2\amp -3\amp 1\\1\amp 3\amp 7
\end{pmatrix}=
\begin{pmatrix}
1\amp 0\amp 0\\\ell_{21}\amp 1\amp 0\\\ell_{31}\amp \ell_{32}\amp 1
\end{pmatrix}
\begin{pmatrix}
u_{11}\amp u_{12}\amp u_{13}\\0\amp u_{22}\amp u_{23}\\0\amp 0\amp u_{33}.
\end{pmatrix}
\end{equation*}
Which implies
\begin{equation*}
\begin{pmatrix}
1\amp -1\amp 5\\2\amp -3\amp 1\\1\amp 3\amp 7
\end{pmatrix}=
\begin{pmatrix}
u_{11}\amp u_{12}\amp u_{13}\\\ell_{21}u_{11}\amp \ell_{21}u_{12}+u_{22}\amp \ell_{21}u_{13}+u_{23}\\
\ell_{31}u_{11}\amp \ell_{31}u_{12}+\ell_{32}u_{22}\amp \ell_{31}u_{13}+\ell_{32}u_{23}+u_{33}
\end{pmatrix}.
\end{equation*}
Equating the two matrices \(A\) and \(LU\) and solving, we get
\begin{gather*}
u_{11}=1, u_{12}=-1, u_{13}=5\\
\ell_{21}u_{11}=2\Rightarrow\ell_{21}=2\\
\ell_{21}u_{12}+u_{22}=-3\Rightarrow u_{22}=-1 \\
\ell_{21}u_{13}+u_{23}=1\Rightarrow u_{23}=-9 \\
\ell_{31}u_{11}=1\Rightarrow\ell_{31}=1 \\
\ell_{31}u_{12}+\ell_{32}u_{22}=3\Rightarrow \ell_{32}=-4 \\
\ell_{31}u_{13}+\ell_{32}u_{23}+u_{33}=7\Rightarrow u_{33}=-34
\end{gather*}
Hence
\begin{equation*}
A=\begin{pmatrix}
1\amp 0\amp 0\\2\amp 1\amp 0\\1\amp -4\amp 1
\end{pmatrix}
\begin{pmatrix}
1\amp -1\amp 5\\0\amp -1\amp -9\\0\amp 0\amp -34
\end{pmatrix}
\end{equation*}
Example 1.6.4.
Let \(A=\begin{pmatrix}1\amp 1\amp 1\\3\amp 1\amp -3\\1\amp -2\amp -5 \end{pmatrix}\text{.}\) Find the LU-decomposition of \(A\) using Crout’s method. Let
\begin{align*}
A=\amp LU=
\begin{pmatrix}
\ell_{11}\amp 0\amp 0\\\ell_{21}\amp \ell_{22}\amp 0\\\ell_{31}\amp \ell_{32}\amp \ell_{33}
\end{pmatrix}
\begin{pmatrix}
1\amp u_{12}\amp u_{13}\\0\amp 1\amp u_{23}\\0\amp 0\amp 1
\end{pmatrix} \\
=\amp \begin{pmatrix}
\ell_{11}\amp \ell_{11}u_{12}\amp \ell_{11}u_{13}\\
\ell_{21}\amp \ell_{21}u_{12}+\ell_{22}\amp \ell_{21}u_{13}+\ell_{22}u_{23}\\
\ell_{31}\amp \ell_{31}u_{12}+\ell_{32}\amp \ell_{31}u_{13}+\ell_{32}u_{23}+\ell_{33}
\end{pmatrix}
\end{align*}
Equating the two matrices \(A\) and \(LU\) and solving, we get
\begin{equation*}
\begin{array}{lll}
\ell_{11}=1,\amp u_{12}=1,\amp u_{13}=1\\
\ell_{21}=3\amp \ell_{22}=-2\amp u_{23}=3\\
\ell_{31}=1\amp \ell_{32}=-3\amp \ell_{33}=3
\end{array}
\end{equation*}
Hence
\begin{equation*}
A=\begin{pmatrix}1\amp 1\amp 1\\3\amp 1\amp -3\\1\amp -2\amp -5 \end{pmatrix}=\begin{pmatrix}
1\amp 0\amp 0\\3\amp -2\amp 0\\1\amp -3\amp 3\end{pmatrix}
\begin{pmatrix}
1\amp 1\amp 1\\0\amp 1\amp 3\\0\amp 0\amp 1
\end{pmatrix}
\end{equation*}
Subsection 1.6.2 Solving system of equations using LU factorization
LU-factorization is very useful in solving system of linear equation. Let \(AX=b\) represents \(n \) equations in \(n \) variables and that \(A=LU \) is a LU factorization of \(A\text{.}\) Then \(AX=LUX=b \) can be written as \(LZ=b \) where \(UX=Z\text{.}\) Now \(LZ=b\) can be solved using forward substitution method and \(UX=z\) can be solved using the backward substitution resulting in solution of the system. LU-factorization method allows to solve a whole lot of linear equations having same coefficient matrix.
Example 1.6.5.
Solve the system of equations \(Ax=b\) using using Doolittle method, where \(A=\begin{pmatrix} 1\amp -1\amp 5\\2\amp -3\amp 1\\1\amp 3\amp 7 \end{pmatrix} \) and \(b=\begin{pmatrix} 7\\-6\\5\end{pmatrix}\)
Solution.
From Example Example 1.6.3, we have \(L=\begin{pmatrix}
1\amp 0\amp 0\\2\amp 1\amp 0\\1\amp -4\amp 1
\end{pmatrix}\) and \(U=\begin{pmatrix}
1\amp -1\amp 5\\0\amp -1\amp -9\\0\amp 0\amp -34
\end{pmatrix}\) Let \(X=\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}\text{,}\) \(Z=\begin{pmatrix}z_1\\z_2\\z_3\end{pmatrix}\text{.}\) \(AX=B\) is equivalent to \(LU=B\) which can be written as \(LZ=b\) where \(Z=UX\text{.}\) First we solve
\begin{equation*}
LZ=b\Longrightarrow \begin{pmatrix}
z_1\\2z_1+z_2\\z_1-4z_2+z_3
\end{pmatrix}=\begin{pmatrix} 7\\-6\\5\end{pmatrix}
\end{equation*}
Using forward substitution we get, \(z_1=7, z_2=-20, z_3=-82\text{.}\) Now, we the required solution by solving \(UX=Z\text{.}\)
\begin{equation*}
UX=Z\Longrightarrow \begin{pmatrix}
x_1-x_2+5x_3\\-x_2-9x_3\\-34x_3
\end{pmatrix}=\begin{pmatrix} 7\\-20\\-82\end{pmatrix}
\end{equation*}
Now solving using the back substitution we get, \(x_1=\frac{115}{17}, x_2=\frac{-29}{17}\) and \(x_3=\frac{41}{17}\text{.}\)
Example 1.6.6.
Solve the following system of equations using LU-decomposition
\begin{equation*}
2x+y+3z=-1; 4x+y+7z=5; -6x-2y-12=-2
\end{equation*}
Solution.
The above system is equivalent to \(AX=B\) where
\begin{equation*}
A=\begin{pmatrix}
2 \amp 1 \amp 3\\ 4 \amp 1 \amp 7 \\ -6 \amp -2 \amp -12
\end{pmatrix}, \quad
X=\begin{pmatrix}
x\\y\\z
\end{pmatrix}
\quad \text {and }
B=\begin{pmatrix}
-1\\5\\-2
\end{pmatrix}
\end{equation*}
First of all let us decompose \(A=LU\) using the elementary row operation.
\begin{align*}
A=\amp\begin{pmatrix}
2 \amp 1 \amp 3\\ 4 \amp 1 \amp 7 \\ -6 \amp -2 \amp -12
\end{pmatrix} \xrightarrow{R_2\to R_2-2R_1}
\begin{pmatrix}
2 \amp 1 \amp 3\\ 0 \amp -1 \amp 1 \\ -6 \amp -2 \amp -12
\end{pmatrix}\\
\amp =\begin{pmatrix}
1 \amp 0 \amp 0\\ -2 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1
\end{pmatrix}A=E_1A \\
\amp \xrightarrow{R_3\to R_3+3R_1}
\begin{pmatrix}
2 \amp 1 \amp 3\\ 0 \amp -1 \amp 1 \\ 0 \amp 1 \amp -3
\end{pmatrix} \\
\amp =\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ 3 \amp 0 \amp 1
\end{pmatrix}E_1A=E_2E_1A\\
\amp \xrightarrow{R_3\to R_2+R_2}
\begin{pmatrix}
2 \amp 1 \amp 3\\ 0 \amp -1 \amp 1 \\ 0 \amp 0 \amp -2
\end{pmatrix}\\
= \amp\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ 0 \amp 1 \amp 1
\end{pmatrix}E_2E_1A=E_3E_2E_1A
\end{align*}
Note that
\begin{equation*}
E_1=
\begin{pmatrix}
1 \amp 0 \amp 0\\ -2 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1
\end{pmatrix},
E_2=\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ 3 \amp 0 \amp 1
\end{pmatrix},
E_3=\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ 0 \amp 1 \amp 1
\end{pmatrix}
\end{equation*}
Their inverses are given by
\begin{equation*}
E^{-1}_1=
\begin{pmatrix}
1 \amp 0 \amp 0\\ 2 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1
\end{pmatrix},
E^{-1}_2=\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ -3 \amp 0 \amp 1
\end{pmatrix},
E^{-1}_3=\begin{pmatrix}
1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0 \\ 0 \amp -1 \amp 1
\end{pmatrix}
\end{equation*}
Thus
\begin{equation*}
L=E^{-1}_1E^{-1}_2E^{-1}_3=\begin{pmatrix}
1 \amp 0 \amp 0\\ 2 \amp 1 \amp 0 \\ -3 \amp -1 \amp 1
\end{pmatrix}
\end{equation*}
and
\begin{equation*}
A=LU=\begin{pmatrix}
1 \amp 0 \amp 0\\ 2 \amp 1 \amp 0 \\ -3 \amp -1 \amp 1
\end{pmatrix} \begin{pmatrix}
2 \amp 1 \amp 3\\ 0 \amp -1 \amp 1 \\ 0 \amp 0 \amp -2
\end{pmatrix}
\end{equation*}
Now the equation \(AX=B\) becomes \(LUX=B\text{.}\) Let \(UX=Y=\begin{pmatrix} u\\v\\w\end{pmatrix}\text{.}\) First we solve \(LY=B\text{.}\) That is
\begin{equation*}
\begin{pmatrix}
1 \amp 0 \amp 0\\ 2 \amp 1 \amp 0 \\ -3 \amp -1 \amp 1
\end{pmatrix} \begin{pmatrix} u\\v\\w\end{pmatrix}=\begin{pmatrix} -1\\5\\-2\end{pmatrix}
\end{equation*}
Solving the above system using the forward substitution we get \(u=-1, v=7, w=2\text{.}\) Now substituting the in equation \(UX=Y\) we get
\begin{equation*}
\begin{pmatrix}
2 \amp 1 \amp 3\\ 0 \amp -1 \amp 1 \\ 0 \amp 0 \amp -2
\end{pmatrix} \begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}-1\\7\\2\end{pmatrix}
\end{equation*}
Solving the above system by back substitution we get the required solution \(x=5, y=-8, z=-1\text{.}\)
Now we shall look at under what conditions a matrix admits LU factorization. The following theorem provides a sufficient condition for ensuring that the algorithm of factorization \(A = LU\) does not break down due to division by zero. If \(A\) is a square matrix, then a sub-matrix \(A^{(k)}\) of \(A\) obtained by taking the first \(k\) rows and first \(k\) columns of \(A\) is called a leading principal minor of \(A\).
Theorem 1.6.7.
Let \(A\) be \(n\times n\) matrix. Suppose that the leading principal minor \(A^{(k)}\) for \(1\leq k\lt n-1\) are nonsingular. Then \(A\) admits \(LU\)-factorization. (Note that \(A\) itself need not be non singular.)
Proof.
Proof of this theorem follows by induction on order of \(A\text{.}\)
Not every matrix can have LU-factorization. Let us consider \(A=\begin{pmatrix}0 \amp 1\\ 1 \amp 0\end{pmatrix}\text{.}\) Suppose \(A\) has LU-factorization, say \(L=\begin{pmatrix}1 \amp 0\\ \ell \amp 1\end{pmatrix}\) and \(U=\begin{pmatrix}u_1 \amp u_2\\ 0 \amp u_3\end{pmatrix}\text{.}\) This implies \(u_1=0\) and \(\ell u_1=1\text{.}\) This is a contradiction. However, if we interchange the first and second row of \(A\text{,}\) then it is identity matrix which has LU-factorization with \(L=U=\begin{pmatrix}1 \amp 0\\ 0 \amp
1\end{pmatrix}\text{.}\)
This leads us to a question, under what conditions, \(LU\)-factorization of a matrix exist. We shall show that even if the matrix \(A\) does not satisfy the conditions of Theorem Theorem 1.6.7, by permuting rows and columns it can be transformed into a new matrix \(\tilde{A}\) of the same size that admits an \(LU\)-factorization. Let first show this result for a \(2\times
2\) matrix.
Example 1.6.8.
Let \(A=\begin{pmatrix}a \amp b\\ c \amp d\end{pmatrix}\text{.}\) Show that there exists a permutation matrix \(P\) of order 2 such that \(PA\) admits the \(LU\)-factorization. If \(a\neq 0\text{.}\) Then by Theorem Theorem 1.6.7, \(A\) admits \(LU\)-factorization. Let \(a=0\) and \(c\neq 0\text{.}\) Define \(P:=\begin{pmatrix}0 \amp 1\\ 1 \amp
0\end{pmatrix}\text{.}\) Then \(PA=\begin{pmatrix}c \amp d\\ 0 \amp b\end{pmatrix}\) which admits \(LU\)-factorization by Theorem Theorem 1.6.7. If \(a=c=0\text{.}\) Then the result is trivial and we have
\begin{equation*}
A=\begin{pmatrix}0 \amp b\\ 0 \amp d\end{pmatrix}=\begin{pmatrix}1 \amp 0\\ 0 \amp
1\end{pmatrix}\begin{pmatrix}0 \amp b\\ 0 \amp d\end{pmatrix}.
\end{equation*}
Theorem 1.6.9.
Given any matrix \(A\) there exists a permutation matrix \(P\text{,}\) a unit lower triangular matrix \(L\) and an upper triangular matrix \(U\) such that \(PA=LU\text{.}\)
The proof of above theorem gives an algorithm for constructing the permutation matrix \(P\text{,}\) and the matrices \(L\) and \(U\text{.}\)
Remark 1.6.10.
If \(P\) is a permutation matrix then \(P^{-1}=P^T\) is also a permutation matrix. If \(PA=LU\text{,}\) then we have \(A=P^TLU\text{.}\) Because of this reason such a factorization is also known as \(PLU\) factorization of \(A\).Subsection 1.6.3 \(LU\)-factorization in Sage
Sage has inbulit method ’LU’ to find \(LU\)-factorization. Suppose that \(A\) is an \(m\times n\) matrix, then an LU decomposition in Sage output is a lower-triangular \(m\times m\) matrix \(L\) with every diagonal element equal to 1, and an upper-triangular \(m\times n\) matrix, \(U\) such that the product \(LU\text{,}\) after a permutation of the rows, is then equal to \(A\text{.}\) For the ’plu’ format the permutation is returned as an m x m permutation matrix \(P\) such that \(A = PLU\text{.}\)
Try to explore help document of LU factorization using ’A.LU?’
Subsection 1.6.4 User defined functions for DooLitlte and Crout’s Methods
Crout’s Method function in Sage
