Let's make up a matrix

Gaussian elimination: towards the wholistic view A=LU through elimination matrices¶

Let's look more closely at the process of Gaussian elimination in matrix form, using the matrix from lecture 1.

Gaussian elimination produces the matrix $U$, which we can compute in Julia as in lecture 1:

Now, let's go through Gaussian elimination in matrix form, by expressing the elimination steps as matrix multiplications. In Gaussian elimination, we make linear combination of rows to cancel elements below the pivot, and we now know that this corresponds to multiplying on the left by some elimination matrix $E$.

The first step is to eliminate in the first column of $A$. The pivot is the 1 in the upper-left-hand corner. For this $A$, we need to:

1. Leave the first row alone.
2. Subtract the first row from the second row to get the new second row.
3. Subtract $3 \times {}$ first frow from the third row to get the new third row.

This corresponds to multiplying $A$ on the left by the matrix E1. As above (in the "row × matrix" picture), the three rows of E1 correspond exactly to the three row operations listed above:

We would love for you to see the above not as a mechanical observation but interpret it as reversing the row operationa above.

As desired, this introduced zeros below the diagonal in the first column. Now, we need to eliminate the 2 below the diagonal in the second column of E1*A. Our new pivot is $-2$ (in the second row), and we just add the second row of E1*A with the third row to make the new third row.

This corresponds to multiplying on the left by the matrix E2, which leaves the first two rows alone and makes the new third row by adding the second and third rows:

As expected, this is upper triangular, and in fact the same as the U matrix returned by the Julia lu function above:

Thus, we have arrived at the formula: $$ \underbrace{E_2 E_1}_E A = U $$ Notice that we multiplied $A$ by the elimination matrices from right to left in the order of the steps: it is $E_2 E_1 A$, not $E_1 E_2 A$. Because matrix multiplication is generally not commutative, $E_2 E_1$ and $E_1 E_2$ give different matrices:

Notice, furthermore, that the matrices $E_1$ and $E_2$ are both lower-triangular matrices with ones on the diagonal. This is a consequence of the structure of Gaussian elimination (assuming no row re-ordering): we always add the pivot row to rows below it, never above it.

The product of lower-triangular matrices is always lower-triangular too. (In homework, you will explore a similar property for upper-triangular matrices) In consequence, the product $E = E_2 E_1$ is lower-triangular, and Gaussian elimination can be viewed as yielding $EA=U$ where $E$ is lower triangular and $U$ is upper triangular.

Inverse elimination: LU factors¶

We can write $EA=U$ as $A= E^{-1} U$, where $E^{-1}$ is the inverse of the matrix $E$. We will have more to say about matrix inverses later in 18.06, but for now we just need to know that it is the matrix that reverses the steps of Gaussian elimination, taking us back from $U$ to $A$. Computing matrix inverses is laborious in general, but in this particular case it is easy. We just need to reverse the steps one by one starting with the last elimination step and working back to the first one.

Hence, we need to reverse (invert) $E_2$ first on $U$, and then reverse (invert) $E_1$: $A = E_1^{-1} E_2^{-1} U$. But reversing an individual elimination step like $E_2$ is easy: we just flip the signs below the diagonal, so that wherever we added the pivot row we subtract and vice-versa. That is: $$ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{pmatrix}^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \end{pmatrix} $$ (The last elimination step was adding the second row to the third row, so we reverse it by subtracting the second row from the third row of $U$.)

Julia can compute matrix inverses for us with the inv function. (It doesn't know the trick of flipping the sign, which only works for very special matrices, but it can compute it the "hard way" so quickly (for such a small matrix) that it doesn't matter.) Of course that gives the same result:

Similarly for $E_1$:

If we didn't make any mistakes, then $E_1^{-1} E_2^{-1} U$ should give $A$, and it does:

We use the letter "L" for the inverse elimination matrix $L = E^{-1} = E_1^{-1} E_2^{-1}$ Since the inverses of each elimination matrix were lower-triangular (with flipped signs), their product $L$ is also lower triangular:

(Think of a mechanical reason that products of lower-triangular are lower-triangular and a conceptual reason.)

As mentioned above, this is the same as the inverse of $E = E_2 E_1$:

The final result, therefore, is that Gaussian elimination (without row swaps) can be viewed as a factorization of the original matrix $A$ $$ A = LU $$ into a product of lower- and upper-triangular matrices. (Furthermore, although we didn't comment on this above, $L$ is always 1 along its diagonal.) This factorization is called the LU factorization of $A$. (It's why we used the lu function in Julia above.) When a computer performs Gaussian elimination, what it computes are the $L$ and $U$ factors.

What this accomplishes is to break a complicated matrix $A$ into much simpler pieces $L$ and $U$. It may not seem at first that $L$ and $U$ are that much simpler than $A$, but they are: lots of operations that are very difficult with $A$, like solving equations or computing the determinant, become easy once you known $L$ and $U$.

E1 and E2 have the negative multipliers below the diagonal

inv(E1) and inv(E2) have the multipliers below the diagonal.

L has the multipliers below the diagonal.

inv(L) does not generally have multipliers anywhere.