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

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]

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]

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]

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]

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]

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]

Get the second column of M:

sage: M[:,1]
[-2]
[ 8]
[ 1]
[ 2]

Get the first row of M:

sage: M[0,:]
[ 1 -2 -1 -1  9]

Get the last row of M (negative numbers count from the end):

sage: M[-1,:]
[-1  2 -2 -1  4]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

We save and load a matrix:

sage: A = matrix(Integers(8),3,range(9))
sage: loads(dumps(A)) == A
True

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

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

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