Misc matrix algorithms

Code goes here mainly when it needs access to the internal structure of several classes, and we want to avoid circular cimports.

NOTE: The whole problem of avoiding circular imports – the reason for existence of this file – is now a non-issue, since some bugs in Cython were fixed. Probably all this code should be moved into the relevant classes and this file deleted.

sage.matrix.misc.cmp_pivots(x, y)

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 <= 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.matrix.misc.cmp_pivots([1,2,3], [4,5,6,7])
-1
sage: sage.matrix.misc.cmp_pivots([1,2,3,5], [4,5,6])
1
sage: sage.matrix.misc.cmp_pivots([1,2,4], [1,2,3])
-1
sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,3])
0
sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,4])
1
sage.matrix.misc.hadamard_row_bound_mpfr(A)

Given a matrix A with entries that coerce to RR, compute the row Hadamard bound on the determinant.

INPUT:

A – a matrix over RR

OUTPUT:

integer – an integer n such that the absolute value of the
determinant of this matrix is at most $10^n$.

EXAMPLES:

We create a very large matrix, compute the row Hadamard bound, and also compute the row Hadamard bound of the transpose, which happens to be sharp.

sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292])
sage: import sage.matrix.misc
sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.change_ring(RR))
13976
sage: len(str(a.det()))
12215
sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.transpose().change_ring(RR))
12215

Note that in the above example using RDF would overflow:

sage: b = a.change_ring(RDF)
sage: b._hadamard_row_bound()
Traceback (most recent call last):
...
OverflowError: cannot convert float infinity to integer
sage.matrix.misc.matrix_integer_dense_rational_reconstruction(A, N)

Given a matrix over the integers and an integer modulus, do rational reconstruction on all entries of the matrix, viewed as numbers mod N. This is done efficiently by assuming there is a large common factor dividing the denominators.

INPUT:

A – matrix N – an integer

EXAMPLES:

sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ)
sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(B, 500)
[ 1/3  2/3    1 -4/3]
[ 7/3  2/3    6    1]
[ 4/3    1  4/3  5/3]
sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, N)

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: 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: sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, 500)
[1/3   2   3  -4]
[  7   2   2   3]
[  4   3   4 5/7]
sage.matrix.misc.matrix_rational_echelon_form_multimodular(self, height_guess=None, proof=None)

Returns 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).

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

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

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

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

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

  6. 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: 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)