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
proof – boolean or None (default: None, see proof.linear_algebra or sage.structure.proof). Note that the global Sage default is proof=True

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

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 input matrix A to have integer entries. 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 the height of the echelon form. E.g., c=1000, e.g., if matrix is very sparse and application is to computing modular symbols.

Let M = n * c * H(A) + 1, where n is the number of columns of A.

List primes p_1, p_2, …, such that the product of the p_i is at least M.

Try to compute the rational reconstruction CRT echelon form of A mod the product of the p_i. If rational reconstruction fails, compute 1 more echelon forms mod the next prime, 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 result is correct by checking that

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

where H denotes the height. If this fails, do step 4 with a few 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