Gaussian Elimination = Matrix multiplications¶

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 row 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:

As desired, this introduced zeros below the diagonal in the first column.

It's fun and perhaps illuminating to try acting E1 on a matrix of symbolic variables, using the Symbolics.jl computer-algebra (symbolic-math) package for Julia.

We'll define a matrix S whose entries are variables instead of numbers, and then act E1 on it.

We could use boring variables like $s_{11}$ etcetera:

But it's hard to visually distinguish all of these subscripts, and nowadays on a computer we have a lot more symbols to choose from.

Let's instead make a symbolic matrix full of emoji:

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 that when we multiply $E_2 E_1$ to get $E$, the elimination steps get "mixed up together", and we get new numbers like -4 (above) that didn't appear as multipliers in the elimination steps. As discussed below, this makes $E$ a bit inconvenient as a way to work with Gaussian elimination, because it requires some pointless effort to compute $E$.)

Elimination without row swaps: U = (lower triangular) * A¶

Notice, furthermore, that the matrices $E_1$ and $E_2$ are both lower-triangular matrices. 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.

Caveat: there is a slight complication if we have to do row swap steps, since row swaps are not lower triangular. We will defer this case until later, since row swaps don't arise in hand calculations unless you are unlucky.

The other way around: A = LU¶

There are two problems with writing $U = EA$.

1. The matrix $E$ is annoying to compute — the different elimination steps are mixed up together, and we don't actually want to multiply those matrices together.

2. We want to think of this process as simplifying A. i.e. we want to replace $A$ with $A = \mbox{simpler stuff}$.

Both of these concerns are addressed by thinking about Gaussian elimination in reverse: How do we get A back from U? Again, this consists of row operations on $U$, and will result in

$$ A = LU $$

where $L$ is a lower-triangular matrix with 1's on the diagonal, and the multipliers from the elimination steps written below the diagonal.

See handwritten notes on:

1. How to get $L$: just a record of the elimination steps, no extra computation!
2. How to use $A = LU$ to solve $Ax=b$ for any right-hand side(s) we want.

Check:

The Identity Matrix¶

A very important special matrix is the identity matrix $I$. Here is a $ 5 \times 5$ identity matrix:

$$ I = \begin{pmatrix} 1&0&0&0&0 \\ 0&1&0&0&0 \\ 0&0&1&0&0 \\ 0&0&0&1&0 \\ 0&0&0&0&1 \end{pmatrix} = \begin{pmatrix} e_1 & e_2 & e_3 & e_4 & e_5 \end{pmatrix} $$

where the columns $e_1 \cdots e_5$ of $I$ are the unit vectors in each component.

The identity matrix, which can be constructed by I(5) in Julia, has the property that $Ix=x$ for any $x$, and hence $IA = A$ for any (here $5 \times 5$) matrix $A$:

Notational ambiguity¶

In principle, "$I$" is a little ambiguous in linear algebra, because there are $m \times m$ identity matrices of every size $m$, and $I$ could refer to any of them. In practice, you usually infer the size of $I$ from context: if you see an expression like $Ix$ or $IA$ or $A + I$, then $I$ refers to the identity matrix of the corresponding size.

There is actually a built-in constant called I in Julia that represents just that: an identity matrix whose size is inferred from context:

When there is an ambiguity, I will sometimes write $I_m$ for the $m \times m$ identity matrix.

For example, if $B$ is an $3\times2$ matrix, then

$$ I_3 B I_2 = B $$

where we need to multiply by a different-size identity matrix on the left and right.

Matrix Inverses¶

It is often conceptually convenient to talk about the inverse $A^{-1}$ of a matrix $A$, which exists for any non-singular square matrix. This is the matrix such that $x = A^{-1}b$ solves $Ax = b$ for any $b$. The inverse is conceptually convenient becuase it allows us to move matrices around in equations almost like numbers (except that matrices don't commute!).

Equivalently, the inverse matrix $A^{-1}$ is the matrix such that $A^{-1} A = A A^{-1} = I$.

Why does this correspond to solving $Ax=b$? Multiplying both sides on the left by $A^{-1}$ (multiplying on the right would make no sense: we can't multiply vector×matrix!), we get

$$ A^{-1}Ax=Ix = x = A^{-1}b $$

In Julia, you can just do inv(A) to get the inverse of a matrix:

(This is the identity matrix up to roundoff errors: computers only keep about 15 significant digits when they are doing most arithmetic.)

Inverses and products¶

Inverses have a special relationship to matrix products: $$ (AB)^{-1} = B^{-1} A^{-1} $$

The reason for this is that we must have $(AB)^{-1} AB = I$, and it is easy to see that $B^{-1} A^{-1}$ does the trick. Equivalently, $AB$ is the matrix that first multiplies by $B$ then by $A$; to invert this, we must reverse the steps: first multiply by the inverse of $A$ and then by the inverse of $B$.

Gaussian Elimination, Inverses, and LU Factorization¶

Recall from the last lecture that we were looking at the Gaussian-elimination steps to convert a matrix $A$ into upper-triangular form via a sequence of row operations:

$$ \underbrace{\begin{pmatrix} \boxed{1} & 3 & 1 \\ 1 & 1 & -1 \\ 3 & 11 & 6 \end{pmatrix}}_A\to \underbrace{\begin{pmatrix} \boxed{1} & 3 & 1 \\ 0 & \boxed{-2} & -2 \\ 0 & 2 & 3 \end{pmatrix}}_{E_1 A}\to \underbrace{\begin{pmatrix} \boxed{1} & 3 & 1 \\ 0 & \boxed{-2} & -2 \\ 0 & 0 & \boxed{1} \end{pmatrix}}_{E_2 E_1 A} $$

Since row operations correspond to multiplying by matrices on the *left*, constructed the "elimination matrices" $E_1$ and $E_2$ corresponding to each of these elimination steps (under the first and second pivots).

Thus, we 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. 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¶

However, in practice, it turns out to be more useful to write this as

$$A= E^{-1} U = E_1^{-1}E_2^{-1}U$$

where $E^{-1}$ is the inverse of the matrix $E$. We will have more to say about how to compute 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, just as for our general rule about inverses of products above, 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 call 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:

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

WARNING: For MOST triangular matrices, negation ≠ inverse¶

Above, we used the trick that we can invert a single elimination matrix $E$ just by flipping the signs below the diagonal (to reverse the steps). This doesn't work for a general lower (or upper) triangular matrix.

For example, even for the $L$ matrix computed above, it doesn't work:

Notice that the 3 in the lower-left corner of $L$ changes to a -4 in the lower-left of $L^{-1}$!

The reason it works for the "elementary" $E$ matrices representing elimination of a single column of $A$, but not for $L$ (or other lower-triangular matrices in general), is essentially that operations on different columns affect one another (the elimination steps don't commute).

In fact, there is no easy trick to invert a general lower/upper triangular matrix quickly. But you shouldn't have to! If you see an expression like $L^{-1} b$, you shouldn't compute $L^{-1}$ explicitly! Instead, solve $Lc = b$ by forward-substitution.

Similarly, if you see $U^{-1} c$ for an upper-triangular matrix, don't try to compute $U^{-1}$. Instead, solve $Ux=c$ by back-substitution.

Finding Matrix Inverses in General¶

How do we find $A^{-1}$? The key is the equation $A A^{-1} = I$, which looks just like $AX=B$ for the right-hand sides consisting of the columns of the identity matrix, i.e. the unit vectors. So, we just solve $Ax=e_i$ for $i=1,\ldots,5$, or equivalently do A \ I in Julia. Of course, Julia comes with a built-in function inv(A) for computing $A^{-1}$ as well:

The difference is just roundoff errors. We can also check approximate equality (ignoring roundoff differences) with ≈ (typed by \approx followed by a tab):

(Again, we get $I$ up to roundoff errors because the computer does arithmetic only to 15–16 significant digits.)

Normally, $AB \ne BA$ for two matrices $A$ and $B$. Why can we multiply $A$ by $A^{-1}$ on either the left or right and get the same answer $I$? It is fairly easy to see why:

$$ A A^{-1} = I \implies A A^{-1} A = IA = A = A (A^{-1} A) $$

Since $A (A^{-1} A) = A$, and $A$ is non-singular (so there is a unique solution to this system of equations), we must have $A^{-1} A = I$.

(Almost) Never Compute Inverses!¶

Matrix inverses are funny, however:

• Inverse matrices are very convenient in analytical manipulations, because they allow you to move matrices from one side to the other of equations easily.

• Inverse matrices are almost never computed in "serious" numerical calculations. Whenever you see $A^{-1} B$ (or $A^{-1} b$), when you go to implement it on a computer you should read $A^{-1} B$ as "solve $AX = B$ by some method." e.g. solve it by A \ B or by first computing the LU factorization of $A$ and then using it to solve $AX = B$.

One reason that you don't usually compute inverse matrices is that it is wasteful: once you have $A=LU$ (later we will generalize this to "$PA = LU$"), you can solve $AX=B$ directly without bothering to find $A^{-1}$, and computing $A^{-1}$ requires much more work if you only have to solve a few right-hand sides.

Another reason is that for many special matrices, there are ways to solve $AX=B$ much more quickly than you can find $A^{-1}$. For example, many large matrices in practice are sparse (mostly zero), and often for sparse matrices you can arrange for $L$ and $U$ to be sparse too. Sparse matrices are much more efficient to work with than general "dense" matrices because you don't have to multiply (or even store) the zeros. Even if $A$ is sparse, however, $A^{-1}$ is usually non-sparse, so you lose the special efficiency of sparsity if you compute the inverse matrix.