Skip to main content

Sage and Linear Algebra Worksheet: FCLA Section PDM

Robert Beezer
Department of Mathematics and Computer Science
University of Puget Sound
Fall 2019

Section 1 LU Decomposition, Triangular Form

This is a topic not covered in our text. You can find a discussion in A Second Course in Linear Algebra at http://linear.ups.edu/scla/html/index.html.

Our goal is to row-reduce a matrix with elementary matrices, track the changes, and arrive at an expression for a square matrix \(A\) as a product of a lower-triangular matrix, \(L\text{,}\) and an upper-triangular matrix, \(U\text{,}\) that is

\begin{equation*} A=LU \end{equation*}

the so-called LU decomposition. I sometimes prefer to call it triangular form.

There are no exercises in this worksheet, but instead there is a careful and detailed exposition of using elementary matrices (row operations) to arrive at a matrix decomposition. There are many kinds of matrix decompositions, such as the singular value decomposition (SVD). Five or six such decompositions form a central part of the linear algebra canon. Again, see A Second Course in Linear Algebra for details on these.

We decompose a \(5\times 5\) matrix. It is most natural to describe an LU decomposition of a square matrix, but the decomposition can be generalized to rectangular matrices.

Elementary matrices to “do” row operations in the first column.

Now in second column, moving to row-echelon form (i.e. not reduced row-echelon form).

The “bottom” of the third column.

And now the penultimate column.

And done.

Clearly, F has determinant 1, since it is an upper triangular matrix with diagonal entries equal to \(1\text{.}\) By tracking the effect of the above manipulations (tantamount to performing row operations) we expect that

\begin{equation*} \det(A) = \left(\frac{1}{-1/6}\right)\left(\frac{1}{-3}\right)\left(\frac{1}{-1/3}\right)\left(\frac{1}{1}\right)\left(\frac{1}{1}\right)\det(F) = -6. \end{equation*}

Let's check.

Yep. But it gets better. F is the product of the “action” matrices on the left of A.

Notice that the elementary matrices we used are all lower triangular (because we just formed zeros below the diagonal of the original matrix as we brought it to row-echelon form, and there were no row swaps). Hence their product is again lower triangular. Now check that we have the correct matrix.

The “total action” matrix is a product of elementary matrices, which are individually nonsingular. So their product is nonsingular. Futhermore, the inverse is again lower triangular.

We reach our goal by rearranging the equality above, writing A as a product of a lower-triangular matrix with an upper-triangular matrix.

Yes! So we have decomposed the original matrix (A) into the product of a lower triangular matrix (inverse of the total action matrix) and an upper triangular matrix with all ones on the diagonal (F, the original matrix in row-echelon form).

This decomposition (the LU decomposition) can be useful for solving systems quickly. You forward solve with \(L\text{,}\) then back solve with \(U\text{.}\)

More specifically, suppose you want to solve \(A\mathbf{x}=\mathbf{b}\) for \(\mathbf{x}\text{,}\) and you have a decomposition \(A=LU\text{.}\) First solve the intermediate system, \(L\mathbf{y}=\mathbf{b}\) for \(\mathbf{y}\text{,}\) which can be accomplished easily by determining the entries of \(\mathbf{y}\) in order, exploiting the lower triangular nature of \(L\text{.}\) This is what is meant by the term forward solve.

With a solution for \(\mathbf{y}\text{,}\) form the system \(U\mathbf{x}=\mathbf{y}\text{.}\) You can check that a solution, \(\mathbf{x}\text{,}\) to this system is also a solution to the original system \(A\mathbf{x}=\mathbf{b}\text{.}\) Further, this solution can be found easily by determining the entries of \(\mathbf{x}\) in reverse order, exploiting the upper triangular nature of \(U\text{.}\) This is what is meant by the term back solve.

We solve two simple systems, but only do half as many row-operations as if we went fully to reduced row-echelon form. If you count the opertions carefully, you will see that this is a big win, roughly reducing computation time by a factor of half for large systems.

This work is Copyright 2016–2019 by Robert A. Beezer. It is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.