Struggling to infer the pseudocode of the LU factorization from the algorithm.

The algorithm that you mention at the top is the explicit solution of a system of linear equations which is illustrated in an example in Wikipedia.

The pseudo code is based on a different idea, it uses the Gaussian elimination method explained in Wikipedia. Gaussian elimination is a fast and efficient algorithm and is usually preferred.

Gaussian elimination relies on the following facts:

• Each step in performing linear combinations of the kth row with every row from the (k + 1)th to the nth is a multiplication (on the left) by an elementary row operation matrix ($E$) that is lower triangular.
• The product of two lower triangular matrices is lower triangular.
• The inverse of a lower triangular matrix is lower triangular.

So what the pseudo code is doing is to convert the original matrix $A$ to an upper triangular matrix ($U$) by multiplying $A$ on the left by a series of lower triangular elementary row operation matrices ($E_1, \dots E_n$). We start with $A$, then compute $E_1A$, then $E_2E_1A$ and so on. After multiplying by $E_k$, the first $k$ columns have become upper triangular. At the end, the entire matrix is upper triangular. In other words, $U=E_n\dots E_1 A$. Multiplying both sides on the left by the inverse of $E_n\dots E_1$, we get $\left(E_1^{-1}\dots E_n^{-1}\right) U=A$. Since products and inverses of lower triangular matrices are lower triangular, it follows that $E_1^{-1}\dots E_n^{-1}$ is a lower triangular matrix $L$. This gives us the factorization $A=LU$.

Note that not all matrices have an $LU$ factorization because the pseudo code sometimes works only with an additional step of permuting rows of the matrix. So in general we only get a decomposition $LPU$ where $P$ is a permutation matrix. And if we want both $L$ and $U$ to have only unity in the diagonal, then we get the decomposition $LPDU$ where $D$ is a diagonal matrix.

Let us observe how to perform LU factorization on a $3\times 3$ matrix. $$A=\begin{pmatrix} \times & \times & \times \\ \times & \times & \times \\ \times & \times & \times \end{pmatrix} .$$

First we multiply $A$ be a matrix $L_1$ on the left to transform it to the following shape: $$ L_1 A= \begin{pmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & \times & \times \end{pmatrix} . $$ Note that this is equivalent to multiply the first row of $A$ by $-\frac{a_{21}}{a_{11}}$ and adding it to the $2$nd row, and then multiply the first row of $A$ by $-\frac{a_{31}}{a_{11}}$ and adding it to the $3$rd row. If we take $L$ to be the corresponding elementary row transformation matrix, then $$ L_1= \begin{pmatrix} 1 & & \\ -\frac{a_{21}}{a_{11}} & 1 & \\ -\frac{a_{31}}{a_{11}} & & 1 \end{pmatrix} $$ with omitted entry being zero.

Now denote $B= L_1 A$. Let $\tilde B$ be the $2\times 2$ matrix consisting of the right-lower $2\times 2$ submatrix of $B$. Then similarly, we can construct a $2\times 2$ matrix $\tilde L_2$ subject to $$ \tilde L_2 \tilde B= \begin{pmatrix} \times & \times \\ 0 & \times \end{pmatrix} . $$ Here $$ \tilde L_2= \begin{pmatrix} 1 & 0 \\ -\frac{\tilde b_{21}}{\tilde b_{11}} & 1 \end{pmatrix}. $$

Consequently $$ L_2B= L_2 L_1 A= \begin{pmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0 & \times \end{pmatrix} $$ where $$ L_2= \text{diag}(1, L_2)= \begin{pmatrix} 1 & & \\ & 1 & \\ & -\frac{\tilde b_{21}}{\tilde b_{11}} & 1 \end{pmatrix}. $$

Consequently $A= L_1^{-1} L_2^{-1} U$ where $U$ is an upper triangular matrix. This is the LU decomposition: indeed, $L_1^{-1} L_2^{-1}$ is a lower triangular matrix. To see this, note that $L_2^{-1}$ is just changing the sign of $(3,2)$ entry of $L_2$, and $L_1^{-1}$ is changing the sign of $(2,1), (3,1)$ entries of $L_1$. Therefore $L_i$ and $L_i^{-1}$ have the same shape. And $L_1^{-1} L_2^{-1}$ is lower triangular. Their product can be computed by considering corresponding elementary row transformation.