Calculate the LU-decomposition $PA=LU$ : How do we calulate $L$?

Question

Calculate the LU-decomposition $PA=LU$ for the matrix $$A=\begin{pmatrix}3 & 1 & -3 & 2 \\ -2 & 1 & 0 & 0 \\ 2 & -2 & 4 & 1 \\ 0 & -1 & -1 & 3\end{pmatrix}$$ with column pivoting.

First we apply Gauss-elimination : \begin{align*}\begin{pmatrix}3 & 1 & -3 & 2 \\ -2 & 1 & 0 & 0 \\ 2 & -2 & 4 & 1 \\ 0 & -1 & -1 & 3\end{pmatrix} & \ \overset{R_2:R_2+\frac{2}{3}\cdot R_1}{\longrightarrow} \ \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & \frac{5}{3} & -2 & \frac{4}{3} \\ 2 & -2 & 4 & 1 \\ 0 & -1 & -1 & 3\end{pmatrix} \overset{R_3:R_3-\frac{2}{3}\cdot R_1}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & \frac{5}{3} & -2 & \frac{4}{3} \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & -1 & -1 & 3\end{pmatrix} \\ & \overset{R_2\leftrightarrow R_3}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & \frac{5}{3} & -2 & \frac{4}{3} \\ 0 & -1 & -1 & 3\end{pmatrix} \overset{R_3:R_3+\frac{5}{8}\cdot R_2}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & 0 & \frac{7}{4} & \frac{9}{8}\\ 0 & -1 & -1 & 3\end{pmatrix} \\ & \overset{R_4:R_4-\frac{3}{8}\cdot R_2}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & 0 & \frac{7}{4} & \frac{9}{8}\\ 0 & 0 & -\frac{13}{4} & \frac{25}{8}\end{pmatrix} \overset{R_3 \leftrightarrow R_4}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & 0 & -\frac{13}{4} & \frac{25}{8} \\ 0 & 0 & \frac{7}{4} & \frac{9}{8} \end{pmatrix} \\ & \overset{R_4 : R_4+\frac{7}{13}\cdot R_3}{\longrightarrow} \begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & 0 & -\frac{13}{4} & \frac{25}{8} \\ 0 & 0 & 0 & \frac{73}{26} \end{pmatrix} \end{align*} After the Gauss elimination we get \begin{equation*}U=\begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & -\frac{8}{3} & 6 & -\frac{1}{3} \\ 0 & 0 & -\frac{13}{4} & \frac{25}{8} \\ 0 & 0 & 0 & \frac{73}{26} \end{pmatrix}\end{equation*} From each step at the Gauss eliminaation we get the matrices \begin{equation*}G_1=\begin{pmatrix}1 & 0 & 0 &0 \\ -\frac{2}{3} & 1 & 0 & 0 \\ \frac{2}{3} & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{pmatrix} \ , \ G_2=\begin{pmatrix}1 & 0 & 0 &0 \\ 0 & 1 & 0 & 0 \\ 0 & -\frac{5}{8} & 1 & 0 \\ 0 & \frac{3}{8} & 0 & 1\end{pmatrix} \ \text{ and } \ G_3=\begin{pmatrix}1 & 0 & 0 &0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -\frac{7}{13} & 1\end{pmatrix}\end{equation*} We have the permutation matrices \begin{equation*}P_0=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \ \text{ und } \ P_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{pmatrix}\end{equation*} We have that\begin{equation*}P=P_1\cdot P_0=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \end{pmatrix}\end{equation*} Now I am confused about how to calculate the matrix $L$.

Is it $L= P\cdot (G_3\cdot P_1\cdot G_2\cdot P_0\cdot G_1)^{-1}$ or $L= (G_3\cdot P_1\cdot G_2\cdot P_0\cdot G_1)^{-1}$ or somehow else?

Answer 1

With an algebraic calculator, we have

A1:rowop(A,2,1,-2/3)

A2:rowop(A1,3,1,2/3)

A3:rowop(A2,3,2,-8/5)

A4:rowop(A3,4,2,-3/5)

A5:rowop(A4,4,3,-11/14)

That gives the matrices:

$G_1=\begin{pmatrix}1 & 0 & 0 & 0\\2/3 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{pmatrix}$

$G_2=\begin{pmatrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\-2/3 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{pmatrix}$

$G_3=\begin{pmatrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 8/5 & 1 & 0\\0 & 0 & 0 & 1\end{pmatrix}$

$G_4=\begin{pmatrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 3/5 & 0 & 1\end{pmatrix}$

$G_5=\begin{pmatrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 11/14 & 1\end{pmatrix}$

We have $G_5 G_4 G_3 G_2 G_1 A = U$

So $L=(G_5 G_4 G_3 G_2 G_1)^{-1}$

Answer 2

You have the row operations right, but you store in $G_1$ the opposite of the row coefficients, that's why you end up with everything inverted ($G_1^{-1}$, etc.). But as you will see, it's the correct approach.

Let's start from $A$ and let's forget about row swaps for a moment.

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

The row ops will be $R_2\leftarrow R_2+\frac23R_1$, $R_3\leftarrow R_3-\frac23R_1$, $R_4\leftarrow R_4+0R_1$ (the last is useless but let's keep it for completeness).

These row operations are exactly equivalent to a left matrix product $A\leftarrow M_1A$ with

$$M_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ \frac23 & 1 & 0 & 0 \\ -\frac23 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{pmatrix}$$

See how the signs in the first column are reversed compared to $G_1$. And notice that the inverse of $M_1$ is obtained by reversing the signs of the coefficients of the first column, below the diagonal. That is, $M_1^{-1}=G_1$.

So, when computing $M_1A$, the first column is annihilated below the pivot:

$$M_1A=\begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & \frac53 & -2 & \frac43 \\ 0 & -\frac83 & 6 & -\frac13 \\ 0 & -1 & -1 & 3\end{pmatrix}$$

Now, let's do the second column: $R_3\leftarrow R_3+\frac85 R_2$, $R_4\leftarrow R_4+\frac35R_1$, and the matrix $M_2$ that describes these row operations is

$$M_2=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & \frac85 & 1 & 0 \\ 0 & \frac35 & 0 & 1\end{pmatrix}$$

And

$$M_2M_1A=\begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & \frac53 & -2 & \frac43 \\ 0 & 0 & \frac{14}5 & \frac95 \\ 0 & 0 & -\frac{11}5 & \frac{19}5\end{pmatrix}$$

Same remark as above about the inverse of $M_2$: change the signs of the coefficients. However, since we didn't swap rows, we don't have the same coefficients up to sign as in your $G_2$, now.

Finally, we need the row operation $R_4\leftarrow R_4+\frac{11}{14}R_3$, and the matrix:

$$M_3=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{11}{14} & 1\end{pmatrix}$$

$$M_3M_2M_1A=\begin{pmatrix}3 & 1 & -3 & 2 \\ 0 & \frac53 & -2 & \frac43 \\ 0 & 0 & \frac{14}5 & \frac95 \\ 0 & 0 & 0 & \frac{73}{14}\end{pmatrix}$$

So far, we have the decomposition:

$$M_3M_2M_1A=U$$

Now, a product of lower triangular matrices with ones on the diagonal has the same form, so $M_3M_2M_1$ is lower triangular, and its inverse is also lower triangular, so we have the decomposition:

$$A=(M_3M_2M_1)^{-1}U=M_1^{-1}M_2^{-1}M_3^{-1}U=G_1G_2G_3U$$

(here $G_2$ and $G_3$, while similar to what you did, are different since we didn't swap rows, but I keep your notation, as I think it makes things clearer)

That is, $A=LU$ with $L=G_1G_2G_3$. Notice how it's similar to your case, except the permutations matrices $P_0,P_1$ do not appear. Theoretically these permutation matrices are only required when you get a zero pivot. For numerical computations it's better to always pick the largest pivot to reduce numerical errors.

But we are not ready to permute rows yet.

This product $G_1G_2G_3$ has an interesting property, that I'll show on an example. Let:

$$T_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ u_2 & 1 & 0 & 0 \\ u_3 & 0 & 1 & 0 \\ u_4 & 0 & 0 & 1\end{pmatrix}$$

$$T_2=\begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & v_3 & 1 & 0 \\ 0 & v_4 & 0 & 1\end{pmatrix}$$

Then

$$T_2T_1=\begin{pmatrix}1 & 0 & 0 & 0 \\ u_2 & 1 & 0 & 0 \\ u_3+v_3u_2 & v_3 & 1 & 0 \\ u_4+v_4u_2 & v_4 & 0 & 1\end{pmatrix}$$

But

$$T_1T_2=\begin{pmatrix}1 & 0 & 0 & 0 \\ u_2 & 1 & 0 & 0 \\ u_3 & v_3 & 1 & 0 \\ u_4 & v_4 & 0 & 1\end{pmatrix}$$

See? The product $T_2T_1$ is relatively difficult to compute, whereas the product $T_1T_2$ is trivial: just put the coefficients together! It's a general result when you multiply these transformation matrices: if the product starts with the leftmost columns and goes to the right, you just put the coefficient together. Otherwise you have to compute the matrix product.

Now remember that $L=G_1G_2G_3$: that is, we just need to put the coefficients together to write $L$. That is,

$$L=\begin{pmatrix}1 & 0 & 0 & 0 \\ -\frac23 & 1 & 0 & 0 \\ \frac23 & -\frac85 & 1 & 0 \\ 0 & -\frac35 & -\frac{11}{14} & 1\end{pmatrix}$$

So, if no row swap is involved, we have an easy way to compute the $LU$ decomposition: write the modified $A$ on one side, the $L$ coeficients on the other, and you can build $L$ and $U$ during the computation. We don't even need to write them apart: when $A$ is transformed to $U$, only the upper triangular part is useful, so the lower triangular array can be used to store $L$. This way, you don't need extra space to compute $L$ and $U$: just overwrite $A$. It's how $LU$ is implemented in practice (with more subtleties if you want to be cache-friendly, but it's another story).

Hence, we compute the $LU$ decomposition as:

$$\begin{pmatrix}3 & 1 & -3 & 2 \\ \color{red}{-\frac23} & \frac53 & -2 & \frac43 \\ \color{red}{\frac23} & \color{red}{-\frac85} & \frac{14}5 & \frac95 \\ \color{red}{0} & \color{red}{-\frac35} & \color{red}{-\frac{11}{14}} & \frac{73}{14}\end{pmatrix}$$

In black: the $U$ part. In red: the $L$ part. The diagonal of $L$ is missing, but we know it's all made of ones.

Now, what happens it swaps intervene? Surely, it must be complicated, with the storage in $A$ not working anymore.

Let's see.

We are doing the same process, but we will choose each time the largest pivot. For the first column, it's already in place, so the first permutation matrix is $P_1=I$, and we compute $M_1P_1A$ as previously. However, we will now store $L$ inside the same matrix. Since this completed matrix is not really equal to the product, I write the formula in brackets to remember we are doing something more.

$$[M_1P_1A]=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{-\frac23} & \frac53 & -2 & \frac43 \\ \color{red}{\frac23} & -\frac83 & 6 & -\frac13 \\ \color{red}{0} & -1 & -1 & 3 \end{pmatrix}$$

In red, the $L$ part, in black, the $A$ part that it transformed to $U$. Remember that the $L$ part has the coefficients of $M_1^{-1}$, that is the coefficients of $G_1$ (signs swapped).

Now the second column: we would like to take $-\frac83$ as pivot, so we need to swap the rows in the black part. How do we keep track of this in the $L$ part? Simple. Swap the $L$ part too!

So, first the swap:

$$P_2=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}$$

$$[P_2M_1P_1A]=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{\frac23} & -\frac83 & 6 & -\frac13 \\ \color{red}{-\frac23} & \frac53 & -2 & \frac43 \\ \color{red}{0} & -1 & -1 & 3 \end{pmatrix}$$

Note how the black bart is identical to what you actually did. Now the pivot, with

$$M_2=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & \frac58 & 1 & 0\\ 0 & -\frac38 & 0 & 1 \end{pmatrix}$$

Now apply this transformation. Wait, to the $A$ part, or also to the $L$ part? Now it's only the $A$ part: we want to swap the rows in $L$ to keep track of the permutation, we don't want to change the values.

$$[M_2P_2M_1P_1A]=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{\frac23} & -\frac83 & 6 & -\frac13 \\ \color{red}{-\frac23} & \color{red}{-\frac58} & \frac74 & \frac98 \\ \color{red}{0} & \color{red}{\frac38} & -\frac{13}4 & \frac{25}8 \end{pmatrix}$$

Let's do it again with the fourth column. This time I don't write $P_3$ and $M_3$, but I write the two steps separately to make it clearer.

First swap rows $3$ and $4$:

$$[P_3M_2P_2M_1P_1A]=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{\frac23} & -\frac83 & 6 & -\frac13 \\ \color{red}{0} & \color{red}{\frac38} & -\frac{13}4 & \frac{25}8 \\ \color{red}{-\frac23} & \color{red}{-\frac58} & \frac74 & \frac98 \end{pmatrix}$$

Now $R_4\leftarrow R_4+\frac7{13}R_3$:

$$[P_3M_2P_2M_1P_1A]=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{\frac23} & -\frac83 & 6 & -\frac13 \\ \color{red}{0} & \color{red}{\frac38} & -\frac{13}4 & \frac{25}8 \\ \color{red}{-\frac23} & \color{red}{-\frac58} & \color{red}{-\frac{7}{13}} & \frac{73}{26} \end{pmatrix}$$

And we are done. We can write down the matrices $L$ and $U$:

$$L=\begin{pmatrix} 1 & 0 & 0 & 0 \\ \color{red}{\frac23} & 1 & 0 & 0 \\ \color{red}{0} & \color{red}{\frac38} & 1 & 0 \\ \color{red}{-\frac23} & \color{red}{-\frac58} & \color{red}{-\frac{7}{13}} & 1 \end{pmatrix}$$

$$U=\begin{pmatrix} 3 & 1 & -3 & 2 \\ \color{red}{0} & -\frac83 & 6 & -\frac13 \\ \color{red}{0} & \color{red}{0} & -\frac{13}4 & \frac{25}8 \\ \color{red}{0} & \color{red}{0} & \color{red}{0} & \frac{73}{26} \end{pmatrix}$$

And now, $PA=LU$, with $P=P_3P_2P_1$.

In practice, $L$ and $U$ are computed in place, in $A$, and the matrix $P$ is not stored, but you store enough information to know the permutation matrix $P$. Since there is a bijection between permutation matrices and permutations of integers, you can simply store the integer vector $P=[1,2,\dots,n]$, and swap elements in $P$ when you swap rows. It's not the only possible storage: if I remember correctly, the storage in LAPACK is a bit different to allow more efficient computations with $L$ and $U$.

For instance, to compute the determinant of $A$, you compute the product of diagonal elements of $U$, and multiply by $\sigma(\pi)$, the parity of the permutation $\pi$. Even though it's not difficult to compute, LAPACK's storage allows for a faster computation of this parity. See https://stackoverflow.com/questions/47315471/compute-determinant-from-lu-decomposition-in-lapack