You should have learned the algorithm of Gaussian elimination. If you took a numerical analysis class, you should also have learned that it can be usefully viewed as a matrix factorization. (To be perfectly honest, that never penetrated my own mind until I started teaching numerical methods classes---it's one thing to watch it on the blackboard, another to program it, and still another to try to teach it.) This fact is worth knowing, so here is an example.

with(linalg):

A := randmatrix(4,4);

In Gaussian elimination, we want to use the pivot **85** to eliminate
everything below it. To do it as an elementary matrix operation,
we create a 4 by 4 upper triangular elementary matrix (i.e. a matrix
that has only one column different from the identity). We choose
that column to be exactly the elements of the first column of **A**,
divided by the pivot .

L1 := matrix(4,4,0):

for i to 4 do L1[i,1] := A[i,1]/A[1,1]; od:

for i to 4 do L1[i,i] := 1; od:

print(L1);

The inverse of this matrix is nearly the same matrix, with just the signs below the diagonals changed (this is true of all elementary matrices).

inverse(L1);

When we do Gaussian elimination we add times the
first row of **A** to the second row; likewise times
the first row of **A** to the third row; and also times
the first row to the last row. Inspection of the matrix multiplication
below shows that this is exactly what we are doing, in matrix
multiplication rather than row operations. This is an important
point (the key idea of the whole works) so you should work it
through by hand until you are quite convinced that this multiplication
really is doing all three row operations.

evalm(" &* A);

A1 := ":Now we do the second column, in a similar way. The pivot is now .

L2 := matrix(4,4,0):

for i from 2 to 4 do L2[i,2] := A1[i,2]/A1[2,2]; od:

L2[1,1] := 1;

for i from 3 to 4 do L2[i,i] := 1; od:

print(L2);

Again the inverse just changes the signs of the elementary matrix, because the inverse row operation is just adding the negative of the 2nd row instead of adding the row.

inverse(L2);

A2 := evalm(" &* A1);

We are now ready for the last one.

L3 := matrix(4,4,0):

for i to 4 do L3[i,i] := 1; od:

L3[4,3] := A2[4,3]/A2[3,3];

print(L3);

evalm( L3^(-1) &* A2);

Note that that is upper triangular; we are done.

U := ":

evalm( L1 &* L2 &* L3 &* U);

evalm("-A);

We have thus factorized . But,

evalm( L1 &* L2 &* L3 );

so the product is itself unit lower triangular.

L := ":

evalm( A - L&*U);

The trick with row exchanges is also important, but completely
unconvincing unless you work it through yourself. You should prove
to yourself that multiplying by a permutation matrix (i.e. the identity
matrix with two rows exchanged) simply exchanges rows. Thus you
can encode the action of row exchanges also as matrix multiplication;
by a rather nifty bit of magic these row exchanges can all be placed
at the front. Work through the details yourself. The end result is
that any nonsingular matrix can be written as **PA = LU** where **P** is some permutation of the identity matrix, **L** is unit lower triangular, and **U** is upper triangular. We will return to this when we study the Row Echelon Decomposition.

Wed Jan 31 11:33:59 EST 1996