How to Perform an LU Factorization

Begin with the matrix equation. Fundamentally, a system of equations can be written in terms of a matrix equation A x = b , {\displaystyle A\mathbf {x} =\mathbf {b} ,} A{\mathbf {x}}={\mathbf {b}}, where matrix A {\displaystyle A} A acts on a vector x {\displaystyle \mathbf {x} } {\mathbf {x}} to output another vector b . {\displaystyle \mathbf {b} .} {\mathbf {b}}. It is often the case that we wish to know x , {\displaystyle \mathbf {x} ,} {\mathbf {x}}, and this is no exception. In LU factorization, we will see that we can define the relation A = L U , {\displaystyle A=LU,} A=LU, where L {\displaystyle L} L and U {\displaystyle U} U are both triangular matrices. ( 2 4 1 − 1 1 3 3 2 5 ) x = ( 2 1 7 ) {\displaystyle {\begin{pmatrix}2&4&1\\-1&1&3\\3&2&5\end{pmatrix}}\mathbf {x} ={\begin{pmatrix}2\\1\\7\end{pmatrix}}} {\begin{pmatrix}2&4&1\\-1&1&3\\3&2&5\end{pmatrix}}{\mathbf {x}}={\begin{pmatrix}2\\1\\7\end{pmatrix}}

Row-reduce A {\displaystyle A} A to row-echelon form. The row-echelon form will become our matrix U . {\displaystyle U.} U. R 2 → 2 R 2 + R 1 , R 3 → 2 R 3 − 3 R 1 {\displaystyle R_{2}\to 2R_{2}+R_{1},\ R_{3}\to 2R_{3}-3R_{1}} R_{{2}}\to 2R_{{2}}+R_{{1}},\ R_{{3}}\to 2R_{{3}}-3R_{{1}} ( 2 4 1 0 6 7 0 − 8 7 ) {\displaystyle {\begin{pmatrix}2&4&1\\0&6&7\\0&-8&7\end{pmatrix}}} {\begin{pmatrix}2&4&1\\0&6&7\\0&-8&7\end{pmatrix}} R 3 → 3 R 3 + 4 R 2 {\displaystyle R_{3}\to 3R_{3}+4R_{2}} R_{{3}}\to 3R_{{3}}+4R_{{2}} ( 2 4 1 0 6 7 0 0 49 ) = U {\displaystyle {\begin{pmatrix}2&4&1\\0&6&7\\0&0&49\end{pmatrix}}=U} {\begin{pmatrix}2&4&1\\0&6&7\\0&0&49\end{pmatrix}}=U The matrix is in row-echelon form now.

Obtain L {\displaystyle L} L by undoing your row-reduction steps. This step may be a bit tricky at first, but we are essentially constructing a matrix by going backwards. Let's look at the most recent row reduction R 3 → 3 R 3 + 4 R 2 . {\displaystyle R_{3}\to 3R_{3}+4R_{2}.} R_{{3}}\to 3R_{{3}}+4R_{{2}}. We found the new row 3 by replacing it with a linear combination of the old rows of the matrix. Now, we wish to find the old row 3, so simply solve. R 3 old → R 3 − 4 R 2 3 {\displaystyle R_{3}^{\text{old}}\to {\frac {R_{3}-4R_{2}}{3}}} R_{{3}}^{{{\text{old}}}}\to {\frac {R_{{3}}-4R_{{2}}}{3}} This undoes the second row-reduction. Now, we put it in matrix form. Let's call this matrix S 2 . {\displaystyle S_{2}.} S_{{2}}. The column vector to the right simply clarifies what we are doing - this matrix we are constructing is a linear transformation that does the same thing as what we just wrote above. Observe that, since we didn't do anything to the top two rows, the resulting elements for the two rows in this matrix are the same as in the identity matrix. Only the third row changes. ( 1 0 0 0 1 0 0 − 4 / 3 1 / 3 ) ( R 1 R 2 R 3 ) {\displaystyle {\begin{pmatrix}1&0&0\\0&1&0\\0&-4/3&1/3\end{pmatrix}}{\begin{pmatrix}R_{1}\\R_{2}\\R_{3}\end{pmatrix}}} {\begin{pmatrix}1&0&0\\0&1&0\\0&-4/3&1/3\end{pmatrix}}{\begin{pmatrix}R_{{1}}\\R_{{2}}\\R_{{3}}\end{pmatrix}} Construct the matrix that undoes the first row-reduction. Similarly, we are solving for the old row 2 and 3. We'll call this matrix S 1 . {\displaystyle S_{1}.} S_{{1}}. R 2 old = R 2 − R 1 2 {\displaystyle R_{2}^{\text{old}}={\frac {R_{2}-R_{1}}{2}}} R_{{2}}^{{{\text{old}}}}={\frac {R_{{2}}-R_{{1}}}{2}} R 3 old = R 3 + 3 R 1 2 {\displaystyle R_{3}^{\text{old}}={\frac {R_{3}+3R_{1}}{2}}} R_{{3}}^{{{\text{old}}}}={\frac {R_{{3}}+3R_{{1}}}{2}} S 1 = ( 1 0 0 − 1 / 2 1 / 2 0 3 / 2 0 1 / 2 ) {\displaystyle S_{1}={\begin{pmatrix}1&0&0\\-1/2&1/2&0\\3/2&0&1/2\end{pmatrix}}} S_{{1}}={\begin{pmatrix}1&0&0\\-1/2&1/2&0\\3/2&0&1/2\end{pmatrix}} Multiply the S {\displaystyle S} S matrices in the order that we found them. This means that S 2 S 1 = L . {\displaystyle S_{2}S_{1}=L.} S_{{2}}S_{{1}}=L. If you did the multiplication correctly, you should get a lower triangular matrix. L = ( 1 0 0 − 1 / 2 1 / 2 0 3 / 2 − 4 / 6 1 / 6 ) {\displaystyle L={\begin{pmatrix}1&0&0\\-1/2&1/2&0\\3/2&-4/6&1/6\end{pmatrix}}} L={\begin{pmatrix}1&0&0\\-1/2&1/2&0\\3/2&-4/6&1/6\end{pmatrix}}

Rewrite the matrix equation A x = b {\displaystyle A\mathbf {x} =\mathbf {b} } A{\mathbf {x}}={\mathbf {b}} in terms of L U {\displaystyle LU} LU. Now that we have both matrices, we can see below that L {\displaystyle L} L acting on the vector U x {\displaystyle U\mathbf {x} } U{\mathbf {x}} outputs b . {\displaystyle \mathbf {b} .} {\mathbf {b}}. L ( U x ) = b {\displaystyle L(U\mathbf {x} )=\mathbf {b} } L(U{\mathbf {x}})={\mathbf {b}} Since U x {\displaystyle U\mathbf {x} } U{\mathbf {x}} is a vector, let y = U x . {\displaystyle \mathbf {y} =U\mathbf {x} .} {\mathbf {y}}=U{\mathbf {x}}. Then, we see that L y = b . {\displaystyle L\mathbf {y} =\mathbf {b} .} L{\mathbf {y}}={\mathbf {b}}. The goal here is to first solve for y , {\displaystyle \mathbf {y} ,} {\mathbf {y}}, then plug into U x = y {\displaystyle U\mathbf {x} =\mathbf {y} } U{\mathbf {x}}={\mathbf {y}} to solve for x . {\displaystyle \mathbf {x} .} {\mathbf {x}}.

Solve for y {\displaystyle \mathbf {y} } {\mathbf {y}}. Because we are dealing with triangular matrices, back-substitution is the way to go. L y = ( y 1 − 1 2 y 1 + 1 2 y 2 3 2 y 1 − 2 3 y 2 + 1 6 y 3 ) = ( 2 1 7 ) {\displaystyle L\mathbf {y} ={\begin{pmatrix}y_{1}\\-{\frac {1}{2}}y_{1}+{\frac {1}{2}}y_{2}\\{\frac {3}{2}}y_{1}-{\frac {2}{3}}y_{2}+{\frac {1}{6}}y_{3}\end{pmatrix}}={\begin{pmatrix}2\\1\\7\end{pmatrix}}} L{\mathbf {y}}={\begin{pmatrix}y_{{1}}\\-{\frac {1}{2}}y_{{1}}+{\frac {1}{2}}y_{{2}}\\{\frac {3}{2}}y_{{1}}-{\frac {2}{3}}y_{{2}}+{\frac {1}{6}}y_{{3}}\end{pmatrix}}={\begin{pmatrix}2\\1\\7\end{pmatrix}} y = ( 2 4 40 ) {\displaystyle \mathbf {y} ={\begin{pmatrix}2\\4\\40\end{pmatrix}}} {\mathbf {y}}={\begin{pmatrix}2\\4\\40\end{pmatrix}}

Solve for x {\displaystyle \mathbf {x} } {\mathbf {x}}. This will again involve back-substitution, because U {\displaystyle U} U is triangular. U x = ( 2 x 1 + 4 x 2 + x 3 6 x 2 + 7 x 3 49 x 3 ) = ( 2 4 40 ) {\displaystyle U\mathbf {x} ={\begin{pmatrix}2x_{1}+4x_{2}+x_{3}\\6x_{2}+7x_{3}\\49x_{3}\end{pmatrix}}={\begin{pmatrix}2\\4\\40\end{pmatrix}}} U{\mathbf {x}}={\begin{pmatrix}2x_{{1}}+4x_{{2}}+x_{{3}}\\6x_{{2}}+7x_{{3}}\\49x_{{3}}\end{pmatrix}}={\begin{pmatrix}2\\4\\40\end{pmatrix}} x = ( 57 / 49 − 2 / 7 40 / 49 ) {\displaystyle \mathbf {x} ={\begin{pmatrix}57/49\\-2/7\\40/49\end{pmatrix}}} {\mathbf {x}}={\begin{pmatrix}57/49\\-2/7\\40/49\end{pmatrix}} Although this method may not seem very efficient to you (and indeed, LU factorization for systems of three equations is no better than row-reduction), computers are well-equipped to perform back-substitution, so the results really show as the number of equations goes up.