2016-05-06.sagews

Math 480: Open Source Mathematical Software

2016-05-06

William Stein

Lectures 18: Linear algebra (part 3 of 3)

Today:

1. Reminder: homework and peer grading is due at 6pm tonight.

2. Start screencast (crossing fingers that projector works this time... if not, maybe I'll switch to this windows machine in front of me)

3. (20 minutes) How multimodular echelon form works

4. (5 minutes) Sync stress test

Multimodular echelon form

I told you before about one way Sage computes the determinant of a large integer (or rational) matrix $A$ by solving $Ax = v$ for a random $v$, and considering the denominator of $x$. Cramer's rule, turned on its head, is a powerful approach to computing determinants.

A = random_matrix(ZZ, 1000, x=-10, y=10);  # random 1000x1000 matrix over ZZ with entries in [-10..10]
%time A.det()   # should take around 10 seconds!

1172096739030457809484099809502450364959655827024062189836728991855081750183075132423763529970516305158219215057496322717866013897829659141110783185872685146901793137119436823833587443678879746455318765807744265187263154690090912038172062035136970301326918259969523822266327573151362777552237306308797036666750581394159902980823251030339946051932903064869487968076712315657784156871364172939435485130388152865428844389843511949932516070314520579871647668754587138646049990790944043279187533812936232782437207620334897651238112723330232144256829680385352397461431002563952039373816367581209085787553298612265130712275732953594321399282753304684043698154992748832454622952999709446991333036448882252727451630334303992724914733432462456658827221777323719282099557870138304304988441229342517993913279084042288736357404682621321427433296150940855089500568452577202145662206330167887388000714559897950468501205879183830015061573922621334469665969239280864858311331656691749763696838003264443026984447674018359318373150159588374658574776076655195737140789265023144502459259251270936405987220136660467081830634242132702460652309450554412906189355111764600003629569436933524145361731811702537132699432453580329282298390762182139697975416273206031858316868288104163366061249710206875294489545101546565632384546478810644467711780122714810077843268772030853876612194881047036348031533609925135367920346883455973883422041580493464604070708394698606687924376110596952755819604668331774129375293056327775787369021947337179627051633773257654825785529152854917529806095393388345690956383164132660331519253121438759756252569186450388803398756314613881053811357567704004264554377624338021382842252309621822581002556038003999425916582718274582715227463057729774319255377260680548549874205974436247058016750237332042830216729807727116667393612034098522994896486361628154001240939240355504132228451453649410052479719873962794867885335572831342826182058145954755196377759051347407675833810406211199958243488267437986515514406856389658577363899327317782477424169573356 CPU time: 8.87 s, Wall time: 9.01 s

Here's another central idea in exact linear algebra: multimodular echelon form.

set_random_seed(42)
A = random_matrix(QQ, 4, 5, num_bound=10, den_bound=1)
show(A)

$\displaystyle \left(\begin{array}{rrrrr} -9 & -7 & -9 & -2 & -4 \\ -3 & 6 & -10 & -4 & -3 \\ -8 & -8 & 2 & -10 & -7 \\ 6 & -2 & 3 & 3 & -7 \end{array}\right)$

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

A[0]*(-3/9) + A[1]

A.add_multiple_of_row(1, 0, -3/9)
show(A)

A.add_multiple_of_row(2, 0, -8/9)
show(A)

A.add_multiple_of_row(3, 0, 6/9)
show(A)

A.rescale_row(0, -1/9)
show(A)

Then clear the second column:

A.rescale_row(1, 3/25)
show(A)

A.add_multiple_of_row(0, 1, -7/9)
show(A)

A.add_multiple_of_row(2, 1, 16/9)
show(A)

A.add_multiple_of_row(3, 1, 20/3)
show(A)

And it gets bad...

show(A.rref())

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.

set_random_seed(42)
A = random_matrix(QQ, 4, 5, num_bound=10, den_bound=1)
show(A)
B = A.change_ring(ZZ).mod(389)
show(B)

B.rescale_row(0, 1/380)
show(B)

B.add_multiple_of_row(1, 0, -386)
show(B)

B.add_multiple_of_row(2, 0, -381)
B.add_multiple_of_row(3, 0, -6)
show(B)

B.rescale_row(1, 1/138)
show(B)

B.add_multiple_of_row(0,1, -44)
B.add_multiple_of_row(2,1,-344)
B.add_multiple_of_row(3,1,-123)
show(B)

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

show(B.rref())

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

show(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}$...

(-3231/1600) % 389

Mod(236, 389).lift()

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.

Mod(236, 389).rational_reconstruction()

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

set_random_seed(42)
A = random_matrix(QQ, 4, 5, num_bound=10, den_bound=1)
show(A)
p = next_prime(10^10)
C = A.change_ring(ZZ).mod(p)
show(C)

C = C.rref()
show(C)

show(A.rref())

for i in range(4):
    print C[i,4].rational_reconstruction()

%md
It worked!

- how do we know for sure?: see Chapter 7 of http://wstein.org/books/modform/stein-modform.pdf

Multimodular!

But this sucks -- 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 an int or long.

• 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)$... https://en.wikipedia.org/wiki/Strassen_algorithm

A = random_matrix(QQ, 500, 501, num_bound=10, den_bound=1)
%time R = A.rref()

show(R[0,500])

Sync stress test

Everybody do something here and see what happens (tmp/2016-05-06-stress.sagews in the shared project):

https://cloud.sagemath.com/projects/b303c6cf-3159-479b-8aa5-c8d2abaed2bd/files/tmp/2016-05-06-stress.sagews