Misc matrix algorithms¶

sage.matrix.misc.cmp_pivots(x, y)[source]¶

Compare two sequences of pivot columns.

If \(x\) is shorter than \(y\), return \(-1\), i.e., \(x < y\), “not as good”. If \(x\) is longer than \(y\), then \(x > y\), so “better” and return \(+1\). If the length is the same, then \(x\) is better, i.e., \(x > y\) if the entries of \(x\) are correspondingly \(\leq\) those of \(y\) with one being strictly less.

INPUT:

x, y – lists or tuples of integers

EXAMPLES:

We illustrate each of the above comparisons.

Sage

sage: from sage.matrix.misc import cmp_pivots
sage: cmp_pivots([1,2,3], [4,5,6,7])
-1
sage: cmp_pivots([1,2,3,5], [4,5,6])
1
sage: cmp_pivots([1,2,4], [1,2,3])
-1
sage: cmp_pivots([1,2,3], [1,2,3])
0
sage: cmp_pivots([1,2,3], [1,2,4])
1

Python

>>> from sage.all import *
>>> from sage.matrix.misc import cmp_pivots
>>> cmp_pivots([Integer(1),Integer(2),Integer(3)], [Integer(4),Integer(5),Integer(6),Integer(7)])
-1
>>> cmp_pivots([Integer(1),Integer(2),Integer(3),Integer(5)], [Integer(4),Integer(5),Integer(6)])
1
>>> cmp_pivots([Integer(1),Integer(2),Integer(4)], [Integer(1),Integer(2),Integer(3)])
-1
>>> cmp_pivots([Integer(1),Integer(2),Integer(3)], [Integer(1),Integer(2),Integer(3)])
0
>>> cmp_pivots([Integer(1),Integer(2),Integer(3)], [Integer(1),Integer(2),Integer(4)])
1

Sage Live

from sage.matrix.misc import cmp_pivots
cmp_pivots([1,2,3], [4,5,6,7])
cmp_pivots([1,2,3,5], [4,5,6])
cmp_pivots([1,2,4], [1,2,3])
cmp_pivots([1,2,3], [1,2,3])
cmp_pivots([1,2,3], [1,2,4])

sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, N)[source]¶

Given a sparse matrix over the integers and an integer modulus, do rational reconstruction on all entries of the matrix, viewed as numbers mod \(N\).

EXAMPLES:

Sage

sage: A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, 7, 2, 2, 3, 4, 3, 4, (5/7)%500], sparse=True)
sage: from sage.matrix.misc import matrix_integer_sparse_rational_reconstruction
sage: matrix_integer_sparse_rational_reconstruction(A, 500)
[1/3   2   3  -4]
[  7   2   2   3]
[  4   3   4 5/7]

Python

>>> from sage.all import *
>>> A = matrix(ZZ, Integer(3), Integer(4), [(Integer(1)/Integer(3))%Integer(500), Integer(2), Integer(3), (-Integer(4))%Integer(500), Integer(7), Integer(2), Integer(2), Integer(3), Integer(4), Integer(3), Integer(4), (Integer(5)/Integer(7))%Integer(500)], sparse=True)
>>> from sage.matrix.misc import matrix_integer_sparse_rational_reconstruction
>>> matrix_integer_sparse_rational_reconstruction(A, Integer(500))
[1/3   2   3  -4]
[  7   2   2   3]
[  4   3   4 5/7]

Sage Live

A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, 7, 2, 2, 3, 4, 3, 4, (5/7)%500], sparse=True)
from sage.matrix.misc import matrix_integer_sparse_rational_reconstruction
matrix_integer_sparse_rational_reconstruction(A, 500)

sage.matrix.misc.matrix_rational_echelon_form_multimodular(self, height_guess=None, proof=None)[source]¶

Return reduced row-echelon form using a multi-modular algorithm. Does not change self.

REFERENCE: Chapter 7 of Stein’s “Explicitly Computing Modular Forms”.

INPUT:

height_guess – integer or None; an estimate for \(H(dE)\), where \(E\) is the output echelon form and \(d\) is its denominator
proof – boolean or None (default: None, see proof.linear_algebra or sage.structure.proof). Note that the global Sage default is proof=True. With proof=False, the algorithm uses an explicitly supplied height_guess to attempt reconstruction early, subject to backend fallback, and checks each candidate first with the usual height certificate; if that fails, it checks the reduced-echelon shape and an exact row-space identity. The early attempt can reduce the number of modular images when the output is much smaller than the input, but failed reconstructions and the exact check can also make it slower

OUTPUT: a pair consisting of a matrix in echelon form and a tuple of pivot positions.

Note

A sparse input that is already close to dense, and wide enough that its echelon form can be far taller than itself, is delegated to dense FLINT: the optimized sparse modular backend is capped at small moduli and cannot always reach the required prime product. The result is still returned sparse, but height_guess and proof have no effect on that path, since FLINT is exact.

ALGORITHM:

The following is a modular algorithm for computing the echelon form. Define the height of a matrix to be the max of the absolute values of the entries.

Given Matrix A with n columns (self).

Rescale each row of the input matrix A independently to a primitive integer row. This does not change echelon form and makes reduction modulo lots of primes significantly easier if there were denominators. Henceforth we assume A has integer entries.

Let c be a guess for H(dE), where E is the output echelon form and d is its denominator. E.g., c=1000 if the matrix is very sparse and the application is to computing modular symbols.

Let M = n * c * H(A) + 1, where n is the number of columns of A. If proof=False and c was supplied explicitly, start instead with M = c + 1 in order to attempt rational reconstruction earlier.

List primes p_1, p_2, …, such that the product of the p_i is at least M. Sparse matrices ordinarily start with the primes supported by the optimized sparse modular backend. If M already exceeds the product of that complete finite range, start genuinely sparse inputs directly with larger primes and generic sparse matrices; otherwise switch to those matrices if the optimized range is exhausted later. When a wide input is already at least three-quarters dense and a heuristic based on the target and a Hadamard bound indicates that the complete optimized prime range may be insufficient, use dense FLINT instead.

Try to compute the rational reconstruction CRT echelon form of A mod the product of the p_i. If rational reconstruction fails, compute more echelon forms modulo subsequent primes, and attempt again. Make sure to keep the result of CRT on the primes from before, so we don’t have to do that computation again. Let E be this matrix.

Compute the denominator d of E. Attempt to prove that the result is correct by checking that

H(d*E)*ncols(A)*H(A) < (prod of reduction primes)

where H denotes the height. With proof=False, if this certificate fails, check the exact identity

d*A = A[:, P] * (d*E)[:r]

where P is the tuple of r pivot columns, after explicitly checking that E has reduced-echelon shape for P. If the applicable check fails, do step 4 with more primes.

EXAMPLES:

Sage

sage: A = matrix(QQ, 3, 7, [1..21])
sage: from sage.matrix.misc import matrix_rational_echelon_form_multimodular
sage: E, pivots = matrix_rational_echelon_form_multimodular(A)
sage: E
[ 1  0 -1 -2 -3 -4 -5]
[ 0  1  2  3  4  5  6]
[ 0  0  0  0  0  0  0]
sage: pivots
(0, 1)

sage: A = matrix(QQ, 3, 4, [0,0] + [1..9] + [-1/2^20])
sage: E, pivots = matrix_rational_echelon_form_multimodular(A)
sage: E
[                1                 0                 0 -10485761/1048576]
[                0                 1                 0  27262979/4194304]
[                0                 0                 1                 2]
sage: pivots
(0, 1, 2)

sage: A.echelon_form()
[                1                 0                 0 -10485761/1048576]
[                0                 1                 0  27262979/4194304]
[                0                 0                 1                 2]
sage: A.pivots()
(0, 1, 2)

Python

>>> from sage.all import *
>>> A = matrix(QQ, Integer(3), Integer(7), (ellipsis_range(Integer(1),Ellipsis,Integer(21))))
>>> from sage.matrix.misc import matrix_rational_echelon_form_multimodular
>>> E, pivots = matrix_rational_echelon_form_multimodular(A)
>>> E
[ 1  0 -1 -2 -3 -4 -5]
[ 0  1  2  3  4  5  6]
[ 0  0  0  0  0  0  0]
>>> pivots
(0, 1)

>>> A = matrix(QQ, Integer(3), Integer(4), [Integer(0),Integer(0)] + (ellipsis_range(Integer(1),Ellipsis,Integer(9))) + [-Integer(1)/Integer(2)**Integer(20)])
>>> E, pivots = matrix_rational_echelon_form_multimodular(A)
>>> E
[                1                 0                 0 -10485761/1048576]
[                0                 1                 0  27262979/4194304]
[                0                 0                 1                 2]
>>> pivots
(0, 1, 2)

>>> A.echelon_form()
[                1                 0                 0 -10485761/1048576]
[                0                 1                 0  27262979/4194304]
[                0                 0                 1                 2]
>>> A.pivots()
(0, 1, 2)

Sage Live

A = matrix(QQ, 3, 7, [1..21])
from sage.matrix.misc import matrix_rational_echelon_form_multimodular
E, pivots = matrix_rational_echelon_form_multimodular(A)
E
pivots
A = matrix(QQ, 3, 4, [0,0] + [1..9] + [-1/2^20])
E, pivots = matrix_rational_echelon_form_multimodular(A)
E
pivots
A.echelon_form()
A.pivots()

A small benchmark, showing that flint fraction-free multimodular algorithm is always faster than the fraction-free multimodular algorithm implemented in Python:

Sage

sage: import copy
sage: def benchmark(num_row, num_col, entry_size, timeout=2, integer_coefficient=True):
....:     A = matrix(QQ, [[
....:         randint(1, 2^entry_size) if integer_coefficient else ZZ(randint(1, 2^entry_size))/randint(1, 2^entry_size)
....:         for col in range(num_col)] for row in range(num_row)])
....:     data=[]
....:     for algorithm in ("flint:fflu", "flint:multimodular", "padic", "multimodular"):
....:         # classical is too slow
....:         B = copy.copy(A)
....:         t = walltime()
....:         alarm(timeout)
....:         try:
....:             B.echelonize(algorithm=algorithm)
....:         except AlarmInterrupt:
....:             pass
....:         finally:
....:             cancel_alarm()
....:         data.append((round(walltime(t), 4), algorithm))
....:     return sorted(data)
sage: benchmark(20, 20, 10000)  # long time
[...'flint:multimodular'...'multimodular'...'flint:fflu'...]
sage: benchmark(39, 40, 200)  # long time
[...'flint:multimodular'...'flint:fflu'...'multimodular'...]

Python

>>> from sage.all import *
>>> import copy
>>> def benchmark(num_row, num_col, entry_size, timeout=Integer(2), integer_coefficient=True):
...     A = matrix(QQ, [[
...         randint(Integer(1), Integer(2)**entry_size) if integer_coefficient else ZZ(randint(Integer(1), Integer(2)**entry_size))/randint(Integer(1), Integer(2)**entry_size)
...         for col in range(num_col)] for row in range(num_row)])
...     data=[]
...     for algorithm in ("flint:fflu", "flint:multimodular", "padic", "multimodular"):
...         # classical is too slow
...         B = copy.copy(A)
...         t = walltime()
...         alarm(timeout)
...         try:
...             B.echelonize(algorithm=algorithm)
...         except AlarmInterrupt:
...             pass
...         finally:
...             cancel_alarm()
...         data.append((round(walltime(t), Integer(4)), algorithm))
...     return sorted(data)
>>> benchmark(Integer(20), Integer(20), Integer(10000))  # long time
[...'flint:multimodular'...'multimodular'...'flint:fflu'...]
>>> benchmark(Integer(39), Integer(40), Integer(200))  # long time
[...'flint:multimodular'...'flint:fflu'...'multimodular'...]

Sage Live

import copy
def benchmark(num_row, num_col, entry_size, timeout=2, integer_coefficient=True):
    A = matrix(QQ, [[
        randint(1, 2^entry_size) if integer_coefficient else ZZ(randint(1, 2^entry_size))/randint(1, 2^entry_size)
        for col in range(num_col)] for row in range(num_row)])
    data=[]
    for algorithm in ("flint:fflu", "flint:multimodular", "padic", "multimodular"):
        # classical is too slow
        B = copy.copy(A)
        t = walltime()
        alarm(timeout)
        try:
            B.echelonize(algorithm=algorithm)
        except AlarmInterrupt:
            pass
        finally:
            cancel_alarm()
        data.append((round(walltime(t), 4), algorithm))
    return sorted(data)
benchmark(20, 20, 10000)  # long time
benchmark(39, 40, 200)  # long time

In older versions of flint before this issue is fixed, algorithm='flint' (automatic choice) may be slower than algorithm='flint:multimodular'.

In this case, there are more columns than rows, which means the resulting matrix has height much higher than the input matrix. We check that the function does not take too long:

Sage

sage: A = matrix(QQ, [[randint(1, 2^500) for col in range(40)] for row in range(20)])
sage: t = walltime()
sage: A.echelonize(algorithm="multimodular")  # long time
sage: t = walltime(t)  # long time
sage: (t < 10, t)  # long time
(True, ...)

Python

>>> from sage.all import *
>>> A = matrix(QQ, [[randint(Integer(1), Integer(2)**Integer(500)) for col in range(Integer(40))] for row in range(Integer(20))])
>>> t = walltime()
>>> A.echelonize(algorithm="multimodular")  # long time
>>> t = walltime(t)  # long time
>>> (t < Integer(10), t)  # long time
(True, ...)

Sage Live

A = matrix(QQ, [[randint(1, 2^500) for col in range(40)] for row in range(20)])
t = walltime()
A.echelonize(algorithm="multimodular")  # long time
t = walltime(t)  # long time
(t < 10, t)  # long time