LU decomposition; do permutation matrices commute?

Question

I have an assignment for my Numerical Methods class to write a function that finds the PA=LU decomposition for a given matrix A and returns P, L, and U.

Nevermind the coding problems for a moment; there is a major mathematical problem I'm having that my professor seems incapable of addressing, and after searching for hours (perhaps inefficiently) I could find no accounting of it. How do we extricate the permutation matrices from the row elimination matrices?

Essentially the idea, if I understand it correctly, is that we perform a series of transformations on a matrix $A$ by applying successive lower triangular matrices that eliminate single elements, thus $L_nL_{n-1}...L_2L_1A = U$. In my understanding, this is computationally useful because lower triangular atomic matrices can be inverted by changing the sign of the off-diagonal element, so $A = L_1^{-1}L_2^{-1}...U$.

That's all fine (assuming I'm correct), but the introduction of pivot matrices between each $L_j$ seems to make the problem intractable. In every accounting I've seen some sorcery occurs that looks like this: $$L_nP_nL_{n-1}...P_3L_2P_2L_1P_1A = U \Rightarrow P_nP_{n-1}...P_2P_1L_nL_{n-1}...L_2L_1A = U$$

And no one bothers to explain how this happens or in fact even states it explicitly.

If possible I would like to know
a) Is this operationally acceptable?

b) What properties of these respective classes of matrices make this kind of willy-nilly commutation legal?

c) Is my understanding of the method and its advantages accurate?

Answer 1

No, in general permutation matrices do not commute.

It seems like you are using the Doolittle algorithm, so I am going to assume that indeed you are.

Let $U_i$ be the $i$:th step in your LU decomposition of $A$, so $$\begin{align} U_0 &= A \\ U_1 &= L_1P_1U_0 \\ \vdots \\ U_n &= L_nP_nU_{n-1} \end{align}$$

This corresponds well to how one would do LU-decomposition by hand; get the largest element as the leading element on the row you are at (i.e. multiply with $P_k$), then eliminate that column on the following rows (i.e. multply with $L_k$).

As you remark, the $L_i$ will be atomic lower triangular matrices, the non-zero elements all being in column $i$. The inverse of $L_i$ can be constructed by negating all off-diagonal elements, see Wikipedia.

The permutation matrix $P_j$ will be a permutation matrix switching row $j$ with a row $k \geq j$, if multiplied on the left of the matrix you want to transform. The inverse to $P_j$ is $P_j$ itself (since $P_j$ switches row $j$ with row $k$, you can undo this by doing the same thing).

Consider the product $P_jL_i$ for $i < j$. $P_j$ will switch two rows of $L_i$, row $j$ and $k \geq j > i$. We switch elements in $L_i$ as follows: $$\begin{align} L_{j,i} &\leftrightarrow L_{k,i} \\ L_{j,j} &\leftrightarrow L_{k,j} \end{align}$$ Let $L_k'$ be the matrix produced by switching just the off-diagonal elements (the first switch above). Note that this is still an atomic lower triangular matrix. We can then produce $P_jL_i$ by just switching column $j$ with column $k$ in $L_k'$, which is multiplying by $P_j$ on the right. Here it is important that $i < j, k$, so column $i$ in $L_i'$ is not changed! In other words: $$P_j L_i = L_i' P_j$$

Thus, you can from your equation $$L_nP_nL_{n-1}...P_3L_2P_2L_1P_1A = U$$ produce $$L_n^S L_{n-1}^S \cdots L_1^S P_n P_{n-1} \cdots P_1A = U$$ which can be transformed to (note that $L_i^S$ is still atomic lower triangular): $$PA = LU$$ taking $$P = P_nP_{n-1} \cdots P_1$$ and $$L = (L_n^S L_{n-1}^S \cdots L_1^S)^{-1}.$$

Here, $L_i^S$ is the matrix made by taking $L_i$ and applying all $P_j$ (on the left) for $j > i$ on the off-diagonal elements.

You do this by doing the following, starting with $$L_nP_nL_{n-1}...P_3L_2P_2L_1P_1A = U$$ move the $P_2$ matrix to the right using $P_2L_1 = L_1' P_2$, producing: $$L_nP_nL_{n-1}...P_3L_2L_1'P_2P_1A = U$$ then do the same for the matrix $P_3$, but this matrix you have to move to the right twice, using $P_3L_2 = L_2'P_3$ and $P_3L_1' = L_1''P_3$: $$L_nP_nL_{n-1}...L_2'L_1''P_3P_2P_1A = U$$ and so on for every $P_j$.

As an example of $P_j L_i = L_i' P_j$ consider $$P = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}$$ $$L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 0 & 1 \end{pmatrix}$$ $$L' = \begin{pmatrix} 1 & 0 & 0 \\ 3 & 1 & 0 \\ 2 & 0 & 1 \end{pmatrix}$$

$$PL = \begin{pmatrix} 1 & 0 & 0 \\ 3 & 0 & 1 \\ 2 & 1 & 0 \end{pmatrix}$$ $$L'P = \begin{pmatrix}1 & 0 & 0 \\ 3 & 0 & 1 \\ 2 & 1 & 0 \end{pmatrix}$$

So, to summarize: These special matrices almost commute, only small changes are needed to swap the matrices. However, all important properties of the matrices are preserved are when swapping them.

Answer 2

I found an answer based on matrix algebra rather than components here. Here's my attempt to reproduce it with some more detail.

Starting from $$L_3 P_3 L_2 P_2 L_1 P_1 A = U$$ Define $$\Lambda_1 = P_3 P_2 L_1 P_2^{-1} P_3^{-1}$$ $$\Lambda_2 = P_3 L_2 P_3^{-1}$$ $$\Lambda_3 = L_3$$

Calculating $$\Lambda_3 \Lambda_2 \Lambda_1 P_3 P_2 P_1 = (L_3)(P_3 L_2 P_3^{-1})(P_3 P_2 L_1 P_2^{-1} P_3^{-1}) P_3 P_2 P_1 = L_3 P_3 L_2 P_2 L_1 P_1$$ We can swap this into our starting formula to get $$\Lambda_3 \Lambda_2 \Lambda_1 P_3 P_2 P_1 A = U$$ And letting $\Lambda = \Lambda_1^{-1} \Lambda_2^{-1} \Lambda_3^{-1}$ and $P = P_3 P_2 P_1$ we have $P A = \Lambda U$.

Now it remains to show that the $\Lambda_i$ are lower triangular and to answer your question about commuting head on. Notice that the $P_i$ only affect rows $i$ and below because the elimination algorithm has completed its work with the rows above $i$. The inverses of the $P_i$ will therefore also only need to affect those rows to unscramble things. Examining $\Lambda_1 = P_3 P_2 L_1 P_2^{-1} P_3^{-1}$ we see that $\Lambda_1$ swaps around rows 2 and below, applies $L_1$ to add multiples of the first row to those lower rows, then unswaps the rows. This is equivalent to adding appropriate multiples of $L_1$ to the original, unswapped rows, which we know is an operation described by a lower triangular matrix. Therefore the matrix $\Lambda_1$ which represents the aggregate operation $P_3 P_2 L_1 P_2^{-1} P_3^{-1}$ must be lower triangular.

In summary, the matrices don't commute as you had asked in your original question, but we can change an $L_i$ into a related vector $\Lambda_i$ by left- and right-multiplying by some permutation matrices. Because the commuting fails, our formula for $\Lambda$ isn't particularly pretty: $$\Lambda = P_3 P_2 L_1^{-1} P_2^{-1} L_2^{-1} P_3^{-1} L_3^{-1}$$ However, we can get a neater conclusion if we return to $$L_3 P_3 L_2 P_2 L_1 P_1 A = U$$ and let $$M = L_3 P_3 L_2 P_2 L_1 P_1$$ so that $M A = U$. Then because $P A = \Lambda U$ we must have $P = \Lambda M$ or $\Lambda = P M^{-1}$. So we can calculate $\Lambda$ from our original $M$ just by applying the composition of our permutations to $M^{-1}$, which isn't quite as good as commuting but isn't half bad either.

Answer 3

$\require{color}$

I found @Calle explanation interesting but have a hard time following the arguments. @kuzooroo explanation fits better my way of thinking but the arguments are too dense. I decided to start a new answer since my arguments will not fit in a comment slot.

Following @Calle: Note that, since each $P_i$ is a one-cycle permutation (one interchange), $P_i^{-1}=P_i$ since doing an interchange twice bring the matrix to the same original form. We now should prove that $\Lambda$ is a lower triangular matrix with ones in the diagonal.

$\Lambda_3$ is, of course, a lower triangular with 1 in the diagonal. To prove that $\Lambda_2$ is lower triangular with ones in the diagonal we take into account the fact that $L_2$ is atomic. That is,

\begin{eqnarray*} L_2 = \begin{pmatrix} 1 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ 0 & 1 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ \vdots & l_{32} & \ddots & \vdots & \vdots &\vdots & \vdots & \vdots & \vdots\\ \vdots & \vdots & \cdots & 1 & 0 & \cdots &\cdots& \cdots & \vdots \\ \vdots & l_{i2} & \cdots & 0 & 1 & 0 &\cdots & \cdots & \vdots \\ \vdots & \vdots & \cdots & \vdots & \ddots & \ddots & \ddots & \cdots & \vdots \\ \vdots & l_{k2} & \vdots & \vdots & \vdots & 0 & 1 & 0 & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\ddots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & 0 & 1 \end{pmatrix} \nonumber. \end{eqnarray*}

Let us assume that $P_3$ interchanges rows $i$ with $k$. (note that both $i$ and $k$ are larger than 2). We have that

\begin{eqnarray*} P_3 L_2 = \begin{pmatrix} 1 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ 0 & 1 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ \vdots & l_{32} & \ddots & \vdots & \vdots &\vdots & \vdots & \vdots & \vdots\\ \vdots & \vdots & \cdots & 1 & 0 & \cdots &\cdots& \cdots & \vdots \\ \vdots & \boxed{l_{k2}} & \vdots & \vdots & \textcolor{green}{0} & 0 & \textcolor{blue}{1} & 0 & \vdots \\ \vdots & \vdots & \cdots & \vdots & \ddots & \ddots & \ddots & \cdots & \vdots \\ \vdots & \boxed{l_{i2}} & \cdots & 0 & \textcolor{green}{1} & 0 & \textcolor{blue}{0} & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\ddots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & 0 & 1 \end{pmatrix} \nonumber. \end{eqnarray*}

Finally $P_3 L_2 P_3^{-1}= P_3 L_2 P_3$ ineterchanges columns $i$ (green) with $k$ (blue) on $P_3 L_2$. The $1$ in column $k$ (blue) takes the place of the 0 in column $i$ (green) . In the same way the $1$ in green in row $k$ takes the position of the blue $0$ on the same row.

That is

\begin{eqnarray*} P_3 L_2 P_3^{-1} = \begin{pmatrix} 1 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ 0 & 1 & \cdots & 0 & \cdots & \cdots &\cdots & \cdots & 0 \\ \vdots & l_{32} & \ddots & \vdots & \vdots &\vdots & \vdots & \vdots & \vdots\\ \vdots & \vdots & \cdots & 1 & 0 & \cdots &\cdots& \cdots & \vdots \\ \vdots & \boxed{l_{k2}} & \vdots & \vdots & \textcolor{green}{1} & 0 & \textcolor{blue}{0} & 0 & \vdots \\ \vdots & \vdots & \cdots & \vdots & \ddots & \ddots & \ddots & \cdots & \vdots \\ \vdots & \boxed{l_{i2}} & \cdots & 0 & \textcolor{green}{0} & 0 & \textcolor{blue}{1} & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\ddots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & \cdots & \cdots &\cdots & 0 & 1 . \end{pmatrix} \nonumber. \end{eqnarray*}e

This is, again, an atomic matrix with the entries $l_{i2}$ and $l_{k2}$ interchanged. In the same way, $P_2 L_1 P_2^{-1}$ is atomic and so it is $\Lambda_1 = P_3 P_2 L_1 P_2^{-1} P_3^{-1}$,

In general for

\begin{eqnarray*} U= L_{n-1} P_{n-1} L_{n-2} P_{n-2} \cdots P_2 L_1 P_1 A, \end{eqnarray*} we define

\begin{eqnarray*} \Lambda_i = P_{n-1} \cdots P_{i+2} P_{i+1} L_i P_{i+1}^{-1} P_{i+2}^{-1} \cdots P_{n-1}^{-1}. \end{eqnarray*} where we assume $P_0=I$. The same procedure shown above guarantee that $\Lambda_i$ is atomic and so $\Lambda=\prod \Lambda_i$ is lower triangular with ones in the diagonal.

Hence this completes what is needed to verify that $PA=LU$.