Matrices over an arbitrary ring¶
AUTHORS:
William Stein
Martin Albrecht: conversion to Pyrex
Jaap Spies: various functions
Gary Zablackis: fixed a sign bug in generic determinant.
William Stein and Robert Bradshaw - complete restructuring.
Rob Beezer - refactor kernel functions.
Elements of matrix spaces are of class Matrix (or a
class derived from Matrix). They can be either sparse or dense, and
can be defined over any base ring.
EXAMPLES:
We create the \(2\times 3\) matrix
\[\begin{split}\left(\begin{matrix} 1&2&3\\4&5&6 \end{matrix}\right)\end{split}\]
as an element of a matrix space over \(\QQ\):
sage: M = MatrixSpace(QQ,2,3)
sage: A = M([1,2,3, 4,5,6]); A
[1 2 3]
[4 5 6]
sage: A.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field
>>> from sage.all import *
>>> M = MatrixSpace(QQ,Integer(2),Integer(3))
>>> A = M([Integer(1),Integer(2),Integer(3), Integer(4),Integer(5),Integer(6)]); A
[1 2 3]
[4 5 6]
>>> A.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field
M = MatrixSpace(QQ,2,3) A = M([1,2,3, 4,5,6]); A A.parent()
Alternatively, we could create A more directly as follows (which would completely avoid having to create the matrix space):
sage: A = matrix(QQ, 2, [1,2,3, 4,5,6]); A
[1 2 3]
[4 5 6]
>>> from sage.all import *
>>> A = matrix(QQ, Integer(2), [Integer(1),Integer(2),Integer(3), Integer(4),Integer(5),Integer(6)]); A
[1 2 3]
[4 5 6]
A = matrix(QQ, 2, [1,2,3, 4,5,6]); A
We next change the top-right entry of \(A\). Note that matrix indexing is \(0\)-based in Sage, so the top right entry is \((0,2)\), which should be thought of as “row number \(0\), column number \(2\)”.
sage: A[0,2] = 389
sage: A
[ 1 2 389]
[ 4 5 6]
>>> from sage.all import *
>>> A[Integer(0),Integer(2)] = Integer(389)
>>> A
[ 1 2 389]
[ 4 5 6]
A[0,2] = 389 A
Also notice how matrices print. All columns have the same width and entries in a given column are right justified. Next we compute the reduced row echelon form of \(A\).
sage: A.rref()
[ 1 0 -1933/3]
[ 0 1 1550/3]
>>> from sage.all import *
>>> A.rref()
[ 1 0 -1933/3]
[ 0 1 1550/3]
A.rref()
Indexing¶
Sage has quite flexible ways of extracting elements or submatrices from a matrix:
sage: m=[(1, -2, -1, -1,9), (1, 8, 6, 2,2), (1, 1, -1, 1,4), (-1, 2, -2, -1,4)] ; M = matrix(m)
sage: M
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
[-1 2 -2 -1 4]
>>> from sage.all import *
>>> m=[(Integer(1), -Integer(2), -Integer(1), -Integer(1),Integer(9)), (Integer(1), Integer(8), Integer(6), Integer(2),Integer(2)), (Integer(1), Integer(1), -Integer(1), Integer(1),Integer(4)), (-Integer(1), Integer(2), -Integer(2), -Integer(1),Integer(4))] ; M = matrix(m)
>>> M
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
[-1 2 -2 -1 4]
m=[(1, -2, -1, -1,9), (1, 8, 6, 2,2), (1, 1, -1, 1,4), (-1, 2, -2, -1,4)] ; M = matrix(m) M
Get the 2 x 2 submatrix of M, starting at row index and column index 1:
sage: M[1:3,1:3]
[ 8 6]
[ 1 -1]
>>> from sage.all import *
>>> M[Integer(1):Integer(3),Integer(1):Integer(3)]
[ 8 6]
[ 1 -1]
M[1:3,1:3]
Get the 2 x 3 submatrix of M starting at row index and column index 1:
sage: M[1:3,[1..3]]
[ 8 6 2]
[ 1 -1 1]
>>> from sage.all import *
>>> M[Integer(1):Integer(3),(ellipsis_range(Integer(1),Ellipsis,Integer(3)))]
[ 8 6 2]
[ 1 -1 1]
M[1:3,[1..3]]
Get the second column of M:
sage: M[:,1]
[-2]
[ 8]
[ 1]
[ 2]
>>> from sage.all import *
>>> M[:,Integer(1)]
[-2]
[ 8]
[ 1]
[ 2]
M[:,1]
Get the first row of M:
sage: M[0,:]
[ 1 -2 -1 -1 9]
>>> from sage.all import *
>>> M[Integer(0),:]
[ 1 -2 -1 -1 9]
M[0,:]
Get the last row of M (negative numbers count from the end):
sage: M[-1,:]
[-1 2 -2 -1 4]
>>> from sage.all import *
>>> M[-Integer(1),:]
[-1 2 -2 -1 4]
M[-1,:]
More examples:
sage: M[range(2),:]
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
sage: M[range(2),4]
[9]
[2]
sage: M[range(3),range(5)]
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
sage: M[3,range(5)]
[-1 2 -2 -1 4]
sage: M[3,:]
[-1 2 -2 -1 4]
sage: M[3,4]
4
sage: M[-1,:]
[-1 2 -2 -1 4]
sage: A = matrix(ZZ,3,4, [3, 2, -5, 0, 1, -1, 1, -4, 1, 0, 1, -3]); A
[ 3 2 -5 0]
[ 1 -1 1 -4]
[ 1 0 1 -3]
>>> from sage.all import *
>>> M[range(Integer(2)),:]
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
>>> M[range(Integer(2)),Integer(4)]
[9]
[2]
>>> M[range(Integer(3)),range(Integer(5))]
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
>>> M[Integer(3),range(Integer(5))]
[-1 2 -2 -1 4]
>>> M[Integer(3),:]
[-1 2 -2 -1 4]
>>> M[Integer(3),Integer(4)]
4
>>> M[-Integer(1),:]
[-1 2 -2 -1 4]
>>> A = matrix(ZZ,Integer(3),Integer(4), [Integer(3), Integer(2), -Integer(5), Integer(0), Integer(1), -Integer(1), Integer(1), -Integer(4), Integer(1), Integer(0), Integer(1), -Integer(3)]); A
[ 3 2 -5 0]
[ 1 -1 1 -4]
[ 1 0 1 -3]
M[range(2),:] M[range(2),4] M[range(3),range(5)] M[3,range(5)] M[3,:] M[3,4] M[-1,:] A = matrix(ZZ,3,4, [3, 2, -5, 0, 1, -1, 1, -4, 1, 0, 1, -3]); A
A series of three numbers, separated by colons, like n:m:s, means
numbers from n up to (but not including) m, in steps of s.
So 0:5:2 means the sequence [0,2,4]:
sage: A[:,0:4:2]
[ 3 -5]
[ 1 1]
[ 1 1]
sage: A[1:,0:4:2]
[1 1]
[1 1]
sage: A[2::-1,:]
[ 1 0 1 -3]
[ 1 -1 1 -4]
[ 3 2 -5 0]
sage: A[1:,3::-1]
[-4 1 -1 1]
[-3 1 0 1]
sage: A[1:,3::-2]
[-4 -1]
[-3 0]
sage: A[2::-1,3:1:-1]
[-3 1]
[-4 1]
[ 0 -5]
>>> from sage.all import *
>>> A[:,Integer(0):Integer(4):Integer(2)]
[ 3 -5]
[ 1 1]
[ 1 1]
>>> A[Integer(1):,Integer(0):Integer(4):Integer(2)]
[1 1]
[1 1]
>>> A[Integer(2)::-Integer(1),:]
[ 1 0 1 -3]
[ 1 -1 1 -4]
[ 3 2 -5 0]
>>> A[Integer(1):,Integer(3)::-Integer(1)]
[-4 1 -1 1]
[-3 1 0 1]
>>> A[Integer(1):,Integer(3)::-Integer(2)]
[-4 -1]
[-3 0]
>>> A[Integer(2)::-Integer(1),Integer(3):Integer(1):-Integer(1)]
[-3 1]
[-4 1]
[ 0 -5]
A[:,0:4:2] A[1:,0:4:2] A[2::-1,:] A[1:,3::-1] A[1:,3::-2] A[2::-1,3:1:-1]
We can also change submatrices using these indexing features:
sage: M=matrix([(1, -2, -1, -1,9), (1, 8, 6, 2,2), (1, 1, -1, 1,4), (-1, 2, -2, -1,4)]); M
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
[-1 2 -2 -1 4]
>>> from sage.all import *
>>> M=matrix([(Integer(1), -Integer(2), -Integer(1), -Integer(1),Integer(9)), (Integer(1), Integer(8), Integer(6), Integer(2),Integer(2)), (Integer(1), Integer(1), -Integer(1), Integer(1),Integer(4)), (-Integer(1), Integer(2), -Integer(2), -Integer(1),Integer(4))]); M
[ 1 -2 -1 -1 9]
[ 1 8 6 2 2]
[ 1 1 -1 1 4]
[-1 2 -2 -1 4]
M=matrix([(1, -2, -1, -1,9), (1, 8, 6, 2,2), (1, 1, -1, 1,4), (-1, 2, -2, -1,4)]); M
Set the 2 x 2 submatrix of M, starting at row index and column index 1:
sage: M[1:3,1:3] = [[1,0],[0,1]]; M
[ 1 -2 -1 -1 9]
[ 1 1 0 2 2]
[ 1 0 1 1 4]
[-1 2 -2 -1 4]
>>> from sage.all import *
>>> M[Integer(1):Integer(3),Integer(1):Integer(3)] = [[Integer(1),Integer(0)],[Integer(0),Integer(1)]]; M
[ 1 -2 -1 -1 9]
[ 1 1 0 2 2]
[ 1 0 1 1 4]
[-1 2 -2 -1 4]
M[1:3,1:3] = [[1,0],[0,1]]; M
Set the 2 x 3 submatrix of M starting at row index and column index 1:
sage: M[1:3,[1..3]] = M[2:4,0:3]; M
[ 1 -2 -1 -1 9]
[ 1 1 0 1 2]
[ 1 -1 2 -2 4]
[-1 2 -2 -1 4]
>>> from sage.all import *
>>> M[Integer(1):Integer(3),(ellipsis_range(Integer(1),Ellipsis,Integer(3)))] = M[Integer(2):Integer(4),Integer(0):Integer(3)]; M
[ 1 -2 -1 -1 9]
[ 1 1 0 1 2]
[ 1 -1 2 -2 4]
[-1 2 -2 -1 4]
M[1:3,[1..3]] = M[2:4,0:3]; M
Set part of the first column of M:
sage: M[1:,0]=[[2],[3],[4]]; M
[ 1 -2 -1 -1 9]
[ 2 1 0 1 2]
[ 3 -1 2 -2 4]
[ 4 2 -2 -1 4]
>>> from sage.all import *
>>> M[Integer(1):,Integer(0)]=[[Integer(2)],[Integer(3)],[Integer(4)]]; M
[ 1 -2 -1 -1 9]
[ 2 1 0 1 2]
[ 3 -1 2 -2 4]
[ 4 2 -2 -1 4]
M[1:,0]=[[2],[3],[4]]; M
Or do a similar thing with a vector:
sage: M[1:,0]=vector([-2,-3,-4]); M
[ 1 -2 -1 -1 9]
[-2 1 0 1 2]
[-3 -1 2 -2 4]
[-4 2 -2 -1 4]
>>> from sage.all import *
>>> M[Integer(1):,Integer(0)]=vector([-Integer(2),-Integer(3),-Integer(4)]); M
[ 1 -2 -1 -1 9]
[-2 1 0 1 2]
[-3 -1 2 -2 4]
[-4 2 -2 -1 4]
M[1:,0]=vector([-2,-3,-4]); M
Or a constant:
sage: M[1:,0]=30; M
[ 1 -2 -1 -1 9]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
>>> from sage.all import *
>>> M[Integer(1):,Integer(0)]=Integer(30); M
[ 1 -2 -1 -1 9]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
M[1:,0]=30; M
Set the first row of M:
sage: M[0,:]=[[20,21,22,23,24]]; M
[20 21 22 23 24]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
sage: M[0,:]=vector([0,1,2,3,4]); M
[ 0 1 2 3 4]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
sage: M[0,:]=-3; M
[-3 -3 -3 -3 -3]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
sage: A = matrix(ZZ,3,4, [3, 2, -5, 0, 1, -1, 1, -4, 1, 0, 1, -3]); A
[ 3 2 -5 0]
[ 1 -1 1 -4]
[ 1 0 1 -3]
>>> from sage.all import *
>>> M[Integer(0),:]=[[Integer(20),Integer(21),Integer(22),Integer(23),Integer(24)]]; M
[20 21 22 23 24]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
>>> M[Integer(0),:]=vector([Integer(0),Integer(1),Integer(2),Integer(3),Integer(4)]); M
[ 0 1 2 3 4]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
>>> M[Integer(0),:]=-Integer(3); M
[-3 -3 -3 -3 -3]
[30 1 0 1 2]
[30 -1 2 -2 4]
[30 2 -2 -1 4]
>>> A = matrix(ZZ,Integer(3),Integer(4), [Integer(3), Integer(2), -Integer(5), Integer(0), Integer(1), -Integer(1), Integer(1), -Integer(4), Integer(1), Integer(0), Integer(1), -Integer(3)]); A
[ 3 2 -5 0]
[ 1 -1 1 -4]
[ 1 0 1 -3]
M[0,:]=[[20,21,22,23,24]]; M M[0,:]=vector([0,1,2,3,4]); M M[0,:]=-3; M A = matrix(ZZ,3,4, [3, 2, -5, 0, 1, -1, 1, -4, 1, 0, 1, -3]); A
We can use the step feature of slices to set every other column:
sage: A[:,0:3:2] = 5; A
[ 5 2 5 0]
[ 5 -1 5 -4]
[ 5 0 5 -3]
sage: A[1:,0:4:2] = [[100,200],[300,400]]; A
[ 5 2 5 0]
[100 -1 200 -4]
[300 0 400 -3]
>>> from sage.all import *
>>> A[:,Integer(0):Integer(3):Integer(2)] = Integer(5); A
[ 5 2 5 0]
[ 5 -1 5 -4]
[ 5 0 5 -3]
>>> A[Integer(1):,Integer(0):Integer(4):Integer(2)] = [[Integer(100),Integer(200)],[Integer(300),Integer(400)]]; A
[ 5 2 5 0]
[100 -1 200 -4]
[300 0 400 -3]
A[:,0:3:2] = 5; A A[1:,0:4:2] = [[100,200],[300,400]]; A
We can also count backwards to flip the matrix upside down:
sage: A[::-1,:]=A; A
[300 0 400 -3]
[100 -1 200 -4]
[ 5 2 5 0]
sage: A[1:,3::-1]=[[2,3,0,1],[9,8,7,6]]; A
[300 0 400 -3]
[ 1 0 3 2]
[ 6 7 8 9]
sage: A[1:,::-2] = A[1:,::2]; A
[300 0 400 -3]
[ 1 3 3 1]
[ 6 8 8 6]
sage: A[::-1,3:1:-1] = [[4,3],[1,2],[-1,-2]]; A
[300 0 -2 -1]
[ 1 3 2 1]
[ 6 8 3 4]
>>> from sage.all import *
>>> A[::-Integer(1),:]=A; A
[300 0 400 -3]
[100 -1 200 -4]
[ 5 2 5 0]
>>> A[Integer(1):,Integer(3)::-Integer(1)]=[[Integer(2),Integer(3),Integer(0),Integer(1)],[Integer(9),Integer(8),Integer(7),Integer(6)]]; A
[300 0 400 -3]
[ 1 0 3 2]
[ 6 7 8 9]
>>> A[Integer(1):,::-Integer(2)] = A[Integer(1):,::Integer(2)]; A
[300 0 400 -3]
[ 1 3 3 1]
[ 6 8 8 6]
>>> A[::-Integer(1),Integer(3):Integer(1):-Integer(1)] = [[Integer(4),Integer(3)],[Integer(1),Integer(2)],[-Integer(1),-Integer(2)]]; A
[300 0 -2 -1]
[ 1 3 2 1]
[ 6 8 3 4]
A[::-1,:]=A; A A[1:,3::-1]=[[2,3,0,1],[9,8,7,6]]; A A[1:,::-2] = A[1:,::2]; A A[::-1,3:1:-1] = [[4,3],[1,2],[-1,-2]]; A
We save and load a matrix:
sage: A = matrix(Integers(8),3,range(9))
sage: loads(dumps(A)) == A
True
>>> from sage.all import *
>>> A = matrix(Integers(Integer(8)),Integer(3),range(Integer(9)))
>>> loads(dumps(A)) == A
True
A = matrix(Integers(8),3,range(9)) loads(dumps(A)) == A
MUTABILITY: Matrices are either immutable or not. When initially
created, matrices are typically mutable, so one can change their
entries. Once a matrix \(A\) is made immutable using
A.set_immutable() the entries of \(A\)
cannot be changed, and \(A\) can never be made mutable again.
However, properties of \(A\) such as its rank, characteristic
polynomial, etc., are all cached so computations involving
\(A\) may be more efficient. Once \(A\) is made
immutable it cannot be changed back. However, one can obtain a
mutable copy of \(A\) using copy(A).
EXAMPLES:
sage: A = matrix(RR,2,[1,10,3.5,2])
sage: A.set_immutable()
sage: copy(A) is A
False
>>> from sage.all import *
>>> A = matrix(RR,Integer(2),[Integer(1),Integer(10),RealNumber('3.5'),Integer(2)])
>>> A.set_immutable()
>>> copy(A) is A
False
A = matrix(RR,2,[1,10,3.5,2]) A.set_immutable() copy(A) is A
The echelon form method always returns immutable matrices with known rank.
EXAMPLES:
sage: A = matrix(Integers(8),3,range(9))
sage: A.determinant()
0
sage: A[0,0] = 5
sage: A.determinant()
1
sage: A.set_immutable()
sage: A[0,0] = 5
Traceback (most recent call last):
...
ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M).
>>> from sage.all import *
>>> A = matrix(Integers(Integer(8)),Integer(3),range(Integer(9)))
>>> A.determinant()
0
>>> A[Integer(0),Integer(0)] = Integer(5)
>>> A.determinant()
1
>>> A.set_immutable()
>>> A[Integer(0),Integer(0)] = Integer(5)
Traceback (most recent call last):
...
ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M).
A = matrix(Integers(8),3,range(9)) A.determinant() A[0,0] = 5 A.determinant() A.set_immutable() A[0,0] = 5
Implementation and Design¶
Class Diagram (an x means that class is currently supported):
x Matrix
x Matrix_sparse
x Matrix_generic_sparse
x Matrix_integer_sparse
x Matrix_rational_sparse
Matrix_cyclo_sparse
x Matrix_modn_sparse
Matrix_RR_sparse
Matrix_CC_sparse
Matrix_RDF_sparse
Matrix_CDF_sparse
x Matrix_dense
x Matrix_generic_dense
x Matrix_integer_dense
x Matrix_rational_dense
Matrix_cyclo_dense -- idea: restrict scalars to QQ, compute charpoly there, then factor
x Matrix_modn_dense
Matrix_RR_dense
Matrix_CC_dense
x Matrix_real_double_dense
x Matrix_complex_double_dense
x Matrix_complex_ball_dense
The corresponding files in the sage/matrix library code directory are named
[matrix] [base ring] [dense or sparse].
New matrices types can only be implemented in Cython.
*********** LEVEL 1 **********
NON-OPTIONAL
For each base field it is *absolutely* essential to completely
implement the following functionality for that base ring:
* __cinit__ -- should use check_allocarray from cysignals.memory
(only needed if allocate memory)
* __init__ -- this signature: 'def __init__(self, parent, entries, copy, coerce)'
* __dealloc__ -- use sig_free (only needed if allocate memory)
* set_unsafe(self, size_t i, size_t j, x) -- doesn't do bounds or any other checks; assumes x is in self._base_ring
* get_unsafe(self, size_t i, size_t j) -- doesn't do checks
* __richcmp__ -- always the same (I don't know why its needed -- bug in PYREX)
Note that the __init__ function must construct the all zero matrix if ``entries == None``.
*********** LEVEL 2 **********
IMPORTANT (and *highly* recommended):
After getting the special class with all level 1 functionality to
work, implement all of the following (they should not change
functionality, except speed (always faster!) in any way):
* def _pickle(self):
return data, version
* def _unpickle(self, data, int version)
reconstruct matrix from given data and version; may assume _parent, _nrows, and _ncols are set.
Use version numbers >= 0 so if you change the pickle strategy then
old objects still unpickle.
* cdef _list -- list of underlying elements (need not be a copy)
* cdef _dict -- sparse dictionary of underlying elements
* cdef _add_ -- add two matrices with identical parents
* _matrix_times_matrix_c_impl -- multiply two matrices with compatible dimensions and
identical base rings (both sparse or both dense)
* cpdef _richcmp_ -- compare two matrices with identical parents
* cdef _lmul_c_impl -- multiply this matrix on the right by a scalar, i.e., self * scalar
* cdef _rmul_c_impl -- multiply this matrix on the left by a scalar, i.e., scalar * self
* __copy__
* __neg__
The list and dict returned by _list and _dict will *not* be changed
by any internal algorithms and are not accessible to the user.
*********** LEVEL 3 **********
OPTIONAL:
* cdef _sub_
* __invert__
* _multiply_classical
* __deepcopy__
Further special support:
* Matrix windows -- to support Strassen multiplication for a given base ring.
* Other functions, e.g., transpose, for which knowing the
specific representation can be helpful.
.. NOTE::
- For caching, use self.fetch and self.cache.
- Any method that can change the matrix should call
``check_mutability()`` first. There are also many fast cdef'd bounds checking methods.
- Kernels of matrices
Implement only a left_kernel() or right_kernel() method, whichever requires
the least overhead (usually meaning little or no transposing). Let the
methods in the matrix2 class handle left, right, generic kernel distinctions.