A \(5\times 7\) matrix, of rank \(5\).

We use row operations, strictly of the third type, to “zero out” entries of columns below the diagonal.

Column one.

Column two.

Column three.

That's it. Done. Everything below the diagonal of this matrix is zero, so the matrix is upper-triangular. Notice that the result is not in reduced row-echelon form and may not even have linearly independent rows.

Let's collect the pieces.

First, the upper-triangular matrix \(U\).

The elementary matrices are square and invertible, so have inverses individually, and so does their product. We will move the product to the other side of the equation by taking an inverse. As a product of elementary matrices of a very specific type, which are all lower-triangular with 1's on the diagonal, the inverse will also be lower-triangular with 1's on the diagonal.

Check our work, we should have \(LU = A\).

Here is the Sage routine.

Read Sage online help about the exact LU method, especially pivoting and a storage efficient return value.

RDF is the “Real Double Field”, 64-bit imitations of real numbers (two 32-bit words). This gives 53 bits of precision in the significand (nee mantissa), 11 bits for the exponent, 1 bit for the sign. Also known as “double-precision floating point” numbers or IEEE 754. Do not confuse this with Sage's RealField which allows you to specify arbitraty precision, and which defaults to 53 bits as RR.

And a rectangular matrix.

To solve \(A\vec{x}=\vec{b}\), first solve \(L\vec{y}=\vec{b}\), by successively solving for the first entries of \(\vec{y}\) first. We will simulate this in Sage.

Now, solve for \(\vec{x}\) from \(U\vec{x}=\vec{y}\), since \[A\vec{x}=LU\vec{x}=L\vec{y}=\vec{b}\]

Notice that this would be “backsolving”, successively getting the entries of \(\vec{x}\) with the last entries first, and in this case by setting the “extra” variables to zero as a convenience.

And check we have a solution: