Matrix factorization

Let $L$ be a lower-triangular matrix ($n \times n$) and $U$ be an upper-triangular matrix (also $n \times n$).

We've discussed how to solve $Ux = y$ with backward substitution:

$$x_i = \frac{\displaystyle y_i - \sum_{j = i+1}^n u_{ij}x_j}{u_{ii}},\quad i = n, n-1, \ldots, 2,1.$$

We also discussed how to solve $Ly = b$ with forward substitution:

$$y_i = \frac{\displaystyle b_i - \sum_{j=1}^{i-1} l_{ij} y_j }{l_{ii}}, \quad i = 1,2,\ldots,n-1,n.$$

Using an LU factorization

Assume that we can express $A = LU$ and we want to solve $Ax=b$:

$$A x = LU x = b = L(Ux) = b.$$

1. Solve for $y = Ux$ which is given by $Ly = b$ by forward substitution.
2. Solve $Ux = y$ for $x$, because $U$ and $y$ are now known by backward substitution.

Operation count

Recall from Lecture 9 that backward substitution requires

$$\frac{n^2+3n}{2} \quad \text{multiplications, and}\\ \frac{n^2+n}{2} \quad \text{additions.}$$

The same is true of forward substitution.

To solve $Ax = b$ when a factorization $A = LU$ is known requires

$${n^2+3n} \quad \text{multiplications, and}\\ {n^2+n} \quad \text{additions.}$$

This should be compared with Gaussian elimination with backward substitution (or worse, matrix inversion) that requires $\sim n^3/3$ operations.

Even if you could compute $A^{-1}$ for free, computing $x = A^{-1}b$ requires (again, see Lecture 9) $n^2$ multiplications and $n^2-n$ additions. This is nearly the same number of operations as the $LU$ method.

Having an LU factorization of a matrix is just as good (if not better) as having the inverse matrix

The LU factorization algorithm

Assume $A$ is an $n\times n$ matrix that can be put in upper-triangular form without row swaps (we'll deal with those next lecture).

Consider the $3\times 3$ case first

$$A = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$ $$R_2 - \frac{a_{21}}{a_{11}}R_1 \to R_1$$

$$R_2 - \frac{a_{21}}{a_{11}}R_1 \to R_1$$ $$A = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

Here $a^{(1)}_{22} = a_{22} - a_{12} \frac{a_{21}}{a_{11}}$ and $a^{(1)}_{23}= a_{23} - a_{13} \frac{a_{21}}{a_{11}}$. The question that will lead us to an $LU$ factorization is: What type of matrix does this row operation correspond to?

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$ $$\begin{bmatrix} 1 & 0 & 0 \\ -a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

Recall that the goal is to express $A = LU$, so we find the inverse of the lower-triangular matrix on the left-hand side:

$$\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ -a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

So, after a single step, we have found

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

This is closer to a lower/upper factorization.

If we perform the row operation $R_{3} - \frac{a_{31}}{a_{11}} R_1 \to R_3$ which corresponds to

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ - a_{31}/a_{11} & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix}$$

where $a^{(1)}_{32} = a_{32} - a_{12} \frac{a_{31}}{a_{11}}$ and $a^{(1)}_{33} = a_{33} - a_{13} \frac{a_{31}}{a_{11}}$

The inverse of the lower-triangular factor can be confirmed

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ - a_{31}/a_{11} & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0& 0 & 1 \end{bmatrix}$$

We use this with our expression from the first step

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix}$$ $$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$$

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix}$$ $$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix}$$

This brings us very close to an $LU$ factorization

For the last step, we need to perform $R_3 - \frac{a^{(1)}_{32}}{a^{(1)}_{22}} R_2 \to R_3$

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & 0 & a^{(2)}_{33} \end{bmatrix}$$

where $a^{(2)}_{33} = a^{(1)}_{33} - a^{(1)}_{23} \frac{a^{(1)}_{32}}{a^{(1)}_{22}}$.

Again, you can check that the inverse of the matrix on the left is simple

$$\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$ $$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & 0 & a^{(2)}_{33} \end{bmatrix}$$

Inserting this into the existing factorization: $$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \end{bmatrix}$$

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ a_{31}/a_{11} & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & 0 & a^{(2)}_{33} \end{bmatrix}$$

This inner factors simplify and we have an LU factorization

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} =\begin{bmatrix} 1 & 0 & 0 \\ a_{21}/a_{11} & 1 & 0 \\ a_{31}/a_{11} & a^{(1)}_{32}/a^{(1)}_{22} & 1 \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\ 0 & 0 & a^{(2)}_{33} \end{bmatrix}$$

Example

Compute the $LU$ factorization of

$$A = \begin{bmatrix} 1 & -1 & 1 \\ 2 & -3 & 4 \\ 1 & 1 & 1 \end{bmatrix}.$$

The LU algorithm (no row interchanges)

INPUT a n x n matrix A

OUTPUT the LU factorization of A, or a message of failure

STEP 1: Set L to be the n x n identity matrix;
STEP 2: For j = 1, 2, ...,n do STEPS 3-4
  STEP 3: If A(j,j) = 0
            OUTPUT('LU Factorization failed')
            STOP.
  STEP 4: For i = j+1, j+2, ... , n do STEPS 5-6
    STEP 5: Set L(i,j) = A(i,j)/A(j,j);
    STEP 6: Perform row operation Ri - L(i,j)*Rj --> Ri on A;
STEP 7: OUTPUT([L,A]);

function [L,A] = LU(A)
   % this function will give two output matrices
end

A = rand(4)
[L,U] = LU(A);
L
U
norm(A-L*U)

A =

    0.7922    0.8491    0.7431    0.7060
    0.9595    0.9340    0.3922    0.0318
    0.6557    0.6787    0.6555    0.2769
    0.0357    0.7577    0.1712    0.0462


L =

    1.0000         0         0         0
    1.2112    1.0000         0         0
    0.8277    0.2554    1.0000         0
    0.0451   -7.6181  -21.9384    1.0000


U =

    0.7922    0.8491    0.7431    0.7060
   -0.0000   -0.0944   -0.5078   -0.8233
    0.0000         0    0.1701   -0.0972
   -0.0000         0         0   -8.3903


ans =

   2.3420e-15