Steps for a Permuted LU Decomposition

As we know, we can use the scipy.linalg.lu command to find the permuted LU decomposition of a matrix. How can we show the steps involved in that process? Can Python help with that?

• edited January 22

Of course we can do this! I'll do so using my personal matrix from the Your personal LU Decomposition problem:

$$A = \left( \begin{array}{cccc} 2 & 4 & 4 & 0 \\ 0 & -1 & 3 & -1 \\ 8 & 16 & 18 & 2 \\ -2 & -3 & 1 & 12 \\ \end{array} \right).$$

The steps shown below are very similar to those shown on pages 28-30 of our textbook. The one distinction is that we are using the numerical pivoting strategy described on pages 55-56 of our text, which suggest that we should always pivot to use the largest term in absolute value below the diagonal in each column.

Note that we can use Python to assist with the computations as shown in this Sage cell.

$$A = \left( \begin{array}{cccc} 2 & 4 & 4 & 0 \\ 0 & -1 & 3 & -1 \\ 8 & 16 & 18 & 2 \\ -2 & -3 & 1 & 12 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right).$$

Swapping the first and third rows yields

$$A = \left( \begin{array}{cccc} 8 & 16 & 18 & 2 \\ 0 & -1 & 3 & -1 \\ 2 & 4 & 4 & 0 \\ -2 & -3 & 1 & 12 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right).$$

Then, eliminating the terms in the first column, I get:

$$A = \left( \begin{array}{cccc} 8 & 16 & 18 & 2 \\ 0 & -1 & 3 & -1 \\ 0 & 0 & -0.5 & -0.5 \\ 0 & 1 & 5.5 & 12.5 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1/4 & 0 & 1 & 0 \\ -1/4 & 0 & 0 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right).$$

Eliminating the term in the second column yields:

$$A = \left( \begin{array}{cccc} 8 & 16 & 18 & 2 \\ 0 & -1 & 3 & -1 \\ 0 & 0 & -0.5 & -0.5 \\ 0 & 0 & 8.5 & 11.5 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1/4 & 0 & 1 & 0 \\ -1/4 & 1 & 0 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right).$$

Looks like we need another pivot. Note that this affects both $L$ and $A$:

$$A = \left( \begin{array}{cccc} 8 & 16 & 18 & 2 \\ 0 & -1 & 3 & -1 \\ 0 & 0 & 8.5 & 11.5 \\ 0 & 0 & -0.5 & -0.5 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -1/4 & 0 & 1 & 0 \\ 1/4 & 1 & 0 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ \end{array} \right).$$

Finally, we eliminate the last term in the third column:

$$A = \left( \begin{array}{cccc} 8 & 16 & 18 & 2 \\ 0 & -1 & 3 & -1 \\ 0 & 0 & 8.5 & 11.5 \\ 0 & 0 & 0 & 0.17647059 \\ \end{array} \right), L = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -1/4 & 0 & 1 & 0 \\ 1/4 & 1 & -0.5/8.5 & 1 \\ \end{array} \right), P = \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ \end{array} \right).$$