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: 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
>>> 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
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: 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]
>>> 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]
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 orNone
proof
– boolean orNone
(default:None
, seeproof.linear_algebra
orsage.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: 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)
>>> 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)
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()