Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 16659

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

%md
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.

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(n3)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.

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!
-13466818614968645388082462498422891015858287714726503794914038833945861569879970893581274912807238067067562025415287158383559222384755019001390597752078698707012624901545748446973090664917103178828516044189989001058357767454396881652920277914269319561217079573570659118407585015540402467663013681892800682249262048905167741256172427519843205754072237734411450470997774246258609866816146985342497454615518377036304588373323896014014153096848684377338304789539220923360962433254751558698140098724370136565339961926983910452181415886121505372040127268243507290846872169186475717059074645991467437859199330798938529547771092330020701868445905292646258531607866600842701587778185479388594835136523597069495593786781542073910866710039961581522597540853906722624425777046752278726134643069795751472843281524652048263868977198422465443386463737170997961988241558325856098868881117286275698821414313043103155336005788593749834207684611189940389383088350198647194034148864729303340278745576761416302981482824843383373238970004540149993657984004685162519644032629699572231001020702562449943359701653788476269539085299139066850745562869137873023540296491189356451161184274022043422995011596989450637487043437899594602518997763901460637176662653607211024812624815151864798193588316950454139880552190319283929815428837956660130798702278897904798459454372741620060528913131832383053253381573068718059313121337056771116385810113530264101468068995995202495314749310543288967277221599107604456092250888496767499725158744437434425576614496008169674006717193123667215583466734391507627512283239296555117935079148978450406125656027756184451621783133121524773894183624658586629927617783715876317946019038678980403124019395962311780888290679041009997618540149072078353486691005953834045515800321231198540111091957873658942367757501944543442129821705696468018641544969784727847193758878349397089010123210025023637171438624785267726671716028471005028755778859426805746461108025364807260717399679397901566492123408820734586154253232736561172758746875406739220710852358873 CPU time: 11.70 s, Wall time: 13.90 s

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.

set_random_seed(42)
A = random_matrix(QQ, 4, 5, num_bound=10, den_bound=1)
show(A)
(−9−7−9−2−4−36−10−4−3−8−82−10−76−233−7)\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]
(0, 25/3, -7, -10/3, -5/3) 
A.add_multiple_of_row(1, 0, -3/9)
show(A)
(−9−7−9−2−40253−7−103−53−8−82−10−76−233−7)\displaystyle \left(\begin{array}{rrrrr} -9 & -7 & -9 & -2 & -4 \\ 0 & \frac{25}{3} & -7 & -\frac{10}{3} & -\frac{5}{3} \\ -8 & -8 & 2 & -10 & -7 \\ 6 & -2 & 3 & 3 & -7 \end{array}\right)
A.add_multiple_of_row(2, 0, -8/9)
show(A)
(−9−7−9−2−40253−7−103−530−16910−749−3196−233−7)\displaystyle \left(\begin{array}{rrrrr} -9 & -7 & -9 & -2 & -4 \\ 0 & \frac{25}{3} & -7 & -\frac{10}{3} & -\frac{5}{3} \\ 0 & -\frac{16}{9} & 10 & -\frac{74}{9} & -\frac{31}{9} \\ 6 & -2 & 3 & 3 & -7 \end{array}\right)
A.add_multiple_of_row(3, 0, 6/9)
show(A)
(−9−7−9−2−40253−7−103−530−16910−749−3190−203−353−293)\displaystyle \left(\begin{array}{rrrrr} -9 & -7 & -9 & -2 & -4 \\ 0 & \frac{25}{3} & -7 & -\frac{10}{3} & -\frac{5}{3} \\ 0 & -\frac{16}{9} & 10 & -\frac{74}{9} & -\frac{31}{9} \\ 0 & -\frac{20}{3} & -3 & \frac{5}{3} & -\frac{29}{3} \end{array}\right)
A.rescale_row(0, -1/9)
show(A)
(179129490253−7−103−530−16910−749−3190−203−353−293)\displaystyle \left(\begin{array}{rrrrr} 1 & \frac{7}{9} & 1 & \frac{2}{9} & \frac{4}{9} \\ 0 & \frac{25}{3} & -7 & -\frac{10}{3} & -\frac{5}{3} \\ 0 & -\frac{16}{9} & 10 & -\frac{74}{9} & -\frac{31}{9} \\ 0 & -\frac{20}{3} & -3 & \frac{5}{3} & -\frac{29}{3} \end{array}\right)

Then clear the second column:

A.rescale_row(1, 3/25)
show(A)
(1791294901−2125−25−150−16910−749−3190−203−353−293)\displaystyle \left(\begin{array}{rrrrr} 1 & \frac{7}{9} & 1 & \frac{2}{9} & \frac{4}{9} \\ 0 & 1 & -\frac{21}{25} & -\frac{2}{5} & -\frac{1}{5} \\ 0 & -\frac{16}{9} & 10 & -\frac{74}{9} & -\frac{31}{9} \\ 0 & -\frac{20}{3} & -3 & \frac{5}{3} & -\frac{29}{3} \end{array}\right)
A.add_multiple_of_row(0, 1, -7/9)
show(A)
(10124758153501−2125−25−150−16910−749−3190−203−353−293)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & \frac{124}{75} & \frac{8}{15} & \frac{3}{5} \\ 0 & 1 & -\frac{21}{25} & -\frac{2}{5} & -\frac{1}{5} \\ 0 & -\frac{16}{9} & 10 & -\frac{74}{9} & -\frac{31}{9} \\ 0 & -\frac{20}{3} & -3 & \frac{5}{3} & -\frac{29}{3} \end{array}\right)
A.add_multiple_of_row(2, 1, 16/9)
show(A)
(10124758153501−2125−25−150063875−13415−1950−203−353−293)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & \frac{124}{75} & \frac{8}{15} & \frac{3}{5} \\ 0 & 1 & -\frac{21}{25} & -\frac{2}{5} & -\frac{1}{5} \\ 0 & 0 & \frac{638}{75} & -\frac{134}{15} & -\frac{19}{5} \\ 0 & -\frac{20}{3} & -3 & \frac{5}{3} & -\frac{29}{3} \end{array}\right)
A.add_multiple_of_row(3, 1, 20/3)
show(A)
(10124758153501−2125−25−150063875−13415−19500−435−1−11)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & \frac{124}{75} & \frac{8}{15} & \frac{3}{5} \\ 0 & 1 & -\frac{21}{25} & -\frac{2}{5} & -\frac{1}{5} \\ 0 & 0 & \frac{638}{75} & -\frac{134}{15} & -\frac{19}{5} \\ 0 & 0 & -\frac{43}{5} & -1 & -11 \end{array}\right)

And it gets bad...



show(A.rref())
(1000−32311600010084596400001014171280000194696400)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & -\frac{3231}{1600} \\ 0 & 1 & 0 & 0 & \frac{8459}{6400} \\ 0 & 0 & 1 & 0 & \frac{1417}{1280} \\ 0 & 0 & 0 & 1 & \frac{9469}{6400} \end{array}\right)

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

  2. Reduce our matrix modulo a prime number pp to get A(modp)A\pmod{p}.

  3. Compute the reduced row echelon form of A(modp)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)
(−9−7−9−2−4−36−10−4−3−8−82−10−76−233−7)\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)
show(B)
B.parent()
(38038238038738538663793853863813812379382638733382)\displaystyle \left(\begin{array}{rrrrr} 380 & 382 & 380 & 387 & 385 \\ 386 & 6 & 379 & 385 & 386 \\ 381 & 381 & 2 & 379 & 382 \\ 6 & 387 & 3 & 3 & 382 \end{array}\right)
Full MatrixSpace of 4 by 5 dense matrices over Ring of integers modulo 389 
B.rescale_row(0, 1/380)
show(B)
(144134630338663793853863813812379382638733382)\displaystyle \left(\begin{array}{rrrrr} 1 & 44 & 1 & 346 & 303 \\ 386 & 6 & 379 & 385 & 386 \\ 381 & 381 & 2 & 379 & 382 \\ 6 & 387 & 3 & 3 & 382 \end{array}\right)
B.add_multiple_of_row?

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.

B.add_multiple_of_row(1, 0, -386)
show(B)
(144134630301383822561283813812379382638733382)\displaystyle \left(\begin{array}{rrrrr} 1 & 44 & 1 & 346 & 303 \\ 0 & 138 & 382 & 256 & 128 \\ 381 & 381 & 2 & 379 & 382 \\ 6 & 387 & 3 & 3 & 382 \end{array}\right)
B.add_multiple_of_row(2, 0, -381)
B.add_multiple_of_row(3, 0, -6)
show(B)
(1441346303013838225612803441035830123386261120)\displaystyle \left(\begin{array}{rrrrr} 1 & 44 & 1 & 346 & 303 \\ 0 & 138 & 382 & 256 & 128 \\ 0 & 344 & 10 & 35 & 83 \\ 0 & 123 & 386 & 261 & 120 \end{array}\right)
B.rescale_row(1, 1/138)
show(B)
(14413463030121723331103441035830123386261120)\displaystyle \left(\begin{array}{rrrrr} 1 & 44 & 1 & 346 & 303 \\ 0 & 1 & 217 & 233 & 311 \\ 0 & 344 & 10 & 35 & 83 \\ 0 & 123 & 386 & 261 & 120 \end{array}\right)
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)
(10178208234012172333110050177400147388378)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 178 & 208 & 234 \\ 0 & 1 & 217 & 233 & 311 \\ 0 & 0 & 50 & 17 & 74 \\ 0 & 0 & 147 & 388 & 378 \end{array}\right)

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

show(A.rref(), B.rref())
(-3231/1600) % 389
(1000−32311600010084596400001014171280000194696400)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & -\frac{3231}{1600} \\ 0 & 1 & 0 & 0 & \frac{8459}{6400} \\ 0 & 0 & 1 & 0 & \frac{1417}{1280} \\ 0 & 0 & 0 & 1 & \frac{9469}{6400} \end{array}\right) (1000236010012100103740001140)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & 236 \\ 0 & 1 & 0 & 0 & 121 \\ 0 & 0 & 1 & 0 & 374 \\ 0 & 0 & 0 & 1 & 140 \end{array}\right)
236

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

show(A.rref())

Obviously, there isn't enough information: −32311600↦236(mod389) -\frac{3231}{1600} \mapsto 236 \pmod{389} but 236 mod 389 doesn't tell you −32311600-\frac{3231}{1600}...

(-3231/1600) % 389
236 
Mod(236, 389).lift()
236

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()
13/5

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)
print(p)
C = A.change_ring(ZZ).mod(p)
show(C)
(−9−7−9−2−4−36−10−4−3−8−82−10−76−233−7)\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)
10000000019
(10000000010100000000121000000001010000000017100000000151000000001661000000000910000000015100000000161000000001110000000011210000000009100000000126100000000173310000000012)\displaystyle \left(\begin{array}{rrrrr} 10000000010 & 10000000012 & 10000000010 & 10000000017 & 10000000015 \\ 10000000016 & 6 & 10000000009 & 10000000015 & 10000000016 \\ 10000000011 & 10000000011 & 2 & 10000000009 & 10000000012 \\ 6 & 10000000017 & 3 & 3 & 10000000012 \end{array}\right)
C = C.rref()
show(C)
(10008431250014010019359375050010310156250700011326562504)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & 8431250014 \\ 0 & 1 & 0 & 0 & 1935937505 \\ 0 & 0 & 1 & 0 & 3101562507 \\ 0 & 0 & 0 & 1 & 1326562504 \end{array}\right)
show(A.rref())
(1000−32311600010084596400001014171280000194696400)\displaystyle \left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & -\frac{3231}{1600} \\ 0 & 1 & 0 & 0 & \frac{8459}{6400} \\ 0 & 0 & 1 & 0 & \frac{1417}{1280} \\ 0 & 0 & 0 & 1 & \frac{9469}{6400} \end{array}\right)
for i in range(4):
    print C[i,4].rational_reconstruction()
-3231/1600 8459/6400 1417/1280 9469/6400 
%md
It worked!

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

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 pip_i to get a matrix with entries modulo M=p1p2⋯pnM = p_1 p_2 \cdots p_n.

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

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

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

A = random_matrix(QQ, 500, 501, num_bound=10, den_bound=1)
%time R = A.rref()
CPU time: 1.97 s, Wall time: 2.69 s 
show(R[0,500])