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?
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.