Math 152: Intro to Mathematical Software

2017-02-06

Kiran Kedlaya; University of California, San Diego

adapted from lectures by William Stein, University of Washington

Lecture 12: Linear algebra (part 3)

As usual, HW due Tuesday at 8pm. Peer evaluations due Thursday at 8pm.

Multimodular echelon form: a peek behind the scenes

How does a computer algebra system compute the determinant of a large matrix?

Definitely not by expansion in minors, nor by using the explicit formula as a sum over permutations. The most reasonable generic method is Gaussian elimination, which for a matrix of size n requires $O(n^3)$ arithmetic operations.

For very large matrices (say n > 50000), there are some more advanced methods that do a bit better, e.g., Strassen's algorithm: https://en.wikipedia.org/wiki/Strassen_algorithm. We won't be encountering such methods here.

However, there is a big difference between doing linear algebra with floating-point real numbers (as is typical in MATLAB) versus doing it with exact numbers like integers or rational numbers. Let me illustrate this now.

First, though, an example to demonstrate what Sage is capable of.

The idea behind this calculation is multimodular echelon form. This uses some ideas from elementary number theory in a way you might not have anticipated.

Imagine with me what it would be like to do Gaussian elimination to put the following matrix $$ into reduced row Echelon form:

Then clear the second column:

And it gets bad...

There is a better way, which seems way more complicated but:

• avoids "coefficient explosion"

• enables many other nice optimizations: parallelization, working mod p, asymptotically fast matrix multiplication, etc.

The basic idea:

1. Choose some random prime $p$.

2. Reduce our matrix modulo a prime number $p$ to get $A\pmod{p}$.

3. Compute the reduced row echelon form of $A \pmod{p}$.

4. (*) Lift the result back to the rational numbers to get the right answer.

This is potentially better, since the numbers all stay small, instead of getting massive/painful.

File: /projects/sage/sage-7.5/src/sage/matrix/matrix0.pyx
Signature : B.add_multiple_of_row(self, i, j, s, start_col=0)
Docstring :
Add s times row j to row i.

EXAMPLES: We add -3 times the first row to the second row of an
integer matrix, remembering to start numbering rows at zero:

   sage: a = matrix(ZZ,2,3,range(6)); a
   [0 1 2]
   [3 4 5]
   sage: a.add_multiple_of_row(1,0,-3)
   sage: a
   [ 0  1  2]
   [ 3  1 -1]

To add a rational multiple, we first need to change the base ring:

   sage: a = a.change_ring(QQ)
   sage: a.add_multiple_of_row(1,0,1/3)
   sage: a
   [   0    1    2]
   [   3  4/3 -1/3]

If not, we get an error message:

   sage: a.add_multiple_of_row(1,0,i)
   Traceback (most recent call last):
   ...
   TypeError: Multiplying row by Symbolic Ring element cannot be done over Rational Field, use change_ring or with_added_multiple_of_row instead.

Notice -- the numbers aren't exploding. (Please imagine doing all this with a 200x201 matrix, say!)

But... how do you get from B.rref() to A.rref()???

Obviously, there isn't enough information: $$-\frac{3231}{1600} \mapsto 236 \pmod{389}$$ but 236 mod 389 doesn't tell you $-\frac{3231}{1600}$...

There is an algorithm called rational reconstruction, which is based on Euclidean GCD:

rational_reconstruction: compute x/y, where x/y is a rational number
in lowest terms such that the reduction of x/y modulo m is equal to
a and the absolute values of x and y are both <= sqrt{m/2}. If
such x/y exists, that pair is unique and this function returns it.
If no such pair exists, this function raises ZeroDivisionError.

Still doesn't work. But if we replace 389 by a bigger prime, it works:

Multimodular!

However, if the matrix is large, the entries are large, and we have to work modulo an enormous prime with thousands of digits. NO.

• Instead, work modulo several small primes (possibly at once, in parallel). Small primes are easier and faster -- they fit in a machine word.

• Use the Chinese Remainder Theorem to combine the rref modulo a bunch of small primes $p_i$ to get a matrix with entries modulo $M = p_1 p_2 \cdots p_n$.

• Use rational reconstruction on that (it works modulo $M$ for any integer $M$).

• And... realize that echelon modulo $p$ can be done via several matrix multiplications! (Permute the matrix and work with blocks...)

• And... that matrix multiplication is way faster that $O(n^3)$ for big $n$... https://en.wikipedia.org/wiki/Strassen_algorithm again.