LU factorization

A major tool in numerical linear algebra is to factor a given matrix into terms that are individually easier to deal with than the original. In this section we derive a means to express a square matrix using triangular factors, which will allow us to solve a linear system using forward and backward substitution.

2.4.1Outer products¶

Our derivation of the factorization hinges on an expression of matrix products in terms of vector outer products. If $\mathbf{u}\in\real^m$ and $\mathbf{v}\in\real^n$, then the outer product of these vectors is the $m\times n$ matrix

We illustrate the connection of outer products to matrix multiplication by a small example.

It is not hard to derive the following generalization of Example 2.4.1 to all matrix products.

2.4.2Triangular product¶

Equation (2.4.4) has some interesting structure for the product $\mathbf{L}\mathbf{U}$, where $\mathbf{L}$ is $n\times n$ and lower triangular (i.e., zero above the main diagonal) and $\mathbf{U}$ is $n\times n$ and upper triangular (zero below the diagonal).

Let the columns of $\mathbf{L}$ be written as $\boldsymbol{\ell}_k$ and the rows of $\mathbf{U}$ be written as $\mathbf{u}_k^T$. Then the first row of $\mathbf{L}\mathbf{U}$ is

Likewise, the first column of $\mathbf{L}\mathbf{U}$ is

These two calculations are enough to derive one of the most important algorithms in scientific computing.

2.4.3Triangular factorization¶

Our goal is to factor a given $n\times n$ matrix $\mathbf{A}$ as the triangular product $\mathbf{A}=\mathbf{L}\mathbf{U}$. It turns out that we have $n^2+n$ total nonzero unknowns in the two triangular matrices, so we set $L_{11}=\cdots = L_{nn}=1$, making $\mathbf{L}$ a unit lower triangular matrix.

We have arrived at the linchpin of solving linear systems.

The outer product algorithm for LU factorization seen in Demo 2.4.3 is coded as Function 2.4.1.

2.4.4Gaussian elimination and linear systems¶

In your first matrix algebra course, you probably learned a triangularization technique called Gaussian elimination or row elimination to solve a linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$. In most presentations, you form an augmented matrix $[\mathbf{A}\;\mathbf{b}]$ and do row operations until the system reaches an upper triangular form, followed by backward substitution. LU factorization is equivalent to Gaussian elimination in which no row swaps are performed, and the elimination procedure produces the factors if you keep track of the row multipliers appropriately.

Like Gaussian elimination, the primary use of LU factorization is to solve a linear system. It reduces a given linear system to two triangular ones. From this, solving $\mathbf{A}\mathbf{x}=\mathbf{b}$ follows immediately from associativity:

Defining $\mathbf{z} = \mathbf{U} \mathbf{x}$ leads to the following.

A key advantage of the factorization point of view is that it depends only on the matrix $\mathbf{A}$. If systems are to be solved for a single $\mathbf{A}$ but multiple different versions of $\mathbf{b}$, then the factorization approach is more efficient, as we’ll see in Efficiency of matrix computations.

As noted in the descriptions of Function 2.4.1 and Algorithm 2.4.2, the LU factorization as we have seen it so far is not stable for all matrices. In fact, it does not always even exist. The missing element is the row swapping allowed in Gaussian elimination. We will address these issues in Row pivoting.

2.4.5Exercises¶

1. ✍ For each matrix, produce an LU factorization by hand.

   (a) $\quad \displaystyle \begin{bmatrix} 2 & 3 & 4 \\ 4 & 5 & 10 \\ 4 & 8 & 2 \end{bmatrix}\qquad$ (b) $\quad \displaystyle \begin{bmatrix} 6 & -2 & -4 & 4\\ 3 & -3 & -6 & 1 \\ -12 & 8 & 21 & -8 \\ -6 & 0 & -10 & 7 \end{bmatrix}$

2. ⌨ The matrices

   $$\mathbf{T}(x,y) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ x & y & 1 \end{bmatrix},\qquad \mathbf{R}(\theta) = \begin{bmatrix} \cos\theta & \sin \theta & 0 \\ -\sin\theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

   (2.4.9)

   are used to represent translations and rotations of plane points in computer graphics. For the following, let

   $$\mathbf{A} = \mathbf{T}(3,-1)\mathbf{R}(\pi/5)\mathbf{T}(-3,1), \qquad \mathbf{z} = \begin{bmatrix} 2 \\ 2 \\ 1 \end{bmatrix}.$$

   (2.4.10)

   (a) Find $\mathbf{b} = \mathbf{A}\mathbf{z}$.

   (b) Use Function 2.4.1 to find the LU factorization of $\mathbf{A}$.

   (c) Use the factors with triangular substitutions in order to solve $\mathbf{A}\mathbf{x}=\mathbf{b}$, and find $\mathbf{x}-\mathbf{z}$.

3. ⌨ Define

   $$\mathbf{A}= \begin{bmatrix} 1 & 0 & 0 & 0 & 10^{12} \\ 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{bmatrix}, \quad \hat{\mathbf{x}} = \begin{bmatrix} 0 \\ 1/3 \\ 2/3 \\ 1 \\ 4/3 \end{bmatrix}, \quad \mathbf{b} = \mathbf{A}\hat{\mathbf{x}}.$$

   (2.4.11)

   (a) Using Function 2.4.1 and triangular substitutions, solve the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$, showing the result. To the nearest integer, how many accurate digits are in the result? (The answer is much less than the full 16 of double precision.)

   (b) Repeat part (a) with 10²⁰ as the element in the upper right corner. (The result is even less accurate. We will study the causes of such low accuracy in Conditioning of linear systems.)

4. ⌨ Let

   $$\mathbf{A} = \begin{bmatrix} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 1 & 0 & 1 \\ 1 & 0 & 0 & 1 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 1 \\ 0 & 1 & 1 & 0 & 0 & 1 \end{bmatrix}.$$

   (2.4.12)

   Verify computationally that if $\mathbf{A}=\mathbf{L}\mathbf{U}$ is the LU factorization, then the elements of $\mathbf{L}$, $\mathbf{U}$, $\mathbf{L}^{-1}$, and $\mathbf{U}^{-1}$ are all integers. Do not rely just on visual inspection of the numbers; perform a more definitive test.

5. ⌨ Function 2.4.1 factors $\mathbf{A}=\mathbf{L}\mathbf{U}$ in such a way that $\mathbf{L}$ is a unit lower triangular matrix—that is, has all ones on the diagonal. It is also possible to define the factorization so that $\mathbf{U}$ is a unit upper triangular matrix instead. Write a function lufact2 that uses Function 2.4.1 without modification to produce this version of the factorization. (Hint: Begin with the standard LU factorization of $\mathbf{A}^T$.) Demonstrate on a nontrivial $4\times 4$ example.

6. When computing the determinant of a matrix by hand, it’s common to use cofactor expansion and apply the definition recursively. But this is terribly inefficient as a function of the matrix size.

   (a) ✍ Explain using determinant properties why, if $\mathbf{A}=\mathbf{L}\mathbf{U}$ is an LU factorization,

   $$\det(\mathbf{A}) = U_{11}U_{22}\cdots U_{nn}=\prod_{i=1}^n U_{ii}.$$

   (2.4.13)

   (b) ⌨ Using the result of part (a), write a function determinant(A) that computes the determinant using Function 2.4.1. Test your function on at least two nontriangular $5\times 5$ matrices, comparing your result to the result of the standard det function.