Speed-up solving linear systems

[rbeezer]

Patch implements solving a linear system of equations in the most naive way possible, just augmenting the matrix and row-reducing.

For fields like ZZ, QQ, and integers mod p, this can be 3% to 10% faster. For fields that use echelon form to get rank, this can be twice as fast since we only row-reduce once, not twice. Matrices full of integers mod p can be a toss-up as the number of columns in the constant matrix is about 10 times greater than for the coefficient matrix.

Timings below and script that produced them is attached.

This has the old doctests, which pass with the new method (except two trivial failures). Old method is included as solve_left_old() for ease of testing timings. This is fully functional, but will need just a bit more work to be ready, so this is up for comments and suggestions at the moment.

Field: Integer Ring

Rows: 100
Cols (coeffs)     100
Cols (constants)  1
Old:
15 loops, best of 2: 30.659 ms per loop
New:
15 loops, best of 2: 26.411 ms per loop

Rows: 100
Cols (coeffs)     500
Cols (constants)  1
Old:
15 loops, best of 2: 1.3159 s per loop
New:
15 loops, best of 2: 1.3377 s per loop

Rows: 100
Cols (coeffs)     100
Cols (constants)  1000
Old:
15 loops, best of 2: 3.3781 s per loop
New:
15 loops, best of 2: 3.3435 s per loop

Rows: 100
Cols (coeffs)     300
Cols (constants)  200
Old:
15 loops, best of 2: 1.7518 s per loop
New:
15 loops, best of 2: 1.346 s per loop
**************************

Field: Rational Field

Rows: 100
Cols (coeffs)     100
Cols (constants)  1
Old:
15 loops, best of 2: 23.646 ms per loop
New:
15 loops, best of 2: 18.902 ms per loop

Rows: 100
Cols (coeffs)     500
Cols (constants)  1
Old:
15 loops, best of 2: 1.0791 s per loop
New:
15 loops, best of 2: 1.0609 s per loop

Rows: 100
Cols (coeffs)     100
Cols (constants)  1000
Old:
15 loops, best of 2: 2.8108 s per loop
New:
15 loops, best of 2: 2.7389 s per loop

Rows: 100
Cols (coeffs)     300
Cols (constants)  200
Old:
15 loops, best of 2: 1.3419 s per loop
New:
15 loops, best of 2: 1.0658 s per loop
**************************

Field: Ring of integers modulo 6337

Rows: 150
Cols (coeffs)     150
Cols (constants)  1
Old:
15 loops, best of 2: 30.246 ms per loop
New:
15 loops, best of 2: 27.992 ms per loop

Rows: 150
Cols (coeffs)     750
Cols (constants)  1
Old:
15 loops, best of 2: 131.07 ms per loop
New:
15 loops, best of 2: 143 ms per loop

Rows: 150
Cols (coeffs)     150
Cols (constants)  1500
Old:
15 loops, best of 2: 585.01 ms per loop
New:
15 loops, best of 2: 636.57 ms per loop

Rows: 150
Cols (coeffs)     450
Cols (constants)  300
Old:
15 loops, best of 2: 479.46 ms per loop
New:
15 loops, best of 2: 204.56 ms per loop
**************************

Field: Finite Field in a of size 5^4

Rows: 50
Cols (coeffs)     50
Cols (constants)  1
Old:
15 loops, best of 2: 35.873 ms per loop
New:
15 loops, best of 2: 17.424 ms per loop

Rows: 50
Cols (coeffs)     250
Cols (constants)  1
Old:
15 loops, best of 2: 179.94 ms per loop
New:
15 loops, best of 2: 139.05 ms per loop

Rows: 50
Cols (coeffs)     50
Cols (constants)  500
Old:
15 loops, best of 2: 346.74 ms per loop
New:
15 loops, best of 2: 339.3 ms per loop

Rows: 50
Cols (coeffs)     150
Cols (constants)  100
Old:
15 loops, best of 2: 360.61 ms per loop
New:
15 loops, best of 2: 139.35 ms per loop
**************************

[dimpase]

For QQ and other infinite exact fields, naive Gauss elimination is not a good strategy all together, as it can lead to an exponential blowup of coefficients.

[rbeezer]

Replying to dimpase:

For QQ and other infinite exact fields, naive Gauss elimination is not a good strategy all together, as it can lead to an exponential blowup of coefficients.

Hi Dima!

Agreed. This is part of the challenge in testing this. But I believe this is the current strategy. In other words, eventually Sage does echelon form, with or without this patch.

QQ gets converted to ZZ, and I do not know if the routines over ZZ control for this. But when I test with matrices containing number fields (cyclotomics), the rational coefficients do get out of hand.

So this patch is an incremental improvement. With or without, we still rely on Gaussian elimination in 100% of cases (or nearly so).

Rob