Since we're interested in the inverses of elimination matrices, we should think about how to consider the inverses of products of matrices.
Let $A$ and $B$ be matrices of the same size. If $A$ and $B$ are invertible, then $AB$ is also invertible. Furthermore, the inverse of $AB$, $(AB)^{-1}$ is $B^{-1} A^{-1}$.
This is important—because the order in which the matrices are multiplied matters, the order in which their inverses are multiplied also matters—we must reverse the order.
If we view matrices as a function that transforms vectors, then we can consider the transformation on a given $\mathbf x$. The product $AB\mathbf x$ can be thought of as a transformation that first creates $B \mathbf x = \mathbf y$ and then computes $A \mathbf y = \mathbf z$. \[\mathbf x \xrightarrow{B} \mathbf y \xrightarrow{A} \mathbf z\]
The procedure to reverse this process becomes clear if we have the inverses $A^{-1}$ and $B^{-1}$: we first need to compute $A^{-1} \mathbf z = \mathbf y$ and then take this vector to get $B^{-1} \mathbf y = \mathbf x$: \[\mathbf z \xrightarrow{A^{-1}} \mathbf y \xrightarrow{B^{-1}} \mathbf x\] This is exactly $(AB)^{-1} \mathbf z = B^{-1} A^{-1} \mathbf z$.
Recall that $E_{21} = \begin{bmatrix}1 & 0 & 0 \\ -4 & 1 & 0 \\0 & 0 & 1 \end{bmatrix}$ subtracts 4 of the first row from the second row and the inverse operation adds 4 of the first row back to the second row: $E_{21}^{-1} = \begin{bmatrix}1 & 0 & 0 \\ 4 & 1 & 0 \\0 & 0 & 1 \end{bmatrix}$.
Now, since we handled the second row, suppose we want to do elimination on the third row by using the second row. We use a matrix like $E_{32} = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\0 & -2 & 1 \end{bmatrix}$. As above, the inverse of this matrix is $E_{32}^{-1} = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\0 & 2 & 1 \end{bmatrix}$.
So suppose we accomplish $E_{32} E_{21} A = U$. What is $E_{32} E_{21}$? We can compute this: \[E_{32} E_{21} = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\0 & -2 & 1 \end{bmatrix}\begin{bmatrix}1 & 0 & 0 \\ -4 & 1 & 0 \\0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -4 & 1 & 0 \\ 8 & -2 & 1 \end{bmatrix}.\]
Now, we went through the trouble of thinking about inverses, so let's consider them. We have \[(E_{32} E_{21})^{-1} = E_{21}^{-1} E_{32}^{-1} = \begin{bmatrix}1 & 0 & 0 \\ 4 & 1 & 0 \\0 & 0 & 1 \end{bmatrix}\begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\0 & 2 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 0 & 2 & 1 \end{bmatrix}.\]
There's something in the "niceness" of the matrices that we get. The inverse of $E$ looks "nicer" or "cleaner" than $E$—the multiples go in the right spot, no fuss, no mess. But in $E$, we have the negatives, yes, but for some reason there's an extra entry. What's going on?
This leads to the following observation. We have been using elimination on a matrix $A$ to transform it into an upper triangular matrix $U$. Elimination can be expressed as a series of elimination matrices. Because elimination changes the bottom rows of a matrix based on the top rows, our elimination matrices are lower triangular matrices, with $1$s along the diagonal.
So if we take all of our elimination matrices, we can multiply them together to get a single matrix $E$. We then have the relationship $E_k E_{k-1} \cdots E_2 E_1 A = EA = U$, where $E_i$ is a single elimination matrix. This is a great result, but as we saw in the example above, the entries are not quite so satisfying. If $\ell_{ij}$ is the corresponding "multiple" from elimination matrix $E_{ij}$, then for a $3 \times 3$ matrix, we have \[E = \begin{bmatrix} 1 & 0 & 0 \\ -\ell_{21} & 1 & 0 \\ \ell_{32} \ell_{21} - \ell_{31} & -\ell_{32} & 1 \end{bmatrix}.\] Note that in this notation, if $E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$, then $\ell_{21} = 2$, not $-2$.
However, $E$ is invertible and to get $E^{-1}$, all we need to do is flip the signs for our elimination matrices and multiply them in reverse order. When we do, we get a very interesting looking matrix: \[E^{-1} = \begin{bmatrix} 1 & 0 & 0 \\ \ell_{21} & 1 & 0 \\ \ell_{31} & \ell_{32} & 1 \end{bmatrix}.\]
What's more is notice that we get \[E^{-1} E A = A = E^{-1} U.\]
Notice that like $E$, $E^{-1}$ is also a lower triangular matrix. Elimination matrices are lower triangular and their inverses are lower triangular and their products are lower triangular. Since $E^{-1}$ is lower triangular, we let $L = E^{-1}$. Then we get the factorization \[A = LU.\]
This is the LU-factorization.
Though $L$ is formally defined as the inverse of $E$, we can informally define it by relating its entries to the corresponding elimination matrices. The items that go in each corresponding entry of the elimination matrices are the additive inverses of entries from $L$. Hence, if $\ell_{21}$ is the $(2,1)$st entry of $L$, then $E_{21}$ holds $-\ell_{21}$ as its $(2,1)$st entry.
There is an intriguing question of why we get such a nice result when we consider the inverse of $E$. Recall that elimination proceeds from the top row and moves progressively down. So which rows are we taking linear combinations of when we perform elimination on, say, the third row?
The answer is that we are taking combinations of rows of $E_{21}A$ (not $A$!). But these rows are actually rows of $U$—they don't get changed after we move on to the third row (assuming we are working systematically). We have \[\text{Row 3 of $U$} = \text{Row 3 of A} - \ell_{31} \times \text{Row 1 of $U$} - \ell_{32} \times \text{Row 2 of $U$}\] Notice that we can take this and describe the corresponding row of $A$ by rearranging into \[\text{Row 3 of $A$} = \ell_{31} \times \text{Row 1 of $U$} + \ell_{32} \times \text{Row 2 of $U$} + \text{Row 3 of $U$}.\] This is what the third row of $L$ is doing in the product $A = LU$: it is describing how to combine rows of $U$ to create rows of $A$. This same argument applies for larger matrices: when we're performing elimination on row $k$, we are taking linear combinations of rows 1 through $k$ of $U$.
But there's more—recall that we can view matrix multiplication as a sum of rank-one matrix products. That is, for $A = LU$, we take the $i$th column of $L$ and multiply it by the $i$th row of $U$. What does this give us?
If we consider the first row of $U$, this is exactly the first row of $A$—we never do anything to this row. But what is the first column of $L$? It describes how we eliminated the entries underneath. So multiplying the first column of $L$ and the first row of $U$ tells how to reconstruct the first column of $A$! \[\begin{bmatrix} 1 \\ \ell_{21} \\ \ell_{31} \\ \ell_{41} \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & \ell_{21} a_{12} & \ell_{21} a_{13} & \ell_{21} a_{14} \\ a_{31} & \ell_{31} a_{12} & \ell_{31} a_{13} & \ell_{31} a_{14} \\ a_{41} & \ell_{41} a_{12} & \ell_{41} a_{13} & \ell_{41} a_{14} \end{bmatrix}. \]
Then we repeat this process for the second, third, and so on rows of $U$ and columns of $L$. For instance, for $i = 2$, we have \[\begin{bmatrix} 0 \\ 1 \\ \ell_{32} \\ \ell_{42} \end{bmatrix} \begin{bmatrix} 0 & u_{22} & u_{23} & u_{24} \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & u_{22} & \ell_{22} u_{23} & \ell_{22} u_{24} \\ 0 & \ell_{32} u_{22} & \ell_{32} u_{23} & \ell_{32} u_{24} \\ 0 & \ell_{42} u_{22} & \ell_{42} u_{23} & \ell_{42} u_{24} \end{bmatrix}. \] But recall that, for example $u_{22} = a_{22} - \ell_{21}a{12}$. So if we add the above two matrices together, we notice that the $2,2$th entry is $u_{22} + \ell_{21}a_{21} = a_{22}$. Similarly, we have \[a_{32} = \ell_{32}u_{22} + \ell_{21} a_{12},\] or in other words, the original $3,2$th entry can be reconstructed by remembering how many times we removed multiples of $a_{12}$ and $u_{22}$ from it.
So the $i$th column of $L$ describes how the entries underneath were eliminated based on the previous rows of $U$. Multiplying these together and adding up all the matrices we get puts things back in place.
Recall our original goal of solving $A \mathbf x = \mathbf b$. Let $A = \begin{bmatrix} 2 & -1 & 5 \\ 4 & 1 & 19 \\ 2 & 8 & 34\end{bmatrix}$. The LU factorization fo $A$ is $U = \begin{bmatrix} 2 & -1 & 5 \\ 0 & 3 & 9 \\ 0 & 0 & 2 \end{bmatrix}$ and $L = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & 3 & 1 \end{bmatrix}$. We can use this to solve $A \mathbf x = \mathbf b = \begin{bmatrix}5 \\ 31 \\ 72 \end{bmatrix}$. To do so, we need to transform $A \mathbf x = \mathbf b$ into a system $U \mathbf x = \mathbf c$. But we know that $\mathbf c = E \mathbf b = L^{-1} \mathbf b = \begin{bmatrix}5 \\ 21 \\ 4 \end{bmatrix}$. Now, we can solve $U \mathbf x = \mathbf c$ to get $\mathbf x$ for $A \mathbf x = \mathbf b$.