Session 3: $LU$ Factorization

Finding the inverse of a matrix, given its factorization $PA = LU$

Although matrix inversion is not used for solving $A\mathbf{x} = \mathbf{b}$, there are a few applications where explicit knowledge of $A^{-1}$ is desirable.

We can use the $LU$ factors of a nonsingular $n$-by-$n$ matrix $A$ to compute $A^{-1}$ efficiently. The strategy is to solve the matrix equation $AX = I$. Recall that $AA^{-1} = I$ implies $A[A^{-1}]_{*j} = \mathbf{e}_j$, where $[A^{-1}]_{*j}$ is the $j^{th}$ column of $A^{-1}$, and $\mathbf{e}_j$ is the $j^{th}$ column of the identity matrix. Thus, $[A^{-1}]_{*j}$ is the solution of a system $A\mathbf{x}_j = \mathbf{e}_j$, for $j = 0, 1, ..., n-1$. Each of these $n$ systems has the same coefficient matrix, so, once the $LU$ factors for $A$ are known, each system $A\mathbf{x}_j = LU\mathbf{x}_j = \mathbf{e}_j$ can be solved by the standard two-step process:

  1. Set $\mathbf{y}_j = U\mathbf{x}_j$, and solve $L\mathbf{y}_j = \mathbf{e}_j$ for $\mathbf{y}_j$ by forward substitution.

  2. Solve $U\mathbf{x}_j = \mathbf{y}_j$ for $\mathbf{x}_j = [A^{-1}]_{*j}$ by back substitution.

The method has at least two advantages: it's efficient, and any code written to solve $A\mathbf{x} = \mathbf{b}$ can also be used to compute $A^{-1}$.

Write a function computeInverse( A ) that uses cs111.LUfactor, cs111.Lsolve, and cs111.Usolve to compute the inverse of the input matrix A by following the above algorithm. Here, however, you must take into account a permutation too. You may assume that A will always be square. Note that you should not call cs111.LUfactor more than once within your code.

Test your function on the following matrix:

$$A = \begin{pmatrix} 1 & 1 & 2 \\ 1 & 1 & 3 \\ 2 & 3 & 4 \\ \end{pmatrix} $$

by computing $AA^{-1}$ and also by comparing your result to the inverse matrix calculated with numpy's builtin function np.linalg.inv( A ). This is the same matrix we saw in lecture, and we know that the two last rows must permute so that $LU$ factorization succeeds.

Some $LU$ factorization theorems and facts

Theorem

An $n$-by-$n$ matrix $A$ has an $LU$ factorization if and only if it is invertible.

Proof:

(No partial pivoting for simplicity here) First, we prove that if $LU$ factorization succeeds, then $A$ is invertible. To do so, recall that an elementary matrix of the form $M = I + \alpha\mathbf{e}_i \mathbf{e}_j^T$ adds a multiple ($\alpha$) of row $j$ to row $i$ when it premultiplies another matrix. For example, if $\alpha = -1$ and

$$ B = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 8 & 9 & 10 \end{pmatrix}, $$

then

$$ M = I + (-1)\mathbf{e}_2\mathbf{e}_0^T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & 0 & 1 \end{pmatrix} $$

and

$$ MB = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 7 & 7 \end{pmatrix}. $$

Further, any elementary matrix of the above type (e.g., $M$) is always invertible and $M^{-1} = I - \alpha\mathbf{e}_i\mathbf{e}_j^T$. For example,

$$ M^{-1} = I - (-1)\mathbf{e}_2\mathbf{e}_0^T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix} \quad \textrm{and} \quad M^{-1}M = MM^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$

Thus, we can use elementary matrices like $M$ above to perform Gaussian elimination and transform $A$ into an upper triangular matrix $U$. If no zero pivot is found along the way (i.e., all diagonal elements in $U$ are nonzero), then

$$ M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0 A = U $$

Since any $M_\ell$ is invertible, we can write

$$ A = M_0^{-1} M_1^{-1} \cdots M_\ell^{-1} \cdots M_{m-1}^{-1} M_m^{-1}U = LU, $$

which proves that $A$ is non singular because $A^{-1} = (LU)^{-1} = U^{-1}L^{-1} = U^{-1} M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0$ exists.

Note: For the explanation involving permutations and partial pivoting, refer to https://piazza.com/class/kt26q57z3e27p3?cid=14.

To prove the converse (if $A$ is invertible, then its $LU$ factorization succeeds), recall that $A$ is invertible if $\textrm{rank}(A) = n$. This implies that all the rows in $A$ span $\mathbb{R}^n$ and no linear combination of $(n-1)$ rows can produce any other row. Therefore, the operations encoded in the elementary matrices $M$ above can never lead us to a zero row (and zero pivot, at least in exact arithmetic). This ensures that $M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0 A = U$, and Gaussian elimination successfuly factors $A$ into $L$ and $U$ (with partial pivoting).

Lemma

If Gaussian elimination (or $LU$ factorization with partial pivoting) fails, then $A$ is singular.

Proof: Suppose $A$ is an $n$-by-$n$ matrix. $LU$ factorization fails if in some step the algorithm produces a $\mathbf{0}$ row. In that case, $\textrm{rank}(A) < n$ or its nullspace has at least one nonzero basis vector. Similarly, if we are solving a system of linear equations, this means $A\mathbf{x} = \mathbf{b}$ has either no solution or infinitely many solutions depending on $\mathbf{b}$. Therefore, $A^{-1}$ does not exist.