Base class for matrices, part 2

For design documentation see matrix/docs.py.

AUTHORS:

  • William Stein: initial version

  • Jaap Spies (2006-02-24): added prod_of_row_sums, permanent, permanental_minor, rook_vector methods

  • Robert Bradshaw (2007-06-14): added subdivide method

  • Jaap Spies (2007-11-14): implemented _binomial, _choose auxiliary functions

  • William Stein (2007-11-18): added _gram_schmidt_noscale method

  • David Loeffler (2008-12-05): added smith_form method

  • David Loeffler (2009-06-01): added _echelon_form_PID method

  • Sebastian Pancratz (2009-06-25): implemented adjoint and charpoly methods; fixed adjoint reflecting the change that _adjoint is now implemented in Matrix; used the division-free algorithm for charpoly

  • Rob Beezer (2009-07-13): added elementwise_product method

  • Miguel Marco (2010-06-19): modified eigenvalues and eigenvectors functions to allow the option extend=False

  • Thierry Monteil (2010-10-05): bugfix for Issue #10063, so that the determinant is computed even for rings for which the is_field method is not implemented.

  • Rob Beezer (2010-12-13): added conjugate_transpose method

  • Rob Beezer (2011-02-05): refactored all of the matrix kernel routines; added extended_echelon_form, right_kernel_matrix, QR, _gram_schmidt_noscale, is_similar methods

  • Moritz Minzlaff (2011-03-17): corrected _echelon_form_PID method for matrices of one row, fixed in Issue #9053

  • Rob Beezer (2011-06-09): added is_normal, is_diagonalizable, LU, cyclic_subspace, zigzag_form, rational_form methods

  • Rob Beezer (2012-05-27): added indefinite_factorization, is_positive_definite, cholesky methods

  • Darij Grinberg (2013-10-01): added first (slow) pfaffian implementation

  • Mario Pernici (2014-07-01): modified rook_vector method

  • Rob Beezer (2015-05-25): modified is_similar method

  • Samuel Lelièvre (2020-09-18): improved method LLL_gram based on a patch by William Stein posted at Issue #5178, moving the method from its initial location in sage.matrix.integer_matrix_dense

  • Michael Jung (2020-10-02): added Bär-Faddeev-LeVerrier algorithm for the Pfaffian

  • Moritz Firsching(2020-10-05): added quantum_determinant

  • Dima Pasechnik (2022-11-08): fixed echelonize for inexact matrices

class sage.matrix.matrix2.Matrix[source]

Bases: Matrix

Base class for matrices, part 2.

C

Return the conjugate matrix.

EXAMPLES:

sage: A = matrix(QQbar, [[     -3,  5 - 3*I, 7 - 4*I],                      # needs sage.rings.number_field
....:                    [7 + 3*I, -1 + 6*I, 3 + 5*I],
....:                    [3 + 3*I, -3 + 6*I, 5 +   I]])
sage: A.C                                                                   # needs sage.rings.number_field
[      -3  5 + 3*I  7 + 4*I]
[ 7 - 3*I -1 - 6*I  3 - 5*I]
[ 3 - 3*I -3 - 6*I  5 - 1*I]
>>> from sage.all import *
>>> A = matrix(QQbar, [[     -Integer(3),  Integer(5) - Integer(3)*I, Integer(7) - Integer(4)*I],                      # needs sage.rings.number_field
...                    [Integer(7) + Integer(3)*I, -Integer(1) + Integer(6)*I, Integer(3) + Integer(5)*I],
...                    [Integer(3) + Integer(3)*I, -Integer(3) + Integer(6)*I, Integer(5) +   I]])
>>> A.C                                                                   # needs sage.rings.number_field
[      -3  5 + 3*I  7 + 4*I]
[ 7 - 3*I -1 - 6*I  3 - 5*I]
[ 3 - 3*I -3 - 6*I  5 - 1*I]
H

Return the conjugate-transpose (Hermitian) matrix.

EXAMPLES:

sage: A = matrix(QQbar, [[     -3,  5 - 3*I, 7 - 4*I],                      # needs sage.rings.number_field
....:                    [7 + 3*I, -1 + 6*I, 3 + 5*I],
....:                    [3 + 3*I, -3 + 6*I, 5 +   I]])
sage: A.H                                                                   # needs sage.rings.number_field
[      -3  7 - 3*I  3 - 3*I]
[ 5 + 3*I -1 - 6*I -3 - 6*I]
[ 7 + 4*I  3 - 5*I  5 - 1*I]
>>> from sage.all import *
>>> A = matrix(QQbar, [[     -Integer(3),  Integer(5) - Integer(3)*I, Integer(7) - Integer(4)*I],                      # needs sage.rings.number_field
...                    [Integer(7) + Integer(3)*I, -Integer(1) + Integer(6)*I, Integer(3) + Integer(5)*I],
...                    [Integer(3) + Integer(3)*I, -Integer(3) + Integer(6)*I, Integer(5) +   I]])
>>> A.H                                                                   # needs sage.rings.number_field
[      -3  7 - 3*I  3 - 3*I]
[ 5 + 3*I -1 - 6*I -3 - 6*I]
[ 7 + 4*I  3 - 5*I  5 - 1*I]
LLL_gram(flag=0)[source]

Return the LLL transformation matrix for this Gram matrix.

That is, the transformation matrix U over ZZ of determinant 1 that transforms the lattice with this matrix as Gram matrix to a lattice that is LLL-reduced.

Always works when self is positive definite, might work in some semidefinite and indefinite cases.

INPUT:

  • self – the Gram matrix of a quadratic form or of a lattice equipped with a bilinear form

  • flag – an optional flag passed to qflllgram According to pari:qflllgram’s documentation the options are:

    • 0 – (default), assume that self has either exact (integral or rational) or real floating point entries. The matrix is rescaled, converted to integers and the behavior is then as in flag=1.

    • 1 – assume that G is integral. Computations involving Gram-Schmidt vectors are approximate, with precision varying as needed.

OUTPUT:

A dense matrix U over the integers with determinant 1 such that U.T * M * U is LLL-reduced.

ALGORITHM:

Calls PARI’s pari:qflllgram.

EXAMPLES:

Create a Gram matrix and LLL-reduce it:

sage: M = Matrix(ZZ, 2, 2, [5, 3, 3, 2])
sage: U = M.LLL_gram()                                                      # needs sage.libs.pari
sage: MM = U.transpose() * M * U                                            # needs sage.libs.pari
sage: M, U, MM                                                              # needs sage.libs.pari
(
[5 3]  [-1  1]  [1 0]
[3 2], [ 1 -2], [0 1]
)
>>> from sage.all import *
>>> M = Matrix(ZZ, Integer(2), Integer(2), [Integer(5), Integer(3), Integer(3), Integer(2)])
>>> U = M.LLL_gram()                                                      # needs sage.libs.pari
>>> MM = U.transpose() * M * U                                            # needs sage.libs.pari
>>> M, U, MM                                                              # needs sage.libs.pari
(
[5 3]  [-1  1]  [1 0]
[3 2], [ 1 -2], [0 1]
)

For a Gram matrix over RR with a length one first vector and a very short second vector, the LLL-reduced basis is obtained by swapping the two basis vectors (and changing sign to preserve orientation).

sage: M = Matrix(RDF, 2, 2, [1, 0, 0, 1e-5])
sage: M.LLL_gram()                                                          # needs sage.libs.pari
[ 0 -1]
[ 1  0]
>>> from sage.all import *
>>> M = Matrix(RDF, Integer(2), Integer(2), [Integer(1), Integer(0), Integer(0), RealNumber('1e-5')])
>>> M.LLL_gram()                                                          # needs sage.libs.pari
[ 0 -1]
[ 1  0]

The algorithm might work for some semidefinite and indefinite forms:

sage: Matrix(ZZ, 2, 2, [2, 6, 6, 3]).LLL_gram()                             # needs sage.libs.pari
[-3 -1]
[ 1  0]
sage: Matrix(ZZ, 2, 2, [1, 0, 0, -1]).LLL_gram()                            # needs sage.libs.pari
[ 0 -1]
[ 1  0]
>>> from sage.all import *
>>> Matrix(ZZ, Integer(2), Integer(2), [Integer(2), Integer(6), Integer(6), Integer(3)]).LLL_gram()                             # needs sage.libs.pari
[-3 -1]
[ 1  0]
>>> Matrix(ZZ, Integer(2), Integer(2), [Integer(1), Integer(0), Integer(0), -Integer(1)]).LLL_gram()                            # needs sage.libs.pari
[ 0 -1]
[ 1  0]

However, it might fail for others, either raising a ValueError:

sage: Matrix(ZZ, 1, 1, [0]).LLL_gram()                                      # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram did not return a square matrix,
perhaps the matrix is not positive definite

sage: Matrix(ZZ, 2, 2, [0, 1, 1, 0]).LLL_gram()                             # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram did not return a square matrix,
perhaps the matrix is not positive definite
>>> from sage.all import *
>>> Matrix(ZZ, Integer(1), Integer(1), [Integer(0)]).LLL_gram()                                      # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram did not return a square matrix,
perhaps the matrix is not positive definite

>>> Matrix(ZZ, Integer(2), Integer(2), [Integer(0), Integer(1), Integer(1), Integer(0)]).LLL_gram()                             # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram did not return a square matrix,
perhaps the matrix is not positive definite

or running forever:

sage: Matrix(ZZ, 2, 2, [-5, -1, -1, -5]).LLL_gram()         # not tested, needs sage.libs.pari
Traceback (most recent call last):
...
RuntimeError: infinite loop while calling qflllgram
>>> from sage.all import *
>>> Matrix(ZZ, Integer(2), Integer(2), [-Integer(5), -Integer(1), -Integer(1), -Integer(5)]).LLL_gram()         # not tested, needs sage.libs.pari
Traceback (most recent call last):
...
RuntimeError: infinite loop while calling qflllgram

Nonreal input leads to a value error:

sage: Matrix(2, 2, [CDF(1, 1), 0, 0, 1]).LLL_gram()                          # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram failed, perhaps the matrix is not positive definite
>>> from sage.all import *
>>> Matrix(Integer(2), Integer(2), [CDF(Integer(1), Integer(1)), Integer(0), Integer(0), Integer(1)]).LLL_gram()                          # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: qflllgram failed, perhaps the matrix is not positive definite
LU(pivot=None, format='plu')[source]

Finds a decomposition into a lower-triangular matrix and an upper-triangular matrix.

INPUT:

  • pivot – pivoting strategy

    • ‘auto’ (default) – see if the matrix entries are ordered (i.e. if they have an absolute value method), and if so, use a the partial pivoting strategy. Otherwise, fall back to the nonzero strategy. This is the best choice for general routines that may call this for matrix entries of a variety of types.

    • 'partial' – each column is examined for the element with the largest absolute value and the row containing this element is swapped into place

    • 'nonzero' – the first nonzero element in a column is located and the row with this element is used

  • format – contents of output, see more discussion below about output

    • ‘plu’ (default) – a triple; matrices P, L and U such that A = P*L*U

    • 'compact' – a pair; row permutation as a tuple, and the matrices L and U combined into one matrix

OUTPUT:

Suppose that \(A\) is an \(m\times n\) matrix, then an LU decomposition is a lower-triangular \(m\times m\) matrix \(L\) with every diagonal element equal to 1, and an upper-triangular \(m\times n\) matrix, \(U\) such that the product \(LU\), after a permutation of the rows, is then equal to \(A\). For the ‘plu’ format the permutation is returned as an \(m\times m\) permutation matrix \(P\) such that

\[A = PLU\]

It is more common to place the permutation matrix just to the left of \(A\). If you desire this version, then use the inverse of \(P\) which is computed most efficiently as its transpose.

If the ‘partial’ pivoting strategy is used, then the non-diagonal entries of \(L\) will be less than or equal to 1 in absolute value. The ‘nonzero’ pivot strategy may be faster, but the growth of data structures for elements of the decomposition might counteract the advantage.

By necessity, returned matrices have a base ring equal to the fraction field of the base ring of the original matrix.

In the ‘compact’ format, the first returned value is a tuple that is a permutation of the rows of \(LU\) that yields \(A\). See the doctest for how you might employ this permutation. Then the matrices \(L\) and \(U\) are merged into one matrix – remove the diagonal of ones in \(L\) and the remaining nonzero entries can replace the entries of \(U\) beneath the diagonal.

The results are cached, only in the compact format, separately for each pivot strategy called. Repeated requests for the ‘plu’ format will require just a small amount of overhead in each call to bust out the compact format to the three matrices. Since only the compact format is cached, the components of the compact format are immutable, while the components of the ‘plu’ format are regenerated, and hence are mutable.

Notice that while \(U\) is similar to row-echelon form and the rows of \(U\) span the row space of \(A\), the rows of \(U\) are not generally linearly independent. Nor are the pivot columns (or rank) immediately obvious. However for rings without specialized echelon form routines, this method is about twice as fast as the generic echelon form routine since it only acts “below the diagonal”, as would be predicted from a theoretical analysis of the algorithms.

Note

This is an exact computation, so limited to exact rings. If you need numerical results, convert the base ring to the field of real double numbers, RDF or the field of complex double numbers, CDF, which will use a faster routine that is careful about numerical subtleties.

ALGORITHM:

“Gaussian Elimination with Partial Pivoting,” Algorithm 21.1 of [TB1997].

EXAMPLES:

Notice the difference in the \(L\) matrix as a result of different pivoting strategies. With partial pivoting, every entry of \(L\) has absolute value 1 or less.

sage: A = matrix(QQ, [[1, -1,  0,  2,  4,  7, -1],
....:                 [2, -1,  0,  6,  4,  8, -2],
....:                 [2,  0,  1,  4,  2,  6,  0],
....:                 [1,  0, -1,  8, -1, -1, -3],
....:                 [1,  1,  2, -2, -1,  1,  3]])
sage: P, L, U = A.LU(pivot='partial')
sage: P
[0 0 0 0 1]
[1 0 0 0 0]
[0 0 0 1 0]
[0 0 1 0 0]
[0 1 0 0 0]
sage: L
[   1    0    0    0    0]
[ 1/2    1    0    0    0]
[ 1/2  1/3    1    0    0]
[   1  2/3  1/5    1    0]
[ 1/2 -1/3 -2/5    0    1]
sage: U
[    2    -1     0     6     4     8    -2]
[    0   3/2     2    -5    -3    -3     4]
[    0     0  -5/3  20/3    -2    -4 -10/3]
[    0     0     0     0   2/5   4/5     0]
[    0     0     0     0   1/5   2/5     0]
sage: A == P*L*U
True
sage: P, L, U = A.LU(pivot='nonzero')
sage: P
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
sage: L
[ 1  0  0  0  0]
[ 2  1  0  0  0]
[ 2  2  1  0  0]
[ 1  1 -1  1  0]
[ 1  2  2  0  1]
sage: U
[ 1 -1  0  2  4  7 -1]
[ 0  1  0  2 -4 -6  0]
[ 0  0  1 -4  2  4  2]
[ 0  0  0  0  1  2  0]
[ 0  0  0  0 -1 -2  0]
sage: A == P*L*U
True
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(1), -Integer(1),  Integer(0),  Integer(2),  Integer(4),  Integer(7), -Integer(1)],
...                 [Integer(2), -Integer(1),  Integer(0),  Integer(6),  Integer(4),  Integer(8), -Integer(2)],
...                 [Integer(2),  Integer(0),  Integer(1),  Integer(4),  Integer(2),  Integer(6),  Integer(0)],
...                 [Integer(1),  Integer(0), -Integer(1),  Integer(8), -Integer(1), -Integer(1), -Integer(3)],
...                 [Integer(1),  Integer(1),  Integer(2), -Integer(2), -Integer(1),  Integer(1),  Integer(3)]])
>>> P, L, U = A.LU(pivot='partial')
>>> P
[0 0 0 0 1]
[1 0 0 0 0]
[0 0 0 1 0]
[0 0 1 0 0]
[0 1 0 0 0]
>>> L
[   1    0    0    0    0]
[ 1/2    1    0    0    0]
[ 1/2  1/3    1    0    0]
[   1  2/3  1/5    1    0]
[ 1/2 -1/3 -2/5    0    1]
>>> U
[    2    -1     0     6     4     8    -2]
[    0   3/2     2    -5    -3    -3     4]
[    0     0  -5/3  20/3    -2    -4 -10/3]
[    0     0     0     0   2/5   4/5     0]
[    0     0     0     0   1/5   2/5     0]
>>> A == P*L*U
True
>>> P, L, U = A.LU(pivot='nonzero')
>>> P
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
>>> L
[ 1  0  0  0  0]
[ 2  1  0  0  0]
[ 2  2  1  0  0]
[ 1  1 -1  1  0]
[ 1  2  2  0  1]
>>> U
[ 1 -1  0  2  4  7 -1]
[ 0  1  0  2 -4 -6  0]
[ 0  0  1 -4  2  4  2]
[ 0  0  0  0  1  2  0]
[ 0  0  0  0 -1 -2  0]
>>> A == P*L*U
True

An example of the compact format.

sage: B = matrix(QQ, [[ 1,  3,  5,  5],
....:                 [ 1,  4,  7,  8],
....:                 [-1, -4, -6, -6],
....:                 [ 0, -2, -5, -8],
....:                 [-2, -6, -6, -2]])
sage: perm, M = B.LU(format='compact')
sage: perm
(4, 3, 0, 1, 2)
sage: M
[  -2   -6   -6   -2]
[   0   -2   -5   -8]
[-1/2    0    2    4]
[-1/2 -1/2  3/4    0]
[ 1/2  1/2 -1/4    0]
>>> from sage.all import *
>>> B = matrix(QQ, [[ Integer(1),  Integer(3),  Integer(5),  Integer(5)],
...                 [ Integer(1),  Integer(4),  Integer(7),  Integer(8)],
...                 [-Integer(1), -Integer(4), -Integer(6), -Integer(6)],
...                 [ Integer(0), -Integer(2), -Integer(5), -Integer(8)],
...                 [-Integer(2), -Integer(6), -Integer(6), -Integer(2)]])
>>> perm, M = B.LU(format='compact')
>>> perm
(4, 3, 0, 1, 2)
>>> M
[  -2   -6   -6   -2]
[   0   -2   -5   -8]
[-1/2    0    2    4]
[-1/2 -1/2  3/4    0]
[ 1/2  1/2 -1/4    0]

We can easily illustrate the relationships between the two formats with a square matrix.

sage: C = matrix(QQ, [[-2,  3, -2, -5],
....:                 [ 1, -2,  1,  3],
....:                 [-4,  7, -3, -8],
....:                 [-3,  8, -1, -5]])
sage: P, L, U = C.LU(format='plu')
sage: perm, M = C.LU(format='compact')
sage: (L - identity_matrix(4)) + U == M
True
sage: p = [perm[i]+1 for i in range(len(perm))]
sage: PP = Permutation(p).to_matrix()
sage: PP == P
True
>>> from sage.all import *
>>> C = matrix(QQ, [[-Integer(2),  Integer(3), -Integer(2), -Integer(5)],
...                 [ Integer(1), -Integer(2),  Integer(1),  Integer(3)],
...                 [-Integer(4),  Integer(7), -Integer(3), -Integer(8)],
...                 [-Integer(3),  Integer(8), -Integer(1), -Integer(5)]])
>>> P, L, U = C.LU(format='plu')
>>> perm, M = C.LU(format='compact')
>>> (L - identity_matrix(Integer(4))) + U == M
True
>>> p = [perm[i]+Integer(1) for i in range(len(perm))]
>>> PP = Permutation(p).to_matrix()
>>> PP == P
True

For a nonsingular matrix, and the ‘nonzero’ pivot strategy there is no need to permute rows, so the permutation matrix will be the identity. Furthermore, it can be shown that then the \(L\) and \(U\) matrices are uniquely determined by requiring \(L\) to have ones on the diagonal.

sage: D = matrix(QQ, [[ 1,  0,  2,  0, -2, -1],
....:                 [ 3, -2,  3, -1,  0,  6],
....:                 [-4,  2, -3,  1, -1, -8],
....:                 [-2,  2, -3,  2,  1,  0],
....:                 [ 0, -1, -1,  0,  2,  5],
....:                 [-1,  2, -4, -1,  5, -3]])
sage: P, L, U = D.LU(pivot='nonzero')
sage: P
[1 0 0 0 0 0]
[0 1 0 0 0 0]
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
sage: L
[   1    0    0    0    0    0]
[   3    1    0    0    0    0]
[  -4   -1    1    0    0    0]
[  -2   -1   -1    1    0    0]
[   0  1/2  1/4  1/2    1    0]
[  -1   -1 -5/2   -2   -6    1]
sage: U
[   1    0    2    0   -2   -1]
[   0   -2   -3   -1    6    9]
[   0    0    2    0   -3   -3]
[   0    0    0    1    0    4]
[   0    0    0    0 -1/4 -3/4]
[   0    0    0    0    0    1]
sage: D == L*U
True
>>> from sage.all import *
>>> D = matrix(QQ, [[ Integer(1),  Integer(0),  Integer(2),  Integer(0), -Integer(2), -Integer(1)],
...                 [ Integer(3), -Integer(2),  Integer(3), -Integer(1),  Integer(0),  Integer(6)],
...                 [-Integer(4),  Integer(2), -Integer(3),  Integer(1), -Integer(1), -Integer(8)],
...                 [-Integer(2),  Integer(2), -Integer(3),  Integer(2),  Integer(1),  Integer(0)],
...                 [ Integer(0), -Integer(1), -Integer(1),  Integer(0),  Integer(2),  Integer(5)],
...                 [-Integer(1),  Integer(2), -Integer(4), -Integer(1),  Integer(5), -Integer(3)]])
>>> P, L, U = D.LU(pivot='nonzero')
>>> P
[1 0 0 0 0 0]
[0 1 0 0 0 0]
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
>>> L
[   1    0    0    0    0    0]
[   3    1    0    0    0    0]
[  -4   -1    1    0    0    0]
[  -2   -1   -1    1    0    0]
[   0  1/2  1/4  1/2    1    0]
[  -1   -1 -5/2   -2   -6    1]
>>> U
[   1    0    2    0   -2   -1]
[   0   -2   -3   -1    6    9]
[   0    0    2    0   -3   -3]
[   0    0    0    1    0    4]
[   0    0    0    0 -1/4 -3/4]
[   0    0    0    0    0    1]
>>> D == L*U
True

The base ring of the matrix may be any field, or a ring which has a fraction field implemented in Sage. The ring needs to be exact (there is a numerical LU decomposition for matrices over RDF and CDF). Matrices returned are over the original field, or the fraction field of the ring. If the field is not ordered (i.e. the absolute value function is not implemented), then the pivot strategy needs to be ‘nonzero’.

sage: A = matrix(RealField(100), 3, 3, range(9))
sage: P, L, U = A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix must be exact,
not Real Field with 100 bits of precision

sage: A = matrix(Integers(6), 3, 2, range(6))
sage: A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix needs a field of fractions,
not Ring of integers modulo 6

sage: R.<y> = PolynomialRing(QQ, 'y')
sage: B = matrix(R, [[y+1, y^2+y], [y^2, y^3]])
sage: P, L, U = B.LU(pivot='partial')
Traceback (most recent call last):
...
TypeError: cannot take absolute value of matrix entries, try 'pivot=nonzero'
sage: P, L, U = B.LU(pivot='nonzero')
sage: P
[1 0]
[0 1]
sage: L
[          1           0]
[y^2/(y + 1)           1]
sage: U
[  y + 1 y^2 + y]
[      0       0]
sage: L.base_ring()
Fraction Field of Univariate Polynomial Ring in y over Rational Field
sage: B == P*L*U
True

sage: # needs sage.rings.finite_rings
sage: F.<a> = FiniteField(5^2)
sage: C = matrix(F, [[a + 3, 4*a + 4, 2, 4*a + 2],
....:                [3, 2*a + 4, 2*a + 4, 2*a + 1],
....:                [3*a + 1, a + 3, 2*a + 4, 4*a + 3],
....:                [a, 3, 3*a + 1, a]])
sage: P, L, U = C.LU(pivot='nonzero')
sage: P
[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
sage: L                                                                     # needs sage.combinat
[      1       0       0       0]
[3*a + 3       1       0       0]
[    2*a 4*a + 2       1       0]
[2*a + 3       2 2*a + 4       1]
sage: U                                                                     # needs sage.combinat
[  a + 3 4*a + 4       2 4*a + 2]
[      0   a + 1   a + 3 2*a + 4]
[      0       0       1 4*a + 2]
[      0       0       0       0]
sage: L.base_ring()                                                         # needs sage.combinat
Finite Field in a of size 5^2
sage: C == P*L*U                                                            # needs sage.combinat
True
>>> from sage.all import *
>>> A = matrix(RealField(Integer(100)), Integer(3), Integer(3), range(Integer(9)))
>>> P, L, U = A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix must be exact,
not Real Field with 100 bits of precision

>>> A = matrix(Integers(Integer(6)), Integer(3), Integer(2), range(Integer(6)))
>>> A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix needs a field of fractions,
not Ring of integers modulo 6

>>> R = PolynomialRing(QQ, 'y', names=('y',)); (y,) = R._first_ngens(1)
>>> B = matrix(R, [[y+Integer(1), y**Integer(2)+y], [y**Integer(2), y**Integer(3)]])
>>> P, L, U = B.LU(pivot='partial')
Traceback (most recent call last):
...
TypeError: cannot take absolute value of matrix entries, try 'pivot=nonzero'
>>> P, L, U = B.LU(pivot='nonzero')
>>> P
[1 0]
[0 1]
>>> L
[          1           0]
[y^2/(y + 1)           1]
>>> U
[  y + 1 y^2 + y]
[      0       0]
>>> L.base_ring()
Fraction Field of Univariate Polynomial Ring in y over Rational Field
>>> B == P*L*U
True

>>> # needs sage.rings.finite_rings
>>> F = FiniteField(Integer(5)**Integer(2), names=('a',)); (a,) = F._first_ngens(1)
>>> C = matrix(F, [[a + Integer(3), Integer(4)*a + Integer(4), Integer(2), Integer(4)*a + Integer(2)],
...                [Integer(3), Integer(2)*a + Integer(4), Integer(2)*a + Integer(4), Integer(2)*a + Integer(1)],
...                [Integer(3)*a + Integer(1), a + Integer(3), Integer(2)*a + Integer(4), Integer(4)*a + Integer(3)],
...                [a, Integer(3), Integer(3)*a + Integer(1), a]])
>>> P, L, U = C.LU(pivot='nonzero')
>>> P
[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
>>> L                                                                     # needs sage.combinat
[      1       0       0       0]
[3*a + 3       1       0       0]
[    2*a 4*a + 2       1       0]
[2*a + 3       2 2*a + 4       1]
>>> U                                                                     # needs sage.combinat
[  a + 3 4*a + 4       2 4*a + 2]
[      0   a + 1   a + 3 2*a + 4]
[      0       0       1 4*a + 2]
[      0       0       0       0]
>>> L.base_ring()                                                         # needs sage.combinat
Finite Field in a of size 5^2
>>> C == P*L*U                                                            # needs sage.combinat
True

With no pivoting strategy given (i.e. pivot=None) the routine will try to use partial pivoting, but then fall back to the nonzero strategy. For the nonsingular matrix below, we see evidence of pivoting when viewed over the rationals, and no pivoting over the integers mod 29.

sage: entries = [3, 20, 11, 7, 16, 28, 5, 15, 21, 23, 22, 18, 8, 23, 15, 2]
sage: A = matrix(Integers(29), 4, 4, entries)
sage: perm, _ = A.LU(format='compact'); perm
(0, 1, 2, 3)
sage: B = matrix(QQ, 4, 4, entries)
sage: perm, _ = B.LU(format='compact'); perm
(2, 0, 1, 3)
>>> from sage.all import *
>>> entries = [Integer(3), Integer(20), Integer(11), Integer(7), Integer(16), Integer(28), Integer(5), Integer(15), Integer(21), Integer(23), Integer(22), Integer(18), Integer(8), Integer(23), Integer(15), Integer(2)]
>>> A = matrix(Integers(Integer(29)), Integer(4), Integer(4), entries)
>>> perm, _ = A.LU(format='compact'); perm
(0, 1, 2, 3)
>>> B = matrix(QQ, Integer(4), Integer(4), entries)
>>> perm, _ = B.LU(format='compact'); perm
(2, 0, 1, 3)

The \(U\) matrix is only guaranteed to be upper-triangular. The rows are not necessarily linearly independent, nor are the pivots columns or rank in evidence.

sage: A = matrix(QQ, [[ 1, -4,  1,  0, -2,  1, 3,  3,  2],
....:                 [-1,  4,  0, -4,  0, -4, 5, -7, -7],
....:                 [ 0,  0,  1, -4, -1, -3, 6, -5, -6],
....:                 [-2,  8, -1, -4,  2, -4, 1, -8, -7],
....:                 [ 1, -4,  2, -4, -3,  2, 5,  6,  4]])
sage: P, L, U = A.LU()
sage: U
[   -2     8    -1    -4     2    -4     1    -8    -7]
[    0     0   1/2    -2    -1    -2   9/2    -3  -7/2]
[    0     0   3/2    -6    -2     0  11/2     2   1/2]
[    0     0     0     0  -1/3    -1   5/3  -5/3  -5/3]
[    0     0     0     0   1/3    -3   7/3 -19/3 -19/3]
sage: A.rref()
[ 1 -4  0  4  0  0 -1 -1 -1]
[ 0  0  1 -4  0  0  1  0 -1]
[ 0  0  0  0  1  0 -2 -1 -1]
[ 0  0  0  0  0  1 -1  2  2]
[ 0  0  0  0  0  0  0  0  0]
sage: A.pivots()
(0, 2, 4, 5)
>>> from sage.all import *
>>> A = matrix(QQ, [[ Integer(1), -Integer(4),  Integer(1),  Integer(0), -Integer(2),  Integer(1), Integer(3),  Integer(3),  Integer(2)],
...                 [-Integer(1),  Integer(4),  Integer(0), -Integer(4),  Integer(0), -Integer(4), Integer(5), -Integer(7), -Integer(7)],
...                 [ Integer(0),  Integer(0),  Integer(1), -Integer(4), -Integer(1), -Integer(3), Integer(6), -Integer(5), -Integer(6)],
...                 [-Integer(2),  Integer(8), -Integer(1), -Integer(4),  Integer(2), -Integer(4), Integer(1), -Integer(8), -Integer(7)],
...                 [ Integer(1), -Integer(4),  Integer(2), -Integer(4), -Integer(3),  Integer(2), Integer(5),  Integer(6),  Integer(4)]])
>>> P, L, U = A.LU()
>>> U
[   -2     8    -1    -4     2    -4     1    -8    -7]
[    0     0   1/2    -2    -1    -2   9/2    -3  -7/2]
[    0     0   3/2    -6    -2     0  11/2     2   1/2]
[    0     0     0     0  -1/3    -1   5/3  -5/3  -5/3]
[    0     0     0     0   1/3    -3   7/3 -19/3 -19/3]
>>> A.rref()
[ 1 -4  0  4  0  0 -1 -1 -1]
[ 0  0  1 -4  0  0  1  0 -1]
[ 0  0  0  0  1  0 -2 -1 -1]
[ 0  0  0  0  0  1 -1  2  2]
[ 0  0  0  0  0  0  0  0  0]
>>> A.pivots()
(0, 2, 4, 5)
QR(full=True)[source]

Return a factorization of self as a unitary matrix and an upper-triangular matrix.

INPUT:

  • full – (default: True) if True then the returned matrices have dimensions as described below. If False the R matrix has no zero rows and the columns of Q are a basis for the column space of self.

OUTPUT:

If self is an \(m\times n\) matrix and full=True then this method returns a pair of matrices: \(Q\) is an \(m\times m\) unitary matrix (meaning its inverse is its conjugate-transpose) and \(R\) is an \(m\times n\) upper-triangular matrix with nonnegative entries on the diagonal. For a matrix of full rank this factorization is unique (due to the restriction to positive entries on the diagonal).

If full=False then \(Q\) has \(m\) rows and the columns form an orthonormal basis for the column space of self. So, in particular, the conjugate-transpose of \(Q\) times \(Q\) will be an identity matrix. The matrix \(R\) will still be upper-triangular but will also have full rank, in particular it will lack the zero rows present in a full factorization of a rank-deficient matrix.

The results obtained when full=True are cached, hence \(Q\) and \(R\) are immutable matrices in this case.

Note

This is an exact computation, so limited to exact rings. Also the base ring needs to have a fraction field implemented in Sage and this field must contain square roots. One example is the field of algebraic numbers, QQbar, as used in the examples below. If you need numerical results, convert the base ring to the field of complex double numbers, CDF, which will use a faster routine that is careful about numerical subtleties.

ALGORITHM:

“Modified Gram-Schmidt,” Algorithm 8.1 of [TB1997].

EXAMPLES:

For a nonsingular matrix, the QR decomposition is unique.

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [[-2, 0, -4, -1, -1],
....:                    [-2, 1, -6, -3, -1],
....:                    [1, 1, 7, 4, 5],
....:                    [3, 0, 8, 3, 3],
....:                    [-1, 1, -6, -6, 5]])
sage: Q, R = A.QR()
sage: Q
[ -0.4588314677411235?  -0.1260506983326509?   0.3812120831224489?   -0.394573711338418?     -0.6874400625964?]
[ -0.4588314677411235?   0.4726901187474409? -0.05198346588033394?   0.7172941251646595?     -0.2209628772631?]
[  0.2294157338705618?   0.6617661662464172?   0.6619227988762521?  -0.1808720937375480?      0.1964114464561?]
[  0.6882472016116853?   0.1890760474989764?  -0.2044682991293135?   0.0966302966543065?     -0.6628886317894?]
[ -0.2294157338705618?   0.5357154679137663?   -0.609939332995919?   -0.536422031427112?      0.0245514308070?]
sage: R
[  4.358898943540674? -0.4588314677411235?   13.07669683062202?   6.194224814505168?   2.982404540317303?]
[                   0   1.670171752907625?  0.5987408170800917?  -1.292019657909672?   6.207996892883057?]
[                   0                    0   5.444401659866974?   5.468660610611130? -0.6827161852283857?]
[                   0                    0                    0   1.027626039419836?  -3.619300149686620?]
[                   0                    0                    0                    0   0.024551430807012?]
sage: Q.conjugate_transpose()*Q
[1.000000000000000?            0.?e-18            0.?e-17            0.?e-16            0.?e-13]
[           0.?e-18 1.000000000000000?            0.?e-17            0.?e-16            0.?e-13]
[           0.?e-17            0.?e-17 1.000000000000000?            0.?e-16            0.?e-13]
[           0.?e-16            0.?e-16            0.?e-16 1.000000000000000?            0.?e-13]
[           0.?e-13            0.?e-13            0.?e-13            0.?e-13   1.0000000000000?]
sage: Q * R == A
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [[-Integer(2), Integer(0), -Integer(4), -Integer(1), -Integer(1)],
...                    [-Integer(2), Integer(1), -Integer(6), -Integer(3), -Integer(1)],
...                    [Integer(1), Integer(1), Integer(7), Integer(4), Integer(5)],
...                    [Integer(3), Integer(0), Integer(8), Integer(3), Integer(3)],
...                    [-Integer(1), Integer(1), -Integer(6), -Integer(6), Integer(5)]])
>>> Q, R = A.QR()
>>> Q
[ -0.4588314677411235?  -0.1260506983326509?   0.3812120831224489?   -0.394573711338418?     -0.6874400625964?]
[ -0.4588314677411235?   0.4726901187474409? -0.05198346588033394?   0.7172941251646595?     -0.2209628772631?]
[  0.2294157338705618?   0.6617661662464172?   0.6619227988762521?  -0.1808720937375480?      0.1964114464561?]
[  0.6882472016116853?   0.1890760474989764?  -0.2044682991293135?   0.0966302966543065?     -0.6628886317894?]
[ -0.2294157338705618?   0.5357154679137663?   -0.609939332995919?   -0.536422031427112?      0.0245514308070?]
>>> R
[  4.358898943540674? -0.4588314677411235?   13.07669683062202?   6.194224814505168?   2.982404540317303?]
[                   0   1.670171752907625?  0.5987408170800917?  -1.292019657909672?   6.207996892883057?]
[                   0                    0   5.444401659866974?   5.468660610611130? -0.6827161852283857?]
[                   0                    0                    0   1.027626039419836?  -3.619300149686620?]
[                   0                    0                    0                    0   0.024551430807012?]
>>> Q.conjugate_transpose()*Q
[1.000000000000000?            0.?e-18            0.?e-17            0.?e-16            0.?e-13]
[           0.?e-18 1.000000000000000?            0.?e-17            0.?e-16            0.?e-13]
[           0.?e-17            0.?e-17 1.000000000000000?            0.?e-16            0.?e-13]
[           0.?e-16            0.?e-16            0.?e-16 1.000000000000000?            0.?e-13]
[           0.?e-13            0.?e-13            0.?e-13            0.?e-13   1.0000000000000?]
>>> Q * R == A
True

An example with complex numbers in QQbar, the field of algebraic numbers.

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [[-8, 4*I + 1, -I + 2, 2*I + 1],
....:                    [1, -2*I - 1, -I + 3, -I + 1],
....:                    [I + 7, 2*I + 1, -2*I + 7, -I + 1],
....:                    [I + 2, 0, I + 12, -1]])
sage: Q, R = A.QR()
sage: Q
[                          -0.7302967433402215?    0.2070566455055649? + 0.5383472783144687?*I    0.2463049809998642? - 0.0764456358723292?*I    0.2381617683194332? - 0.1036596032779695?*I]
[                           0.0912870929175277?   -0.2070566455055649? - 0.3778783780476559?*I    0.3786559533863033? - 0.1952221495524667?*I     0.701244450214469? - 0.3643711650986595?*I]
[   0.6390096504226938? + 0.0912870929175277?*I    0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I    0.3140171085506764? - 0.0825191718705412?*I]
[   0.1825741858350554? + 0.0912870929175277?*I  -0.03623491296347385? + 0.0724698259269477?*I   0.8632284069415110? + 0.06322839976356195?*I   -0.4499694867611521? - 0.0116119181208918?*I]
sage: R
[                          10.95445115010333?               0.?e-18 - 1.917028951268082?*I    5.385938482134133? - 2.190890230020665?*I  -0.2738612787525831? - 2.190890230020665?*I]
[                                           0               4.829596256417300? + 0.?e-18*I   -0.869637911123373? - 5.864879483945125?*I   0.993871898426712? - 0.3054085521207082?*I]
[                                           0                                            0               12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I]
[                                           0                                            0                                            0               1.942963944258992? + 0.?e-16*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-19*I            0.?e-18 + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-17 + 0.?e-17*I            0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-16*I]
sage: Q*R - A
[            0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[            0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [[-Integer(8), Integer(4)*I + Integer(1), -I + Integer(2), Integer(2)*I + Integer(1)],
...                    [Integer(1), -Integer(2)*I - Integer(1), -I + Integer(3), -I + Integer(1)],
...                    [I + Integer(7), Integer(2)*I + Integer(1), -Integer(2)*I + Integer(7), -I + Integer(1)],
...                    [I + Integer(2), Integer(0), I + Integer(12), -Integer(1)]])
>>> Q, R = A.QR()
>>> Q
[                          -0.7302967433402215?    0.2070566455055649? + 0.5383472783144687?*I    0.2463049809998642? - 0.0764456358723292?*I    0.2381617683194332? - 0.1036596032779695?*I]
[                           0.0912870929175277?   -0.2070566455055649? - 0.3778783780476559?*I    0.3786559533863033? - 0.1952221495524667?*I     0.701244450214469? - 0.3643711650986595?*I]
[   0.6390096504226938? + 0.0912870929175277?*I    0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I    0.3140171085506764? - 0.0825191718705412?*I]
[   0.1825741858350554? + 0.0912870929175277?*I  -0.03623491296347385? + 0.0724698259269477?*I   0.8632284069415110? + 0.06322839976356195?*I   -0.4499694867611521? - 0.0116119181208918?*I]
>>> R
[                          10.95445115010333?               0.?e-18 - 1.917028951268082?*I    5.385938482134133? - 2.190890230020665?*I  -0.2738612787525831? - 2.190890230020665?*I]
[                                           0               4.829596256417300? + 0.?e-18*I   -0.869637911123373? - 5.864879483945125?*I   0.993871898426712? - 0.3054085521207082?*I]
[                                           0                                            0               12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I]
[                                           0                                            0                                            0               1.942963944258992? + 0.?e-16*I]
>>> Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-19*I            0.?e-18 + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-17 + 0.?e-17*I            0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-16*I]
>>> Q*R - A
[            0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[            0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]

A rank-deficient rectangular matrix, with both values of the full keyword.

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [[2, -3, 3],
....:                    [-1, 1, -1],
....:                    [-1, 3, -3],
....:                    [-5, 1, -1]])
sage: Q, R = A.QR()
sage: Q
[  0.3592106040535498?  -0.5693261797050169?   0.7239227659930268?   0.1509015305256380?]
[ -0.1796053020267749?   0.1445907757980996?                     0   0.9730546968377341?]
[ -0.1796053020267749?   0.7048800320157352?    0.672213996993525?  -0.1378927778941174?]
[ -0.8980265101338745?  -0.3976246334447737?   0.1551263069985058? -0.10667177157846818?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
[                  0                   0                   0]
[                  0                   0                   0]
sage: Q.conjugate_transpose() * Q
[                 1            0.?e-18            0.?e-18            0.?e-18]
[           0.?e-18                  1            0.?e-18            0.?e-18]
[           0.?e-18            0.?e-18 1.000000000000000?            0.?e-18]
[           0.?e-18            0.?e-18            0.?e-18 1.000000000000000?]

sage: # needs sage.rings.number_field
sage: Q, R = A.QR(full=False)
sage: Q
[ 0.3592106040535498? -0.5693261797050169?]
[-0.1796053020267749?  0.1445907757980996?]
[-0.1796053020267749?  0.7048800320157352?]
[-0.8980265101338745? -0.3976246334447737?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
sage: Q.conjugate_transpose()*Q
[      1 0.?e-18]
[0.?e-18       1]
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [[Integer(2), -Integer(3), Integer(3)],
...                    [-Integer(1), Integer(1), -Integer(1)],
...                    [-Integer(1), Integer(3), -Integer(3)],
...                    [-Integer(5), Integer(1), -Integer(1)]])
>>> Q, R = A.QR()
>>> Q
[  0.3592106040535498?  -0.5693261797050169?   0.7239227659930268?   0.1509015305256380?]
[ -0.1796053020267749?   0.1445907757980996?                     0   0.9730546968377341?]
[ -0.1796053020267749?   0.7048800320157352?    0.672213996993525?  -0.1378927778941174?]
[ -0.8980265101338745?  -0.3976246334447737?   0.1551263069985058? -0.10667177157846818?]
>>> R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
[                  0                   0                   0]
[                  0                   0                   0]
>>> Q.conjugate_transpose() * Q
[                 1            0.?e-18            0.?e-18            0.?e-18]
[           0.?e-18                  1            0.?e-18            0.?e-18]
[           0.?e-18            0.?e-18 1.000000000000000?            0.?e-18]
[           0.?e-18            0.?e-18            0.?e-18 1.000000000000000?]

>>> # needs sage.rings.number_field
>>> Q, R = A.QR(full=False)
>>> Q
[ 0.3592106040535498? -0.5693261797050169?]
[-0.1796053020267749?  0.1445907757980996?]
[-0.1796053020267749?  0.7048800320157352?]
[-0.8980265101338745? -0.3976246334447737?]
>>> R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
>>> Q.conjugate_transpose()*Q
[      1 0.?e-18]
[0.?e-18       1]

Another rank-deficient rectangular matrix, with complex entries, as a reduced decomposition.

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [[-3*I - 3, I - 3, -12*I + 1, -2],
....:                    [-I - 1, -2, 5*I - 1, -I - 2],
....:                    [-4*I - 4, I - 5, -7*I, -I - 4]])
sage: Q, R = A.QR(full=False)
sage: Q
[ -0.4160251471689219? - 0.4160251471689219?*I   0.5370861555295747? + 0.1790287185098583?*I]
[ -0.1386750490563073? - 0.1386750490563073?*I  -0.7519206177414046? - 0.2506402059138015?*I]
[ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I]
sage: R
[                        7.211102550927979?  3.328201177351375? - 5.269651864139676?*I   7.904477796209515? + 8.45917799243475?*I  4.021576422632911? - 2.634825932069838?*I]
[                                         0                         1.074172311059150?  -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I]
sage: Q.conjugate_transpose()*Q
[1 0]
[0 1]
sage: Q*R - A
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [[-Integer(3)*I - Integer(3), I - Integer(3), -Integer(12)*I + Integer(1), -Integer(2)],
...                    [-I - Integer(1), -Integer(2), Integer(5)*I - Integer(1), -I - Integer(2)],
...                    [-Integer(4)*I - Integer(4), I - Integer(5), -Integer(7)*I, -I - Integer(4)]])
>>> Q, R = A.QR(full=False)
>>> Q
[ -0.4160251471689219? - 0.4160251471689219?*I   0.5370861555295747? + 0.1790287185098583?*I]
[ -0.1386750490563073? - 0.1386750490563073?*I  -0.7519206177414046? - 0.2506402059138015?*I]
[ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I]
>>> R
[                        7.211102550927979?  3.328201177351375? - 5.269651864139676?*I   7.904477796209515? + 8.45917799243475?*I  4.021576422632911? - 2.634825932069838?*I]
[                                         0                         1.074172311059150?  -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I]
>>> Q.conjugate_transpose()*Q
[1 0]
[0 1]
>>> Q*R - A
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]

Results of full decompositions are cached and thus returned immutable.

sage: # needs sage.rings.number_field
sage: A = random_matrix(QQbar, 2, 2)
sage: Q, R = A.QR()
sage: Q.is_mutable()
False
sage: R.is_mutable()
False
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = random_matrix(QQbar, Integer(2), Integer(2))
>>> Q, R = A.QR()
>>> Q.is_mutable()
False
>>> R.is_mutable()
False

Trivial cases return trivial results of the correct size, and we check \(Q\) itself in one case.

sage: # needs sage.rings.number_field
sage: A = zero_matrix(QQbar, 0, 10)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(0, 0)
sage: R.nrows(), R.ncols()
(0, 10)
sage: A = zero_matrix(QQbar, 3, 0)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(3, 3)
sage: R.nrows(), R.ncols()
(3, 0)
sage: Q
[1 0 0]
[0 1 0]
[0 0 1]
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = zero_matrix(QQbar, Integer(0), Integer(10))
>>> Q, R = A.QR()
>>> Q.nrows(), Q.ncols()
(0, 0)
>>> R.nrows(), R.ncols()
(0, 10)
>>> A = zero_matrix(QQbar, Integer(3), Integer(0))
>>> Q, R = A.QR()
>>> Q.nrows(), Q.ncols()
(3, 3)
>>> R.nrows(), R.ncols()
(3, 0)
>>> Q
[1 0 0]
[0 1 0]
[0 0 1]
T

Return the transpose of a matrix.

EXAMPLES:

sage: A = matrix(QQ, 5, range(25))
sage: A.T
[ 0  5 10 15 20]
[ 1  6 11 16 21]
[ 2  7 12 17 22]
[ 3  8 13 18 23]
[ 4  9 14 19 24]
>>> from sage.all import *
>>> A = matrix(QQ, Integer(5), range(Integer(25)))
>>> A.T
[ 0  5 10 15 20]
[ 1  6 11 16 21]
[ 2  7 12 17 22]
[ 3  8 13 18 23]
[ 4  9 14 19 24]
adjoint(*args, **kwds)[source]

Deprecated: Use adjugate() instead. See Issue #10501 for details.

adjoint_classical()[source]

Return the adjugate matrix of self (that is, the transpose of the matrix of cofactors).

Let \(M\) be an \(n \times n\)-matrix. The adjugate matrix of \(M\) is the \(n \times n\)-matrix \(N\) whose \((i, j)\)-th entry is \((-1)^{i + j} \det(M_{j, i})\), where \(M_{j,i}\) is the matrix \(M\) with its \(j\)-th row and \(i\)-th column removed. It is known to satisfy \(NM = MN = \det(M)I\).

EXAMPLES:

sage: M = Matrix(ZZ,2,2,[5,2,3,4]); M
[5 2]
[3 4]
sage: N = M.adjugate(); N
[ 4 -2]
[-3  5]
sage: M * N
[14  0]
[ 0 14]
sage: N * M
[14  0]
[ 0 14]
sage: M = Matrix(QQ, 2, 2, [5/3,2/56, 33/13,41/10]); M
[  5/3  1/28]
[33/13 41/10]
sage: N = M.adjugate(); N                                                   # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
sage: M * N                                                                 # needs sage.libs.pari
[7363/1092         0]
[        0 7363/1092]
>>> from sage.all import *
>>> M = Matrix(ZZ,Integer(2),Integer(2),[Integer(5),Integer(2),Integer(3),Integer(4)]); M
[5 2]
[3 4]
>>> N = M.adjugate(); N
[ 4 -2]
[-3  5]
>>> M * N
[14  0]
[ 0 14]
>>> N * M
[14  0]
[ 0 14]
>>> M = Matrix(QQ, Integer(2), Integer(2), [Integer(5)/Integer(3),Integer(2)/Integer(56), Integer(33)/Integer(13),Integer(41)/Integer(10)]); M
[  5/3  1/28]
[33/13 41/10]
>>> N = M.adjugate(); N                                                   # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
>>> M * N                                                                 # needs sage.libs.pari
[7363/1092         0]
[        0 7363/1092]

An alias is adjoint_classical(), which replaces the deprecated adjoint() method:

sage: M.adjoint()                                                           # needs sage.libs.pari
...: DeprecationWarning: adjoint is deprecated. Please use adjugate instead.
See https://github.com/sagemath/sage/issues/10501 for details.
[ 41/10  -1/28]
[-33/13    5/3]
sage: M.adjoint_classical()                                                 # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
>>> from sage.all import *
>>> M.adjoint()                                                           # needs sage.libs.pari
...: DeprecationWarning: adjoint is deprecated. Please use adjugate instead.
See https://github.com/sagemath/sage/issues/10501 for details.
[ 41/10  -1/28]
[-33/13    5/3]
>>> M.adjoint_classical()                                                 # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]

ALGORITHM:

Use PARI whenever the method self._adjugate is included to do so in an inheriting class. Otherwise, use a generic division-free algorithm that computes the adjugate matrix from the characteristic polynomial.

The result is cached.

adjugate()[source]

Return the adjugate matrix of self (that is, the transpose of the matrix of cofactors).

Let \(M\) be an \(n \times n\)-matrix. The adjugate matrix of \(M\) is the \(n \times n\)-matrix \(N\) whose \((i, j)\)-th entry is \((-1)^{i + j} \det(M_{j, i})\), where \(M_{j,i}\) is the matrix \(M\) with its \(j\)-th row and \(i\)-th column removed. It is known to satisfy \(NM = MN = \det(M)I\).

EXAMPLES:

sage: M = Matrix(ZZ,2,2,[5,2,3,4]); M
[5 2]
[3 4]
sage: N = M.adjugate(); N
[ 4 -2]
[-3  5]
sage: M * N
[14  0]
[ 0 14]
sage: N * M
[14  0]
[ 0 14]
sage: M = Matrix(QQ, 2, 2, [5/3,2/56, 33/13,41/10]); M
[  5/3  1/28]
[33/13 41/10]
sage: N = M.adjugate(); N                                                   # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
sage: M * N                                                                 # needs sage.libs.pari
[7363/1092         0]
[        0 7363/1092]
>>> from sage.all import *
>>> M = Matrix(ZZ,Integer(2),Integer(2),[Integer(5),Integer(2),Integer(3),Integer(4)]); M
[5 2]
[3 4]
>>> N = M.adjugate(); N
[ 4 -2]
[-3  5]
>>> M * N
[14  0]
[ 0 14]
>>> N * M
[14  0]
[ 0 14]
>>> M = Matrix(QQ, Integer(2), Integer(2), [Integer(5)/Integer(3),Integer(2)/Integer(56), Integer(33)/Integer(13),Integer(41)/Integer(10)]); M
[  5/3  1/28]
[33/13 41/10]
>>> N = M.adjugate(); N                                                   # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
>>> M * N                                                                 # needs sage.libs.pari
[7363/1092         0]
[        0 7363/1092]

An alias is adjoint_classical(), which replaces the deprecated adjoint() method:

sage: M.adjoint()                                                           # needs sage.libs.pari
...: DeprecationWarning: adjoint is deprecated. Please use adjugate instead.
See https://github.com/sagemath/sage/issues/10501 for details.
[ 41/10  -1/28]
[-33/13    5/3]
sage: M.adjoint_classical()                                                 # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]
>>> from sage.all import *
>>> M.adjoint()                                                           # needs sage.libs.pari
...: DeprecationWarning: adjoint is deprecated. Please use adjugate instead.
See https://github.com/sagemath/sage/issues/10501 for details.
[ 41/10  -1/28]
[-33/13    5/3]
>>> M.adjoint_classical()                                                 # needs sage.libs.pari
[ 41/10  -1/28]
[-33/13    5/3]

ALGORITHM:

Use PARI whenever the method self._adjugate is included to do so in an inheriting class. Otherwise, use a generic division-free algorithm that computes the adjugate matrix from the characteristic polynomial.

The result is cached.

apply_map(phi, R=None, sparse=None)[source]

Apply the given map phi (an arbitrary Python function or callable object) to this dense matrix. If R is not given, automatically determine the base ring of the resulting matrix.

INPUT:

  • sparse – boolean (default: False); ``True to make the output a sparse matrix

  • phi – arbitrary Python function or callable object

  • R – (optional) ring

OUTPUT: a matrix over R

EXAMPLES:

sage: m = matrix(ZZ, 3, 3, range(9))
sage: k.<a> = GF(9)                                                         # needs sage.rings.finite_rings
sage: f = lambda x: k(x)
sage: n = m.apply_map(f); n                                                 # needs sage.rings.finite_rings
[0 1 2]
[0 1 2]
[0 1 2]
sage: n.parent()                                                            # needs sage.rings.finite_rings
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field in a of size 3^2
>>> from sage.all import *
>>> m = matrix(ZZ, Integer(3), Integer(3), range(Integer(9)))
>>> k = GF(Integer(9), names=('a',)); (a,) = k._first_ngens(1)# needs sage.rings.finite_rings
>>> f = lambda x: k(x)
>>> n = m.apply_map(f); n                                                 # needs sage.rings.finite_rings
[0 1 2]
[0 1 2]
[0 1 2]
>>> n.parent()                                                            # needs sage.rings.finite_rings
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field in a of size 3^2

In this example, we explicitly specify the codomain.

sage: s = GF(3)
sage: f = lambda x: s(x)
sage: n = m.apply_map(f, k); n                                              # needs sage.rings.finite_rings
[0 1 2]
[0 1 2]
[0 1 2]
sage: n.parent()                                                            # needs sage.rings.finite_rings
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field in a of size 3^2
>>> from sage.all import *
>>> s = GF(Integer(3))
>>> f = lambda x: s(x)
>>> n = m.apply_map(f, k); n                                              # needs sage.rings.finite_rings
[0 1 2]
[0 1 2]
[0 1 2]
>>> n.parent()                                                            # needs sage.rings.finite_rings
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field in a of size 3^2

If self is subdivided, the result will be as well:

sage: m = matrix(2, 2, srange(4))
sage: m.subdivide(None, 1); m
[0|1]
[2|3]
sage: m.apply_map(lambda x: x*x)
[0|1]
[4|9]
>>> from sage.all import *
>>> m = matrix(Integer(2), Integer(2), srange(Integer(4)))
>>> m.subdivide(None, Integer(1)); m
[0|1]
[2|3]
>>> m.apply_map(lambda x: x*x)
[0|1]
[4|9]

If the matrix is sparse, the result will be as well:

sage: m = matrix(ZZ,100,100,sparse=True)
sage: m[18,32] = -6
sage: m[1,83] = 19
sage: n = m.apply_map(abs, R=ZZ)
sage: n.dict()
{(1, 83): 19, (18, 32): 6}
sage: n.is_sparse()
True
>>> from sage.all import *
>>> m = matrix(ZZ,Integer(100),Integer(100),sparse=True)
>>> m[Integer(18),Integer(32)] = -Integer(6)
>>> m[Integer(1),Integer(83)] = Integer(19)
>>> n = m.apply_map(abs, R=ZZ)
>>> n.dict()
{(1, 83): 19, (18, 32): 6}
>>> n.is_sparse()
True

If the map sends most of the matrix to zero, then it may be useful to get the result as a sparse matrix.

sage: m = matrix(ZZ, 3, 3, range(1, 10))
sage: n = m.apply_map(lambda x: 1//x, sparse=True); n
[1 0 0]
[0 0 0]
[0 0 0]
sage: n.parent()
Full MatrixSpace of 3 by 3 sparse matrices over Integer Ring
>>> from sage.all import *
>>> m = matrix(ZZ, Integer(3), Integer(3), range(Integer(1), Integer(10)))
>>> n = m.apply_map(lambda x: Integer(1)//x, sparse=True); n
[1 0 0]
[0 0 0]
[0 0 0]
>>> n.parent()
Full MatrixSpace of 3 by 3 sparse matrices over Integer Ring
apply_morphism(phi)[source]

Apply the morphism phi to the coefficients of this dense matrix.

The resulting matrix is over the codomain of phi.

INPUT:

  • phi – a morphism, so phi is callable and phi.domain() and phi.codomain() are defined. The codomain must be a ring

OUTPUT: a matrix over the codomain of phi

EXAMPLES:

sage: m = matrix(ZZ, 3, 3, range(9))
sage: phi = ZZ.hom(GF(5))
sage: m.apply_morphism(phi)
[0 1 2]
[3 4 0]
[1 2 3]
sage: parent(m.apply_morphism(phi))
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field of size 5
>>> from sage.all import *
>>> m = matrix(ZZ, Integer(3), Integer(3), range(Integer(9)))
>>> phi = ZZ.hom(GF(Integer(5)))
>>> m.apply_morphism(phi)
[0 1 2]
[3 4 0]
[1 2 3]
>>> parent(m.apply_morphism(phi))
Full MatrixSpace of 3 by 3 dense matrices
 over Finite Field of size 5

We apply a morphism to a matrix over a polynomial ring:

sage: R.<x,y> = QQ[]
sage: m = matrix(2, [x,x^2 + y, 2/3*y^2-x, x]); m
[          x     x^2 + y]
[2/3*y^2 - x           x]
sage: phi = R.hom([y,x])
sage: m.apply_morphism(phi)
[          y     y^2 + x]
[2/3*x^2 - y           y]
>>> from sage.all import *
>>> R = QQ['x, y']; (x, y,) = R._first_ngens(2)
>>> m = matrix(Integer(2), [x,x**Integer(2) + y, Integer(2)/Integer(3)*y**Integer(2)-x, x]); m
[          x     x^2 + y]
[2/3*y^2 - x           x]
>>> phi = R.hom([y,x])
>>> m.apply_morphism(phi)
[          y     y^2 + x]
[2/3*x^2 - y           y]
as_bipartite_graph()[source]

Construct a bipartite graph B representing the matrix uniquely.

Vertices are labeled 1 to nrows on the left and nrows + 1 to nrows + ncols on the right, representing rows and columns correspondingly. Each row is connected to each column with an edge weighted by the value of the corresponding matrix entry.

This graph is a helper for calculating automorphisms of a matrix under row and column permutations. See automorphisms_of_rows_and_columns().

OUTPUT: a bipartite graph

EXAMPLES:

sage: M = matrix(QQ, [[1/3, 7], [6, 1/4], [8, -5]])
sage: M
[1/3   7]
[  6 1/4]
[  8  -5]

sage: # needs sage.graphs
sage: B = M.as_bipartite_graph(); B
Bipartite graph on 5 vertices
sage: B.edges(sort=True)
[(1, 4, 1/3), (1, 5, 7), (2, 4, 6), (2, 5, 1/4), (3, 4, 8), (3, 5, -5)]
sage: len(B.left) == M.nrows()
True
sage: len(B.right) == M.ncols()
True
>>> from sage.all import *
>>> M = matrix(QQ, [[Integer(1)/Integer(3), Integer(7)], [Integer(6), Integer(1)/Integer(4)], [Integer(8), -Integer(5)]])
>>> M
[1/3   7]
[  6 1/4]
[  8  -5]

>>> # needs sage.graphs
>>> B = M.as_bipartite_graph(); B
Bipartite graph on 5 vertices
>>> B.edges(sort=True)
[(1, 4, 1/3), (1, 5, 7), (2, 4, 6), (2, 5, 1/4), (3, 4, 8), (3, 5, -5)]
>>> len(B.left) == M.nrows()
True
>>> len(B.right) == M.ncols()
True
as_sum_of_permutations()[source]

Return the current matrix as a sum of permutation matrices.

According to the Birkhoff-von Neumann Theorem, any bistochastic matrix can be written as a positive sum of permutation matrices, which also means that the polytope of bistochastic matrices is integer.

As a non-bistochastic matrix can obviously not be written as a sum of permutations, this theorem is an equivalence.

This function, given a bistochastic matrix, returns the corresponding decomposition.

See also

EXAMPLES:

We create a bistochastic matrix from a convex sum of permutations, then try to deduce the decomposition from the matrix

sage: L = []
sage: L.append((9, Permutation([4, 1, 3, 5, 2])))
sage: L.append((6, Permutation([5, 3, 4, 1, 2])))
sage: L.append((3, Permutation([3, 1, 4, 2, 5])))
sage: L.append((2, Permutation([1, 4, 2, 3, 5])))
sage: M = sum([c * p.to_matrix() for c, p in L])
sage: from sage.combinat.permutation import bistochastic_as_sum_of_permutations
sage: decomp = bistochastic_as_sum_of_permutations(M)                       # needs sage.combinat sage.graphs
sage: print(decomp)                                                         # needs sage.combinat sage.graphs
2*B[[1, 4, 2, 3, 5]] + 3*B[[3, 1, 4, 2, 5]] + 9*B[[4, 1, 3, 5, 2]] + 6*B[[5, 3, 4, 1, 2]]
>>> from sage.all import *
>>> L = []
>>> L.append((Integer(9), Permutation([Integer(4), Integer(1), Integer(3), Integer(5), Integer(2)])))
>>> L.append((Integer(6), Permutation([Integer(5), Integer(3), Integer(4), Integer(1), Integer(2)])))
>>> L.append((Integer(3), Permutation([Integer(3), Integer(1), Integer(4), Integer(2), Integer(5)])))
>>> L.append((Integer(2), Permutation([Integer(1), Integer(4), Integer(2), Integer(3), Integer(5)])))
>>> M = sum([c * p.to_matrix() for c, p in L])
>>> from sage.combinat.permutation import bistochastic_as_sum_of_permutations
>>> decomp = bistochastic_as_sum_of_permutations(M)                       # needs sage.combinat sage.graphs
>>> print(decomp)                                                         # needs sage.combinat sage.graphs
2*B[[1, 4, 2, 3, 5]] + 3*B[[3, 1, 4, 2, 5]] + 9*B[[4, 1, 3, 5, 2]] + 6*B[[5, 3, 4, 1, 2]]

An exception is raised when the matrix is not bistochastic:

sage: M = Matrix([[2,3],[2,2]])
sage: decomp = bistochastic_as_sum_of_permutations(M)                       # needs sage.graphs
Traceback (most recent call last):
...
ValueError: The matrix is not bistochastic
>>> from sage.all import *
>>> M = Matrix([[Integer(2),Integer(3)],[Integer(2),Integer(2)]])
>>> decomp = bistochastic_as_sum_of_permutations(M)                       # needs sage.graphs
Traceback (most recent call last):
...
ValueError: The matrix is not bistochastic
automorphisms_of_rows_and_columns()[source]

Return the automorphisms of self under permutations of rows and columns as a list of pairs of PermutationGroupElement objects.

EXAMPLES:

sage: # needs sage.graphs sage.groups
sage: M = matrix(ZZ,[[1,0],[1,0],[0,1]]); M
[1 0]
[1 0]
[0 1]
sage: A = M.automorphisms_of_rows_and_columns(); A
[((), ()), ((1,2), ())]
sage: M = matrix(ZZ, [[1,1,1,1],[1,1,1,1]])
sage: A = M.automorphisms_of_rows_and_columns()
sage: len(A)
48
>>> from sage.all import *
>>> # needs sage.graphs sage.groups
>>> M = matrix(ZZ,[[Integer(1),Integer(0)],[Integer(1),Integer(0)],[Integer(0),Integer(1)]]); M
[1 0]
[1 0]
[0 1]
>>> A = M.automorphisms_of_rows_and_columns(); A
[((), ()), ((1,2), ())]
>>> M = matrix(ZZ, [[Integer(1),Integer(1),Integer(1),Integer(1)],[Integer(1),Integer(1),Integer(1),Integer(1)]])
>>> A = M.automorphisms_of_rows_and_columns()
>>> len(A)
48

One can now apply these automorphisms to M to show that it leaves it invariant:

sage: all(M.with_permuted_rows_and_columns(*i) == M for i in A)             # needs sage.graphs sage.groups
True
>>> from sage.all import *
>>> all(M.with_permuted_rows_and_columns(*i) == M for i in A)             # needs sage.graphs sage.groups
True

Check that Issue #25426 is fixed:

sage: j = matrix([(3, 2, 1, 0, 0),
....:             (2, 2, 0, 1, 0),
....:             (1, 0, 3, 0, 2),
....:             (0, 1, 0, 2, 1),
....:             (0, 0, 2, 1, 2)])
sage: j.automorphisms_of_rows_and_columns()                                 # needs sage.graphs sage.groups
[((), ()), ((1,3)(2,5), (1,3)(2,5))]
>>> from sage.all import *
>>> j = matrix([(Integer(3), Integer(2), Integer(1), Integer(0), Integer(0)),
...             (Integer(2), Integer(2), Integer(0), Integer(1), Integer(0)),
...             (Integer(1), Integer(0), Integer(3), Integer(0), Integer(2)),
...             (Integer(0), Integer(1), Integer(0), Integer(2), Integer(1)),
...             (Integer(0), Integer(0), Integer(2), Integer(1), Integer(2))])
>>> j.automorphisms_of_rows_and_columns()                                 # needs sage.graphs sage.groups
[((), ()), ((1,3)(2,5), (1,3)(2,5))]
block_ldlt(classical=False)[source]

Compute a block-\(LDL^{T}\) factorization of a Hermitian matrix.

The standard \(LDL^{T}\) factorization of a positive-definite matrix \(A\) factors it as \(A = LDL^{T}\) where \(L\) is unit-lower-triangular and \(D\) is diagonal. If one allows row/column swaps via a permutation matrix \(P\), then this factorization can be extended to many positive-semidefinite matrices \(A\) via the factorization \(P^{T}AP = LDL^{T}\) that places the zeros at the bottom of \(D\) to avoid division by zero. These factorizations extend easily to complex Hermitian matrices when one replaces the transpose by the conjugate-transpose.

However, we can go one step further. If, in addition, we allow \(D\) to potentially contain \(2 \times 2\) blocks on its diagonal, then every real or complex Hermitian matrix \(A\) can be factored as \(A = PLDL^{*}P^{T}\). When the row/column swaps are made intelligently, this process is numerically stable over inexact rings like RDF. Bunch and Kaufman describe such a “pivot” scheme that is suitable for the solution of Hermitian systems, and that is how we choose our row and column swaps.

INPUT:

  • classical – boolean (default: False); whether or not to attempt a classical non-block \(LDL^{T}\) factorization with no row/column swaps.

Warning

Not all matrices have a classical \(LDL^{T}\) factorization. Set classical=True at your own risk, preferably after verifying that your matrix is positive-definite and (over inexact rings) not ill-conditioned.

OUTPUT:

If the input matrix is not Hermitian, the output from this function is undefined. Otherwise, we return a triple \((P,L,D)\) such that \(A = PLDL^{*}P^{T}\) and

  • \(P\) is a permutation matrix,

  • \(L\) is unit lower-triangular,

  • \(D\) is a block-diagonal matrix whose blocks are of size one or two.

With classical=True, the permutation matrix \(P\) is always an identity matrix and the diagonal blocks are always one-by-one. A ValueError is raised if the matrix has no classical \(LDL^{T}\) factorization.

ALGORITHM:

We essentially follow “Algorithm A” in the paper by Bunch and Kaufman [BK1977] that describes the stable pivoting strategy. The same scheme is described by Higham [Hig2002].

REFERENCES:

EXAMPLES:

This three-by-three real symmetric matrix has one positive, one negative, and one zero eigenvalue – so it is not any flavor of (semi)definite, yet we can still factor it:

sage: A =  matrix(QQ, [[0, 1, 0],
....:                  [1, 1, 2],
....:                  [0, 2, 0]])
sage: P,L,D = A.block_ldlt()
sage: P
[0 0 1]
[1 0 0]
[0 1 0]
sage: L
[  1   0   0]
[  2   1   0]
[  1 1/2   1]
sage: D
[ 1| 0| 0]
[--+--+--]
[ 0|-4| 0]
[--+--+--]
[ 0| 0| 0]
sage: P.transpose()*A*P == L*D*L.transpose()
True
>>> from sage.all import *
>>> A =  matrix(QQ, [[Integer(0), Integer(1), Integer(0)],
...                  [Integer(1), Integer(1), Integer(2)],
...                  [Integer(0), Integer(2), Integer(0)]])
>>> P,L,D = A.block_ldlt()
>>> P
[0 0 1]
[1 0 0]
[0 1 0]
>>> L
[  1   0   0]
[  2   1   0]
[  1 1/2   1]
>>> D
[ 1| 0| 0]
[--+--+--]
[ 0|-4| 0]
[--+--+--]
[ 0| 0| 0]
>>> P.transpose()*A*P == L*D*L.transpose()
True

This two-by-two matrix has no classical factorization, but it constitutes its own block-factorization:

sage: A = matrix(QQ, [ [0,1],
....:                  [1,0] ])
sage: A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
sage: A.block_ldlt()
(
[1 0]  [1 0]  [0 1]
[0 1], [0 1], [1 0]
)
>>> from sage.all import *
>>> A = matrix(QQ, [ [Integer(0),Integer(1)],
...                  [Integer(1),Integer(0)] ])
>>> A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
>>> A.block_ldlt()
(
[1 0]  [1 0]  [0 1]
[0 1], [0 1], [1 0]
)

The same is true of the following complex Hermitian matrix:

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [ [ 0,I],
....:                     [-I,0] ])
sage: A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
sage: A.block_ldlt()
(
[1 0]  [1 0]  [ 0  I]
[0 1], [0 1], [-I  0]
)
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [ [ Integer(0),I],
...                     [-I,Integer(0)] ])
>>> A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
>>> A.block_ldlt()
(
[1 0]  [1 0]  [ 0  I]
[0 1], [0 1], [-I  0]
)

Complete diagonal pivoting could cause problems for the following matrix, since the diagonal entries are small compared to the off-diagonals that must be zeroed; however, the block algorithm refuses to factor it:

sage: A = matrix(RDF, 2, 2, [ [1e-10, 1    ],
....:                         [1    , 2e-10] ])
sage: _,L,D = A.block_ldlt(classical=True)
sage: L*D*L.T
[1e-10   1.0]
[  1.0   0.0]
sage: A.block_ldlt()                                                        # needs scipy
(
[1.0 0.0]  [1.0 0.0]  [1e-10   1.0]
[0.0 1.0], [0.0 1.0], [  1.0 2e-10]
)
>>> from sage.all import *
>>> A = matrix(RDF, Integer(2), Integer(2), [ [RealNumber('1e-10'), Integer(1)    ],
...                         [Integer(1)    , RealNumber('2e-10')] ])
>>> _,L,D = A.block_ldlt(classical=True)
>>> L*D*L.T
[1e-10   1.0]
[  1.0   0.0]
>>> A.block_ldlt()                                                        # needs scipy
(
[1.0 0.0]  [1.0 0.0]  [1e-10   1.0]
[0.0 1.0], [0.0 1.0], [  1.0 2e-10]
)

The factorization over an inexact ring is necessarily inexact, but \(P^{T}AP\) will ideally be close to \(LDL^{*}\) in the metric induced by the norm:

sage: # needs scipy sage.rings.complex_double sage.symbolic
sage: A = matrix(CDF, 2, 2, [ [-1.1933, -0.3185 - 1.3553*I],
....:                         [-0.3185 + 1.3553*I, 1.5729 ] ])
sage: P,L,D = A.block_ldlt()
sage: P.T*A*P == L*D*L.H
False
sage: (P.T*A*P - L*D*L.H).norm() < 1e-10
True
>>> from sage.all import *
>>> # needs scipy sage.rings.complex_double sage.symbolic
>>> A = matrix(CDF, Integer(2), Integer(2), [ [-RealNumber('1.1933'), -RealNumber('0.3185') - RealNumber('1.3553')*I],
...                         [-RealNumber('0.3185') + RealNumber('1.3553')*I, RealNumber('1.5729') ] ])
>>> P,L,D = A.block_ldlt()
>>> P.T*A*P == L*D*L.H
False
>>> (P.T*A*P - L*D*L.H).norm() < RealNumber('1e-10')
True

This matrix has a singular three-by-three leading principal submatrix, and therefore has no classical factorization:

sage: A = matrix(QQ, [[21, 15, 12, -2],
....:                 [15, 12,  9,  6],
....:                 [12,  9,  7,  3],
....:                 [-2,  6,  3,  8]])
sage: A[0:3,0:3].det() == 0
True
sage: A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
sage: A.block_ldlt()
(
[1 0 0 0]  [     1      0      0      0]
[0 0 1 0]  [ -2/21      1      0      0]
[0 0 0 1]  [   5/7  39/41      1      0]
[0 1 0 0], [   4/7 87/164  48/79      1],

[     21|      0|      0|      0]
[-------+-------+-------+-------]
[      0| 164/21|      0|      0]
[-------+-------+-------+-------]
[      0|      0|-237/41|      0]
[-------+-------+-------+-------]
[      0|      0|      0| 25/316]
)
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(21), Integer(15), Integer(12), -Integer(2)],
...                 [Integer(15), Integer(12),  Integer(9),  Integer(6)],
...                 [Integer(12),  Integer(9),  Integer(7),  Integer(3)],
...                 [-Integer(2),  Integer(6),  Integer(3),  Integer(8)]])
>>> A[Integer(0):Integer(3),Integer(0):Integer(3)].det() == Integer(0)
True
>>> A.block_ldlt(classical=True)
Traceback (most recent call last):
...
ValueError: matrix has no classical LDL^T factorization
>>> A.block_ldlt()
(
[1 0 0 0]  [     1      0      0      0]
[0 0 1 0]  [ -2/21      1      0      0]
[0 0 0 1]  [   5/7  39/41      1      0]
[0 1 0 0], [   4/7 87/164  48/79      1],
<BLANKLINE>
[     21|      0|      0|      0]
[-------+-------+-------+-------]
[      0| 164/21|      0|      0]
[-------+-------+-------+-------]
[      0|      0|-237/41|      0]
[-------+-------+-------+-------]
[      0|      0|      0| 25/316]
)

An indefinite symmetric matrix that happens to have a classical factorization:

sage: A = matrix(QQ, [[ 3,  -6,   9,   6,  -9],
....:                 [-6,  11, -16, -11,  17],
....:                 [ 9, -16,  28,  16, -40],
....:                 [ 6, -11,  16,   9, -19],
....:                 [-9,  17, -40, -19,  68]])
sage: A.block_ldlt(classical=True)[1:]
(
                  [ 3| 0| 0| 0| 0]
                  [--+--+--+--+--]
                  [ 0|-1| 0| 0| 0]
                  [--+--+--+--+--]
[ 1  0  0  0  0]  [ 0| 0| 5| 0| 0]
[-2  1  0  0  0]  [--+--+--+--+--]
[ 3 -2  1  0  0]  [ 0| 0| 0|-2| 0]
[ 2 -1  0  1  0]  [--+--+--+--+--]
[-3  1 -3  1  1], [ 0| 0| 0| 0|-1]
)
>>> from sage.all import *
>>> A = matrix(QQ, [[ Integer(3),  -Integer(6),   Integer(9),   Integer(6),  -Integer(9)],
...                 [-Integer(6),  Integer(11), -Integer(16), -Integer(11),  Integer(17)],
...                 [ Integer(9), -Integer(16),  Integer(28),  Integer(16), -Integer(40)],
...                 [ Integer(6), -Integer(11),  Integer(16),   Integer(9), -Integer(19)],
...                 [-Integer(9),  Integer(17), -Integer(40), -Integer(19),  Integer(68)]])
>>> A.block_ldlt(classical=True)[Integer(1):]
(
                  [ 3| 0| 0| 0| 0]
                  [--+--+--+--+--]
                  [ 0|-1| 0| 0| 0]
                  [--+--+--+--+--]
[ 1  0  0  0  0]  [ 0| 0| 5| 0| 0]
[-2  1  0  0  0]  [--+--+--+--+--]
[ 3 -2  1  0  0]  [ 0| 0| 0|-2| 0]
[ 2 -1  0  1  0]  [--+--+--+--+--]
[-3  1 -3  1  1], [ 0| 0| 0| 0|-1]
)

An indefinite Hermitian matrix that happens to have a classical factorization:

sage: F.<I> = QuadraticField(-1)                                            # needs sage.rings.number_field
sage: A = matrix(F, [[      2, 4 - 2*I, 2 + 2*I],                           # needs sage.rings.number_field
....:                [4 + 2*I,       8,    10*I],
....:                [2 - 2*I,   -10*I,      -3]])
sage: A.block_ldlt(classical=True)[1:]                                      # needs sage.rings.number_field
(
                           [ 2| 0| 0]
                           [--+--+--]
[      1       0       0]  [ 0|-2| 0]
[  I + 2       1       0]  [--+--+--]
[ -I + 1 2*I + 1       1], [ 0| 0| 3]
)
>>> from sage.all import *
>>> F = QuadraticField(-Integer(1), names=('I',)); (I,) = F._first_ngens(1)# needs sage.rings.number_field
>>> A = matrix(F, [[      Integer(2), Integer(4) - Integer(2)*I, Integer(2) + Integer(2)*I],                           # needs sage.rings.number_field
...                [Integer(4) + Integer(2)*I,       Integer(8),    Integer(10)*I],
...                [Integer(2) - Integer(2)*I,   -Integer(10)*I,      -Integer(3)]])
>>> A.block_ldlt(classical=True)[Integer(1):]                                      # needs sage.rings.number_field
(
                           [ 2| 0| 0]
                           [--+--+--]
[      1       0       0]  [ 0|-2| 0]
[  I + 2       1       0]  [--+--+--]
[ -I + 1 2*I + 1       1], [ 0| 0| 3]
)
characteristic_polynomial(*args, **kwds)[source]

Synonym for self.charpoly(…).

EXAMPLES:

sage: a = matrix(QQ, 2,2, [1,2,3,4]); a
[1 2]
[3 4]
sage: a.characteristic_polynomial('T')                                      # needs sage.libs.pari
T^2 - 5*T - 2
>>> from sage.all import *
>>> a = matrix(QQ, Integer(2),Integer(2), [Integer(1),Integer(2),Integer(3),Integer(4)]); a
[1 2]
[3 4]
>>> a.characteristic_polynomial('T')                                      # needs sage.libs.pari
T^2 - 5*T - 2
charpoly(var='x', algorithm=None)[source]

Return the characteristic polynomial of self, as a polynomial over the base ring.

ALGORITHM:

If the base ring has a method \(_matrix_charpoly\), we use it.

In the generic case of matrices over a ring (commutative and with unity), there is a division-free algorithm, which can be accessed using 'df', with complexity \(O(n^4)\). Alternatively, by specifying 'hessenberg', this method computes the Hessenberg form of the matrix and then reads off the characteristic polynomial. Moreover, for matrices over number fields, this method can use PARI’s charpoly implementation instead.

The method’s logic is as follows: If no algorithm is specified, first check if the base ring is a number field (and then use PARI), otherwise check if the base ring is the ring of integers modulo n (in which case compute the characteristic polynomial of a lift of the matrix to the integers, and then coerce back to the base), next check if the base ring is an exact field (and then use the Hessenberg form), or otherwise, use the generic division-free algorithm. If an algorithm is specified explicitly, if algorithm == "hessenberg", use the Hessenberg form, or otherwise use the generic division-free algorithm.

The result is cached.

INPUT:

  • var – a variable name (default: 'x')

  • algorithm – string; one of
    • 'df' – generic \(O(n^4)\) division-free algorithm

    • 'hessenberg' – use the Hessenberg form of the matrix

EXAMPLES:

First a matrix over \(\ZZ\):

sage: A = MatrixSpace(ZZ,2)( [1,2,  3,4] )
sage: f = A.charpoly('x')
sage: f
x^2 - 5*x - 2
sage: f.parent()
Univariate Polynomial Ring in x over Integer Ring
sage: f(A)
[0 0]
[0 0]
>>> from sage.all import *
>>> A = MatrixSpace(ZZ,Integer(2))( [Integer(1),Integer(2),  Integer(3),Integer(4)] )
>>> f = A.charpoly('x')
>>> f
x^2 - 5*x - 2
>>> f.parent()
Univariate Polynomial Ring in x over Integer Ring
>>> f(A)
[0 0]
[0 0]

An example over \(\QQ\):

sage: A = MatrixSpace(QQ, 3)(range(9))
sage: A.charpoly('x')                                                       # needs sage.libs.pari
x^3 - 12*x^2 - 18*x
sage: A.trace()
12
sage: A.determinant()
0
>>> from sage.all import *
>>> A = MatrixSpace(QQ, Integer(3))(range(Integer(9)))
>>> A.charpoly('x')                                                       # needs sage.libs.pari
x^3 - 12*x^2 - 18*x
>>> A.trace()
12
>>> A.determinant()
0

We compute the characteristic polynomial of a matrix over the polynomial ring \(\ZZ[a]\):

sage: R.<a> = PolynomialRing(ZZ)
sage: M = MatrixSpace(R, 2)([a,1,  a,a+1]); M
[    a     1]
[    a a + 1]
sage: f = M.charpoly('x'); f
x^2 + (-2*a - 1)*x + a^2
sage: f.parent()
Univariate Polynomial Ring in x
 over Univariate Polynomial Ring in a over Integer Ring
sage: M.trace()
2*a + 1
sage: M.determinant()
a^2
>>> from sage.all import *
>>> R = PolynomialRing(ZZ, names=('a',)); (a,) = R._first_ngens(1)
>>> M = MatrixSpace(R, Integer(2))([a,Integer(1),  a,a+Integer(1)]); M
[    a     1]
[    a a + 1]
>>> f = M.charpoly('x'); f
x^2 + (-2*a - 1)*x + a^2
>>> f.parent()
Univariate Polynomial Ring in x
 over Univariate Polynomial Ring in a over Integer Ring
>>> M.trace()
2*a + 1
>>> M.determinant()
a^2

We compute the characteristic polynomial of a matrix over the multi-variate polynomial ring \(\ZZ[x,y]\):

sage: R.<x,y> = PolynomialRing(ZZ,2)
sage: A = MatrixSpace(R, 2)([x, y, x^2, y^2])
sage: f = A.charpoly('x'); f
x^2 + (-y^2 - x)*x - x^2*y + x*y^2
>>> from sage.all import *
>>> R = PolynomialRing(ZZ,Integer(2), names=('x', 'y',)); (x, y,) = R._first_ngens(2)
>>> A = MatrixSpace(R, Integer(2))([x, y, x**Integer(2), y**Integer(2)])
>>> f = A.charpoly('x'); f
x^2 + (-y^2 - x)*x - x^2*y + x*y^2

It’s a little difficult to distinguish the variables. To fix this, we temporarily view the indeterminate as \(Z\):

sage: with localvars(f.parent(), 'Z'): print(f)
Z^2 + (-y^2 - x)*Z - x^2*y + x*y^2
>>> from sage.all import *
>>> with localvars(f.parent(), 'Z'): print(f)
Z^2 + (-y^2 - x)*Z - x^2*y + x*y^2

We could also compute f in terms of Z from the start:

sage: A.charpoly('Z')
Z^2 + (-y^2 - x)*Z - x^2*y + x*y^2
>>> from sage.all import *
>>> A.charpoly('Z')
Z^2 + (-y^2 - x)*Z - x^2*y + x*y^2

Here is an example over a number field:

sage: # needs sage.rings.number_field
sage: x = QQ['x'].gen()
sage: K.<a> = NumberField(x^2 - 2)
sage: m = matrix(K, [[a-1, 2], [a, a+1]])
sage: m.charpoly('Z')
Z^2 - 2*a*Z - 2*a + 1
sage: m.charpoly('a')(m) == 0
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> x = QQ['x'].gen()
>>> K = NumberField(x**Integer(2) - Integer(2), names=('a',)); (a,) = K._first_ngens(1)
>>> m = matrix(K, [[a-Integer(1), Integer(2)], [a, a+Integer(1)]])
>>> m.charpoly('Z')
Z^2 - 2*a*Z - 2*a + 1
>>> m.charpoly('a')(m) == Integer(0)
True

Over integers modulo \(n\) with composite \(n\):

sage: A = Mat(Integers(6), 3, 3)(range(9))
sage: A.charpoly()
x^3
>>> from sage.all import *
>>> A = Mat(Integers(Integer(6)), Integer(3), Integer(3))(range(Integer(9)))
>>> A.charpoly()
x^3

Here is an example over a general commutative ring, that is to say, as of version 4.0.2, Sage does not even positively determine that S in the following example is an integral domain. But the computation of the characteristic polynomial succeeds as follows:

sage: # needs sage.libs.singular
sage: R.<a,b> = QQ[]
sage: S.<x,y> = R.quo((b^3))
sage: A = matrix(S, [[x*y^2, 2*x], [2, x^10*y]]); A
[ x*y^2    2*x]
[     2 x^10*y]
sage: A.charpoly('T')
T^2 + (-x^10*y - x*y^2)*T - 4*x
>>> from sage.all import *
>>> # needs sage.libs.singular
>>> R = QQ['a, b']; (a, b,) = R._first_ngens(2)
>>> S = R.quo((b**Integer(3)), names=('x', 'y',)); (x, y,) = S._first_ngens(2)
>>> A = matrix(S, [[x*y**Integer(2), Integer(2)*x], [Integer(2), x**Integer(10)*y]]); A
[ x*y^2    2*x]
[     2 x^10*y]
>>> A.charpoly('T')
T^2 + (-x^10*y - x*y^2)*T - 4*x
cholesky()[source]

Return the Cholesky decomposition of a Hermitian matrix.

Applies to a positive-definite matrix. Generally, the base ring for the entries of the matrix needs to be a subfield of the algebraic numbers (QQbar). Examples include the rational numbers (QQ), some number fields, and real algebraic numbers and the algebraic numbers themselves. Symbolic matrices can also occasionally be factored.

OUTPUT:

For a matrix \(A\) the routine returns a lower triangular matrix \(L\) such that,

\[A = LL^\ast\]

where \(L^\ast\) is the conjugate-transpose. If the matrix is not positive-definite (for example, if it is not Hermitian) then a ValueError results.

If possible, the output matrix will be over the fraction field of the base ring of the input matrix. If that fraction field is missing the requisite square roots but if no imaginaries are encountered, then the algebraic-reals will be used. Otherwise, the algebraic closure of the fraction field (typically QQbar) will be used.

ALGORITHM:

First we ensure that the matrix \(A\) is_hermitian(). Afterwards, we attempt to compute a classical block_ldlt() factorization, \(A = LDL^{*}\), of the matrix. If that fails, then the matrix was not positive-definite and an error is raised. Otherwise we take the entrywise square-root \(\sqrt{D}\) of the diagonal matrix \(D\) (whose entries are the positive eigenvalues of the original matrix) to obtain the Cholesky factorization \(A = \left(L\sqrt{D}\right)\left(L\sqrt{D}\right)^{*}\). If the necessary square roots cannot be taken in the fraction field of original base ring, then we move to either its algebraic closure or the algebraic reals, depending on whether or not imaginary numbers are required.

EXAMPLES:

This simple example has a result with entries that remain in the field of rational numbers:

sage: A = matrix(QQ, [[ 4, -2,  4,  2],
....:                 [-2, 10, -2, -7],
....:                 [ 4, -2,  8,  4],
....:                 [ 2, -7,  4,  7]])
sage: A.is_symmetric()
True
sage: L = A.cholesky(); L
[ 2  0  0  0]
[-1  3  0  0]
[ 2  0  2  0]
[ 1 -2  1  1]
sage: L.parent()
Full MatrixSpace of 4 by 4 dense matrices over Rational Field
sage: L*L.transpose() == A
True
>>> from sage.all import *
>>> A = matrix(QQ, [[ Integer(4), -Integer(2),  Integer(4),  Integer(2)],
...                 [-Integer(2), Integer(10), -Integer(2), -Integer(7)],
...                 [ Integer(4), -Integer(2),  Integer(8),  Integer(4)],
...                 [ Integer(2), -Integer(7),  Integer(4),  Integer(7)]])
>>> A.is_symmetric()
True
>>> L = A.cholesky(); L
[ 2  0  0  0]
[-1  3  0  0]
[ 2  0  2  0]
[ 1 -2  1  1]
>>> L.parent()
Full MatrixSpace of 4 by 4 dense matrices over Rational Field
>>> L*L.transpose() == A
True

This seemingly simple example requires first moving to the rational numbers for field operations, and then square roots necessitate that the result has entries in the field of algebraic numbers:

sage: A = matrix(ZZ, [[ 78, -30, -37,  -2],
....:                 [-30, 102, 179, -18],
....:                 [-37, 179, 326, -38],
....:                 [ -2, -18, -38,  15]])
sage: A.is_symmetric()
True
sage: L = A.cholesky(); L                                                   # needs sage.rings.number_field
[   8.83176086632785?                    0                    0                    0]
[ -3.396831102433787?    9.51112708681461?                    0                    0]
[ -4.189425026335004?   17.32383862241232?   2.886751345948129?                    0]
[-0.2264554068289192?  -1.973397116652010?  -1.649572197684645?   2.886751345948129?]
sage: L.parent()                                                            # needs sage.rings.number_field
Full MatrixSpace of 4 by 4 dense matrices over Algebraic Real Field
sage: L*L.transpose() == A                                                  # needs sage.rings.number_field
True
>>> from sage.all import *
>>> A = matrix(ZZ, [[ Integer(78), -Integer(30), -Integer(37),  -Integer(2)],
...                 [-Integer(30), Integer(102), Integer(179), -Integer(18)],
...                 [-Integer(37), Integer(179), Integer(326), -Integer(38)],
...                 [ -Integer(2), -Integer(18), -Integer(38),  Integer(15)]])
>>> A.is_symmetric()
True
>>> L = A.cholesky(); L                                                   # needs sage.rings.number_field
[   8.83176086632785?                    0                    0                    0]
[ -3.396831102433787?    9.51112708681461?                    0                    0]
[ -4.189425026335004?   17.32383862241232?   2.886751345948129?                    0]
[-0.2264554068289192?  -1.973397116652010?  -1.649572197684645?   2.886751345948129?]
>>> L.parent()                                                            # needs sage.rings.number_field
Full MatrixSpace of 4 by 4 dense matrices over Algebraic Real Field
>>> L*L.transpose() == A                                                  # needs sage.rings.number_field
True

Some subfields of the complex numbers, such as this number field of complex numbers with rational real and imaginary parts, allow for this computation:

sage: # needs sage.rings.number_field
sage: C.<I> = QuadraticField(-1)
sage: A = matrix(C, [[        23,  17*I + 3,  24*I + 25,     21*I],
....:                [ -17*I + 3,        38, -69*I + 89, 7*I + 15],
....:                [-24*I + 25, 69*I + 89,        976, 24*I + 6],
....:                [     -21*I, -7*I + 15,  -24*I + 6,       28]])
sage: A.is_hermitian()
True
sage: L = A.cholesky(); L
[                4.79...?                         0                       0        0]
[   0.62...? - 3.54...?*I                  5.00...?                       0        0]
[   5.21...? - 5.00...?*I   13.58...? + 10.72...?*I               24.98...?        0]
[             -4.37...?*I   -0.10...? -  0.85...?*I  -0.21...? + 0.37...?*I 2.81...?]
sage: L.parent()
Full MatrixSpace of 4 by 4 dense matrices over Algebraic Field
sage: (L*L.conjugate_transpose() - A.change_ring(QQbar)).norm() < 10^-10
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> C = QuadraticField(-Integer(1), names=('I',)); (I,) = C._first_ngens(1)
>>> A = matrix(C, [[        Integer(23),  Integer(17)*I + Integer(3),  Integer(24)*I + Integer(25),     Integer(21)*I],
...                [ -Integer(17)*I + Integer(3),        Integer(38), -Integer(69)*I + Integer(89), Integer(7)*I + Integer(15)],
...                [-Integer(24)*I + Integer(25), Integer(69)*I + Integer(89),        Integer(976), Integer(24)*I + Integer(6)],
...                [     -Integer(21)*I, -Integer(7)*I + Integer(15),  -Integer(24)*I + Integer(6),       Integer(28)]])
>>> A.is_hermitian()
True
>>> L = A.cholesky(); L
[                4.79...?                         0                       0        0]
[   0.62...? - 3.54...?*I                  5.00...?                       0        0]
[   5.21...? - 5.00...?*I   13.58...? + 10.72...?*I               24.98...?        0]
[             -4.37...?*I   -0.10...? -  0.85...?*I  -0.21...? + 0.37...?*I 2.81...?]
>>> L.parent()
Full MatrixSpace of 4 by 4 dense matrices over Algebraic Field
>>> (L*L.conjugate_transpose() - A.change_ring(QQbar)).norm() < Integer(10)**-Integer(10)
True

The field of algebraic numbers is an ideal setting for this computation:

sage: # needs sage.rings.number_field
sage: A = matrix(QQbar, [[        2,   4 + 2*I,   6 - 4*I],
....:                    [ -2*I + 4,        11, 10 - 12*I],
....:                    [  4*I + 6, 10 + 12*I,        37]])
sage: A.is_hermitian()
True
sage: L = A.cholesky()
sage: L
[                       1.414213562373095?          0                    0]
[2.828427124746190? - 1.414213562373095?*I          1                    0]
[4.242640687119285? + 2.828427124746190?*I   -2*I + 2   1.732050807568878?]
sage: L.parent()
Full MatrixSpace of 3 by 3 dense matrices over Algebraic Field
sage: L*L.conjugate_transpose() == A
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQbar, [[        Integer(2),   Integer(4) + Integer(2)*I,   Integer(6) - Integer(4)*I],
...                    [ -Integer(2)*I + Integer(4),        Integer(11), Integer(10) - Integer(12)*I],
...                    [  Integer(4)*I + Integer(6), Integer(10) + Integer(12)*I,        Integer(37)]])
>>> A.is_hermitian()
True
>>> L = A.cholesky()
>>> L
[                       1.414213562373095?          0                    0]
[2.828427124746190? - 1.414213562373095?*I          1                    0]
[4.242640687119285? + 2.828427124746190?*I   -2*I + 2   1.732050807568878?]
>>> L.parent()
Full MatrixSpace of 3 by 3 dense matrices over Algebraic Field
>>> L*L.conjugate_transpose() == A
True

Results are cached, hence immutable. Use the copy function if you need to make a change:

sage: A = matrix(QQ, [[ 4, -2,  4,  2],
....:                 [-2, 10, -2, -7],
....:                 [ 4, -2,  8,  4],
....:                 [ 2, -7,  4,  7]])
sage: L = A.cholesky()
sage: L.is_immutable()
True
sage: from copy import copy
sage: LC = copy(L)
sage: LC[0,0] = 1000
sage: LC
[1000    0    0    0]
[  -1    3    0    0]
[   2    0    2    0]
[   1   -2    1    1]
>>> from sage.all import *
>>> A = matrix(QQ, [[ Integer(4), -Integer(2),  Integer(4),  Integer(2)],
...                 [-Integer(2), Integer(10), -Integer(2), -Integer(7)],
...                 [ Integer(4), -Integer(2),  Integer(8),  Integer(4)],
...                 [ Integer(2), -Integer(7),  Integer(4),  Integer(7)]])
>>> L = A.cholesky()
>>> L.is_immutable()
True
>>> from copy import copy
>>> LC = copy(L)
>>> LC[Integer(0),Integer(0)] = Integer(1000)
>>> LC
[1000    0    0    0]
[  -1    3    0    0]
[   2    0    2    0]
[   1   -2    1    1]

The base ring need not be exact, although you should expect the result to be inexact (correct only in the norm) as well in that case:

sage: F = RealField(100)
sage: A = A = matrix(F, [[1.0, 2.0], [2.0, 6.0]])
sage: L = A.cholesky(); L
[ 1.000... 0.000...]
[ 2.000... 1.414...]
sage: (L*L.transpose() - A).norm() < 1e-10                                  # needs scipy
True
>>> from sage.all import *
>>> F = RealField(Integer(100))
>>> A = A = matrix(F, [[RealNumber('1.0'), RealNumber('2.0')], [RealNumber('2.0'), RealNumber('6.0')]])
>>> L = A.cholesky(); L
[ 1.000... 0.000...]
[ 2.000... 1.414...]
>>> (L*L.transpose() - A).norm() < RealNumber('1e-10')                                  # needs scipy
True

Even symbolic matrices can sometimes be factored:

sage: A = matrix(SR, [[pi,0], [0,pi]])                                      # needs sage.symbolic
sage: A.cholesky()                                                          # needs sage.symbolic
[sqrt(pi)        0]
[       0 sqrt(pi)]
>>> from sage.all import *
>>> A = matrix(SR, [[pi,Integer(0)], [Integer(0),pi]])                                      # needs sage.symbolic
>>> A.cholesky()                                                          # needs sage.symbolic
[sqrt(pi)        0]
[       0 sqrt(pi)]

There are a variety of situations which will prevent the computation of a Cholesky decomposition.

The base ring may not be able to be viewed as a subset of the complex numbers, implying that “Hermitian” is meaningless:

sage: A = matrix(Integers(6), [[2, 0], [0, 4]])
sage: A.cholesky()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.finite_rings.integer_mod.IntegerMod_int'
object has no attribute 'conjugate'
>>> from sage.all import *
>>> A = matrix(Integers(Integer(6)), [[Integer(2), Integer(0)], [Integer(0), Integer(4)]])
>>> A.cholesky()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.finite_rings.integer_mod.IntegerMod_int'
object has no attribute 'conjugate'

The matrix may not be Hermitian:

sage: F.<a> = FiniteField(5^4)                                              # needs sage.rings.finite_rings
sage: A = matrix(F, [[2+a^3, 3], [3, 3]])                                   # needs sage.rings.finite_rings
sage: A.cholesky()                                                          # needs sage.rings.finite_rings
Traceback (most recent call last):
...
ValueError: matrix is not Hermitian
>>> from sage.all import *
>>> F = FiniteField(Integer(5)**Integer(4), names=('a',)); (a,) = F._first_ngens(1)# needs sage.rings.finite_rings
>>> A = matrix(F, [[Integer(2)+a**Integer(3), Integer(3)], [Integer(3), Integer(3)]])                                   # needs sage.rings.finite_rings
>>> A.cholesky()                                                          # needs sage.rings.finite_rings
Traceback (most recent call last):
...
ValueError: matrix is not Hermitian

The matrix may not be positive-definite:

sage: # needs sage.rings.number_field
sage: C.<I> = QuadraticField(-1)
sage: B = matrix(C, [[      2, 4 - 2*I, 2 + 2*I],
....:                [4 + 2*I,       8,    10*I],
....:                [2 - 2*I,   -10*I,      -3]])
sage: B.is_positive_definite()
False
sage: B.cholesky()
Traceback (most recent call last):
...
ValueError: matrix is not positive definite
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> C = QuadraticField(-Integer(1), names=('I',)); (I,) = C._first_ngens(1)
>>> B = matrix(C, [[      Integer(2), Integer(4) - Integer(2)*I, Integer(2) + Integer(2)*I],
...                [Integer(4) + Integer(2)*I,       Integer(8),    Integer(10)*I],
...                [Integer(2) - Integer(2)*I,   -Integer(10)*I,      -Integer(3)]])
>>> B.is_positive_definite()
False
>>> B.cholesky()
Traceback (most recent call last):
...
ValueError: matrix is not positive definite

sage: A = matrix(QQ, [[21, 15, 12, -3],
....:                 [15, 12,  9,  12],
....:                 [12,  9,  7,  3],
....:                 [-3,  12,  3,  8]])
sage: A.is_positive_definite()
False
sage: A.cholesky()
Traceback (most recent call last):
...
ValueError: matrix is not positive definite
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(21), Integer(15), Integer(12), -Integer(3)],
...                 [Integer(15), Integer(12),  Integer(9),  Integer(12)],
...                 [Integer(12),  Integer(9),  Integer(7),  Integer(3)],
...                 [-Integer(3),  Integer(12),  Integer(3),  Integer(8)]])
>>> A.is_positive_definite()
False
>>> A.cholesky()
Traceback (most recent call last):
...
ValueError: matrix is not positive definite
column_module()[source]

Return the free module over the base ring spanned by the columns of this matrix.

EXAMPLES:

sage: t = matrix(QQ, 3, 3, range(9)); t
[0 1 2]
[3 4 5]
[6 7 8]
sage: t.column_module()
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -1]
[ 0  1  2]
>>> from sage.all import *
>>> t = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); t
[0 1 2]
[3 4 5]
[6 7 8]
>>> t.column_module()
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -1]
[ 0  1  2]
column_space()[source]

Return the vector space over the base ring spanned by the columns of this matrix.

EXAMPLES:

sage: M = MatrixSpace(QQ, 3, 3)
sage: A = M([1,9,-7, 4/5,4,3, 6,4,3])
sage: A.column_space()
Vector space of degree 3 and dimension 3 over Rational Field
Basis matrix:
[1 0 0]
[0 1 0]
[0 0 1]

sage: # needs sage.rings.real_mpfr sage.symbolic
sage: W = MatrixSpace(CC, 2, 2)
sage: B = W([1, 2 + 3*I, 4 + 5*I, 9]); B
[                     1.00000000000000 2.00000000000000 + 3.00000000000000*I]
[4.00000000000000 + 5.00000000000000*I                      9.00000000000000]
sage: B.column_space()
Vector space of degree 2 and dimension 2
 over Complex Field with 53 bits of precision
 Basis matrix:
 [ 1.00000000000000 0.000000000000000]
 [0.000000000000000  1.00000000000000]
>>> from sage.all import *
>>> M = MatrixSpace(QQ, Integer(3), Integer(3))
>>> A = M([Integer(1),Integer(9),-Integer(7), Integer(4)/Integer(5),Integer(4),Integer(3), Integer(6),Integer(4),Integer(3)])
>>> A.column_space()
Vector space of degree 3 and dimension 3 over Rational Field
Basis matrix:
[1 0 0]
[0 1 0]
[0 0 1]

>>> # needs sage.rings.real_mpfr sage.symbolic
>>> W = MatrixSpace(CC, Integer(2), Integer(2))
>>> B = W([Integer(1), Integer(2) + Integer(3)*I, Integer(4) + Integer(5)*I, Integer(9)]); B
[                     1.00000000000000 2.00000000000000 + 3.00000000000000*I]
[4.00000000000000 + 5.00000000000000*I                      9.00000000000000]
>>> B.column_space()
Vector space of degree 2 and dimension 2
 over Complex Field with 53 bits of precision
 Basis matrix:
 [ 1.00000000000000 0.000000000000000]
 [0.000000000000000  1.00000000000000]
conjugate()[source]

Return the conjugate of self, i.e. the matrix whose entries are the conjugates of the entries of self.

EXAMPLES:

sage: # needs sage.rings.complex_double sage.symbolic
sage: A = matrix(CDF, [[1+I,1],[0,2*I]])
sage: A.conjugate()
[1.0 - 1.0*I         1.0]
[        0.0      -2.0*I]
>>> from sage.all import *
>>> # needs sage.rings.complex_double sage.symbolic
>>> A = matrix(CDF, [[Integer(1)+I,Integer(1)],[Integer(0),Integer(2)*I]])
>>> A.conjugate()
[1.0 - 1.0*I         1.0]
[        0.0      -2.0*I]

A matrix over a not-totally-real number field:

sage: x = polygen(ZZ, 'x')
sage: K.<j> = NumberField(x^2 + 5)                                          # needs sage.rings.number_field
sage: M = matrix(K, [[1+j,1], [0,2*j]])                                     # needs sage.rings.number_field
sage: M.conjugate()                                                         # needs sage.rings.number_field
[-j + 1      1]
[     0   -2*j]
>>> from sage.all import *
>>> x = polygen(ZZ, 'x')
>>> K = NumberField(x**Integer(2) + Integer(5), names=('j',)); (j,) = K._first_ngens(1)# needs sage.rings.number_field
>>> M = matrix(K, [[Integer(1)+j,Integer(1)], [Integer(0),Integer(2)*j]])                                     # needs sage.rings.number_field
>>> M.conjugate()                                                         # needs sage.rings.number_field
[-j + 1      1]
[     0   -2*j]

There is a shortcut for the conjugate:

sage: M.C                                                                   # needs sage.rings.number_field
[-j + 1      1]
[     0   -2*j]
>>> from sage.all import *
>>> M.C                                                                   # needs sage.rings.number_field
[-j + 1      1]
[     0   -2*j]

There is also a shortcut for the conjugate transpose, or “Hermitian transpose”:

sage: M.H                                                                   # needs sage.rings.number_field
[-j + 1      0]
[     1   -2*j]
>>> from sage.all import *
>>> M.H                                                                   # needs sage.rings.number_field
[-j + 1      0]
[     1   -2*j]

Conjugates work (trivially) for matrices over rings that embed canonically into the real numbers:

sage: M = random_matrix(ZZ, 2)
sage: M == M.conjugate()
True
sage: M = random_matrix(QQ, 3)
sage: M == M.conjugate()
True
sage: M = random_matrix(RR, 2)
sage: M == M.conjugate()
True
>>> from sage.all import *
>>> M = random_matrix(ZZ, Integer(2))
>>> M == M.conjugate()
True
>>> M = random_matrix(QQ, Integer(3))
>>> M == M.conjugate()
True
>>> M = random_matrix(RR, Integer(2))
>>> M == M.conjugate()
True
conjugate_transpose()[source]

Return the transpose of self after each entry has been converted to its complex conjugate.

Note

This function is sometimes known as the “adjoint” of a matrix, though there is substantial variation and some confusion with the use of that term.

OUTPUT:

A matrix formed by taking the complex conjugate of every entry of self and then transposing the resulting matrix.

Complex conjugation is implemented for many subfields of the complex numbers. See the examples below, or more at conjugate().

EXAMPLES:

sage: M = matrix(SR, 2, 2, [[2-I, 3+4*I], [9-6*I, 5*I]])                    # needs sage.symbolic
sage: M.base_ring()                                                         # needs sage.symbolic
Symbolic Ring
sage: M.conjugate_transpose()                                               # needs sage.symbolic
[   I + 2  6*I + 9]
[-4*I + 3     -5*I]

sage: # needs sage.rings.real_mpfr sage.symbolic
sage: P = matrix(CC, 3, 2, [0.95-0.63*I, 0.84+0.13*I,
....:                       0.94+0.23*I, 0.23+0.59*I,
....:                       0.52-0.41*I, -0.50+0.90*I])
sage: P.base_ring()
Complex Field with 53 bits of precision
sage: P.conjugate_transpose()
[ 0.950... + 0.630...*I  0.940... - 0.230...*I  0.520... + 0.410...*I]
[ 0.840... - 0.130...*I  0.230... - 0.590...*I -0.500... - 0.900...*I]
>>> from sage.all import *
>>> M = matrix(SR, Integer(2), Integer(2), [[Integer(2)-I, Integer(3)+Integer(4)*I], [Integer(9)-Integer(6)*I, Integer(5)*I]])                    # needs sage.symbolic
>>> M.base_ring()                                                         # needs sage.symbolic
Symbolic Ring
>>> M.conjugate_transpose()                                               # needs sage.symbolic
[   I + 2  6*I + 9]
[-4*I + 3     -5*I]

>>> # needs sage.rings.real_mpfr sage.symbolic
>>> P = matrix(CC, Integer(3), Integer(2), [RealNumber('0.95')-RealNumber('0.63')*I, RealNumber('0.84')+RealNumber('0.13')*I,
...                       RealNumber('0.94')+RealNumber('0.23')*I, RealNumber('0.23')+RealNumber('0.59')*I,
...                       RealNumber('0.52')-RealNumber('0.41')*I, -RealNumber('0.50')+RealNumber('0.90')*I])
>>> P.base_ring()
Complex Field with 53 bits of precision
>>> P.conjugate_transpose()
[ 0.950... + 0.630...*I  0.940... - 0.230...*I  0.520... + 0.410...*I]
[ 0.840... - 0.130...*I  0.230... - 0.590...*I -0.500... - 0.900...*I]

There is also a shortcut for the conjugate transpose, or “Hermitian transpose”:

sage: M.H                                                                   # needs sage.symbolic
[   I + 2  6*I + 9]
[-4*I + 3     -5*I]
>>> from sage.all import *
>>> M.H                                                                   # needs sage.symbolic
[   I + 2  6*I + 9]
[-4*I + 3     -5*I]

Matrices over base rings that can be embedded in the real numbers will behave as expected.

sage: P = random_matrix(QQ, 3, 4)
sage: P.conjugate_transpose() == P.transpose()
True
>>> from sage.all import *
>>> P = random_matrix(QQ, Integer(3), Integer(4))
>>> P.conjugate_transpose() == P.transpose()
True

The conjugate of a matrix is formed by taking conjugates of all the entries. Some specialized subfields of the complex numbers are implemented in Sage and complex conjugation can be applied. (Matrices over quadratic number fields are another class of examples.)

sage: # needs sage.rings.number_field
sage: C = CyclotomicField(5)
sage: a = C.gen(); a
zeta5
sage: CC(a)
0.309016994374947 + 0.951056516295154*I
sage: M = matrix(C, 1, 2, [a^2, a+a^3])
sage: M.conjugate_transpose()
[             zeta5^3]
[-zeta5^3 - zeta5 - 1]
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> C = CyclotomicField(Integer(5))
>>> a = C.gen(); a
zeta5
>>> CC(a)
0.309016994374947 + 0.951056516295154*I
>>> M = matrix(C, Integer(1), Integer(2), [a**Integer(2), a+a**Integer(3)])
>>> M.conjugate_transpose()
[             zeta5^3]
[-zeta5^3 - zeta5 - 1]

Furthermore, this method can be applied to matrices over quadratic extensions of finite fields:

sage: F.<a> = GF(9,'a')                                                     # needs sage.rings.finite_rings
sage: N = matrix(F, 2, [0,a,-a,1]); N                                       # needs sage.rings.finite_rings
[  0   a]
[2*a   1]
sage: N.conjugate_transpose()                                               # needs sage.rings.finite_rings
[      0   a + 2]
[2*a + 1       1]
>>> from sage.all import *
>>> F = GF(Integer(9),'a', names=('a',)); (a,) = F._first_ngens(1)# needs sage.rings.finite_rings
>>> N = matrix(F, Integer(2), [Integer(0),a,-a,Integer(1)]); N                                       # needs sage.rings.finite_rings
[  0   a]
[2*a   1]
>>> N.conjugate_transpose()                                               # needs sage.rings.finite_rings
[      0   a + 2]
[2*a + 1       1]

Conjugation does not make sense over rings not containing complex numbers or finite fields which are not a quadratic extension:

sage: N = matrix(GF(5), 2, [0,1,2,3])
sage: N.conjugate_transpose()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.finite_rings.integer_mod.IntegerMod_int' object
has no attribute 'conjugate'...
>>> from sage.all import *
>>> N = matrix(GF(Integer(5)), Integer(2), [Integer(0),Integer(1),Integer(2),Integer(3)])
>>> N.conjugate_transpose()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.finite_rings.integer_mod.IntegerMod_int' object
has no attribute 'conjugate'...
cyclic_subspace(v, var=None, basis='echelon')[source]

Create a cyclic subspace for a vector, and optionally, a minimal polynomial for the iterated powers.

These subspaces are also known as Krylov subspaces. They are spanned by the vectors

\[\{v, Av, A^2v, A^3v, \dots \}\]

INPUT:

  • self – a square matrix with entries from a field

  • v – a vector with a degree equal to the size of the matrix and entries compatible with the entries of the matrix

  • var – (default: None) if specified as a string or a generator of a polynomial ring, then this will be used to construct a polynomial reflecting a relation of linear dependence on the powers \(A^iv\) and this will cause the polynomial to be returned along with the subspace. A generator must create polynomials with coefficients from the same field as the matrix entries.

  • basis – (default: echelon) the basis for the subspace is “echelonized” by default, but the keyword ‘iterates’ will return a subspace with a user basis equal to the largest linearly independent set \(\{v, Av, A^2v, A^3v, \dots, A^{k-1}v \}\).

OUTPUT:

Suppose \(k\) is the smallest power such that \(\{v, Av, A^2v, A^3v, \dots, A^{k}v \}\) is linearly dependent. Then the subspace returned will have dimension \(k\) and be spanned by the powers \(0\) through \(k-1\).

If a polynomial is requested through the use of the var keyword, then a pair is returned, with the polynomial first and the subspace second. The polynomial is the unique monic polynomial whose coefficients provide a relation of linear dependence on the first \(k\) powers.

For less convenient, but more flexible output, see the helper method “_cyclic_subspace” in this module.

EXAMPLES:

sage: A = matrix(QQ, [[5,4,2,1],[0,1,-1,-1],[-1,-1,3,0],[1,1,-1,2]])
sage: v = vector(QQ, [0,1,0,0])
sage: E = A.cyclic_subspace(v); E
Vector space of degree 4 and dimension 3 over Rational Field
Basis matrix:
[ 1  0  0  0]
[ 0  1  0  0]
[ 0  0  1 -1]
sage: F = A.cyclic_subspace(v, basis='iterates'); F
Vector space of degree 4 and dimension 3 over Rational Field
User basis matrix:
[ 0  1  0  0]
[ 4  1 -1  1]
[23  1 -8  8]
sage: E == F
True
sage: p, S = A.cyclic_subspace(v, var='T'); p
T^3 - 9*T^2 + 24*T - 16
sage: p.degree() == E.dimension()
True
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(5),Integer(4),Integer(2),Integer(1)],[Integer(0),Integer(1),-Integer(1),-Integer(1)],[-Integer(1),-Integer(1),Integer(3),Integer(0)],[Integer(1),Integer(1),-Integer(1),Integer(2)]])
>>> v = vector(QQ, [Integer(0),Integer(1),Integer(0),Integer(0)])
>>> E = A.cyclic_subspace(v); E
Vector space of degree 4 and dimension 3 over Rational Field
Basis matrix:
[ 1  0  0  0]
[ 0  1  0  0]
[ 0  0  1 -1]
>>> F = A.cyclic_subspace(v, basis='iterates'); F
Vector space of degree 4 and dimension 3 over Rational Field
User basis matrix:
[ 0  1  0  0]
[ 4  1 -1  1]
[23  1 -8  8]
>>> E == F
True
>>> p, S = A.cyclic_subspace(v, var='T'); p
T^3 - 9*T^2 + 24*T - 16
>>> p.degree() == E.dimension()
True

The polynomial has coefficients that yield a non-trivial relation of linear dependence on the iterates. Or, equivalently, evaluating the polynomial with the matrix will create a matrix that annihilates the vector.

sage: A = matrix(QQ, [[15, 37/3, -16, -104/3, -29, -7/3, 35, 2/3, -29/3, -1/3],
....:                 [ 2, 9, -1, -6, -6, 0, 7, 0, -2, 0],
....:                 [24, 74/3, -29, -208/3, -58, -14/3, 70, 4/3, -58/3, -2/3],
....:                 [-6, -19, 3, 21, 19, 0, -21, 0, 6, 0],
....:                 [2, 6, -1, -6, -3, 0, 7, 0, -2, 0],
....:                 [-96, -296/3, 128, 832/3, 232, 65/3, -279, -16/3, 232/3, 8/3],
....:                 [0, 0, 0, 0, 0, 0, 3, 0, 0, 0],
....:                 [20, 26/3, -30, -199/3, -42, -14/3, 70, 13/3, -55/3, -2/3],
....:                 [18, 57, -9, -54, -57, 0, 63, 0, -15, 0],
....:                 [0, 0, 0, 0, 0, 0, 0, 0, 0, 3]])
sage: u = zero_vector(QQ, 10); u[0] = 1
sage: p, S = A.cyclic_subspace(u, var='t', basis='iterates')
sage: S
Vector space of degree 10 and dimension 3 over Rational Field
User basis matrix:
[   1    0    0    0    0    0    0    0    0    0]
[  15    2   24   -6    2  -96    0   20   18    0]
[  79   12  140  -36   12 -560    0  116  108    0]
sage: p
t^3 - 9*t^2 + 27*t - 27
sage: k = p.degree()
sage: coeffs = p.list()
sage: iterates = S.basis() + [A^k*u]
sage: sum(coeffs[i]*iterates[i] for i in range(k+1))
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
sage: u in p(A).right_kernel()
True
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(15), Integer(37)/Integer(3), -Integer(16), -Integer(104)/Integer(3), -Integer(29), -Integer(7)/Integer(3), Integer(35), Integer(2)/Integer(3), -Integer(29)/Integer(3), -Integer(1)/Integer(3)],
...                 [ Integer(2), Integer(9), -Integer(1), -Integer(6), -Integer(6), Integer(0), Integer(7), Integer(0), -Integer(2), Integer(0)],
...                 [Integer(24), Integer(74)/Integer(3), -Integer(29), -Integer(208)/Integer(3), -Integer(58), -Integer(14)/Integer(3), Integer(70), Integer(4)/Integer(3), -Integer(58)/Integer(3), -Integer(2)/Integer(3)],
...                 [-Integer(6), -Integer(19), Integer(3), Integer(21), Integer(19), Integer(0), -Integer(21), Integer(0), Integer(6), Integer(0)],
...                 [Integer(2), Integer(6), -Integer(1), -Integer(6), -Integer(3), Integer(0), Integer(7), Integer(0), -Integer(2), Integer(0)],
...                 [-Integer(96), -Integer(296)/Integer(3), Integer(128), Integer(832)/Integer(3), Integer(232), Integer(65)/Integer(3), -Integer(279), -Integer(16)/Integer(3), Integer(232)/Integer(3), Integer(8)/Integer(3)],
...                 [Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(3), Integer(0), Integer(0), Integer(0)],
...                 [Integer(20), Integer(26)/Integer(3), -Integer(30), -Integer(199)/Integer(3), -Integer(42), -Integer(14)/Integer(3), Integer(70), Integer(13)/Integer(3), -Integer(55)/Integer(3), -Integer(2)/Integer(3)],
...                 [Integer(18), Integer(57), -Integer(9), -Integer(54), -Integer(57), Integer(0), Integer(63), Integer(0), -Integer(15), Integer(0)],
...                 [Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0), Integer(3)]])
>>> u = zero_vector(QQ, Integer(10)); u[Integer(0)] = Integer(1)
>>> p, S = A.cyclic_subspace(u, var='t', basis='iterates')
>>> S
Vector space of degree 10 and dimension 3 over Rational Field
User basis matrix:
[   1    0    0    0    0    0    0    0    0    0]
[  15    2   24   -6    2  -96    0   20   18    0]
[  79   12  140  -36   12 -560    0  116  108    0]
>>> p
t^3 - 9*t^2 + 27*t - 27
>>> k = p.degree()
>>> coeffs = p.list()
>>> iterates = S.basis() + [A**k*u]
>>> sum(coeffs[i]*iterates[i] for i in range(k+Integer(1)))
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
>>> u in p(A).right_kernel()
True
decomposition(algorithm='spin', is_diagonalizable=False, dual=False)[source]

Return the decomposition of the free module on which this matrix A acts from the right (i.e., the action is x goes to x A), along with whether this matrix acts irreducibly on each factor. The factors are guaranteed to be sorted in the same way as the corresponding factors of the characteristic polynomial.

Let A be the matrix acting from the on the vector space V of column vectors. Assume that A is square. This function computes maximal subspaces W_1, …, W_n corresponding to Galois conjugacy classes of eigenvalues of A. More precisely, let \(f(X)\) be the characteristic polynomial of A. This function computes the subspace \(W_i = ker(g_(A)^n)\), where \(g_i(X)\) is an irreducible factor of \(f(X)\) and \(g_i(X)\) exactly divides \(f(X)\). If the optional parameter is_diagonalizable is True, then we let \(W_i = ker(g(A))\), since then we know that \(ker(g(A)) = ker(g(A)^n)\).

INPUT:

  • self – a matrix

  • algorithm – string (default: 'spin'); 'spin': involves iterating the action of self on a vector. 'kernel': naively just compute \(ker(f_i(A))\) for each factor \(f_i\).

  • dual – boolean (default: False); if True, also returns the corresponding decomposition of V under the action of the transpose of A. The factors are guaranteed to correspond.

  • is_diagonalizable – if the matrix is known to be diagonalizable, set this to True, which might speed up the algorithm in some cases.

Note

If the base ring is not a field, the kernel algorithm is used.

OUTPUT:

  • Sequence – list of pairs (V,t), where V is a vector spaces and t is a boolean, and t is True exactly when the charpoly of self on V is irreducible.

  • (optional) list – list of pairs (W,t), where W is a vector space and t is a boolean, and t is True exactly when the charpoly of the transpose of self on W is irreducible.

EXAMPLES:

sage: A = matrix(ZZ, 4, [3,4,5,6, 7,3,8,10, 14,5,6,7, 2,2,10,9])
sage: B = matrix(QQ, 6, 6, range(36))
sage: B*11
[  0  11  22  33  44  55]
[ 66  77  88  99 110 121]
[132 143 154 165 176 187]
[198 209 220 231 242 253]
[264 275 286 297 308 319]
[330 341 352 363 374 385]
sage: A.decomposition()                                                     # needs sage.libs.pari
[ (Ambient free module of rank 4
    over the principal ideal domain Integer Ring,
   True) ]
sage: B.decomposition()                                                     # needs sage.libs.pari
[ (Vector space of degree 6 and dimension 2 over Rational Field
    Basis matrix:
    [ 1  0 -1 -2 -3 -4]
    [ 0  1  2  3  4  5],
   True),
  (Vector space of degree 6 and dimension 4 over Rational Field
    Basis matrix:
    [ 1  0  0  0 -5  4]
    [ 0  1  0  0 -4  3]
    [ 0  0  1  0 -3  2]
    [ 0  0  0  1 -2  1],
   False) ]
>>> from sage.all import *
>>> A = matrix(ZZ, Integer(4), [Integer(3),Integer(4),Integer(5),Integer(6), Integer(7),Integer(3),Integer(8),Integer(10), Integer(14),Integer(5),Integer(6),Integer(7), Integer(2),Integer(2),Integer(10),Integer(9)])
>>> B = matrix(QQ, Integer(6), Integer(6), range(Integer(36)))
>>> B*Integer(11)
[  0  11  22  33  44  55]
[ 66  77  88  99 110 121]
[132 143 154 165 176 187]
[198 209 220 231 242 253]
[264 275 286 297 308 319]
[330 341 352 363 374 385]
>>> A.decomposition()                                                     # needs sage.libs.pari
[ (Ambient free module of rank 4
    over the principal ideal domain Integer Ring,
   True) ]
>>> B.decomposition()                                                     # needs sage.libs.pari
[ (Vector space of degree 6 and dimension 2 over Rational Field
    Basis matrix:
    [ 1  0 -1 -2 -3 -4]
    [ 0  1  2  3  4  5],
   True),
  (Vector space of degree 6 and dimension 4 over Rational Field
    Basis matrix:
    [ 1  0  0  0 -5  4]
    [ 0  1  0  0 -4  3]
    [ 0  0  1  0 -3  2]
    [ 0  0  0  1 -2  1],
   False) ]
decomposition_of_subspace(M, check_restrict=True, **kwds)[source]

Suppose the right action of self on M leaves M invariant. Return the decomposition of M as a list of pairs (W, is_irred) where is_irred is True if the charpoly of self acting on the factor W is irreducible.

Additional inputs besides M are passed onto the decomposition command.

INPUT:

  • M – a subspace of the free module self acts on

  • check_restrict – boolean (default: True); call restrict with or without check

  • kwds – keywords that will be forwarded to decomposition()

EXAMPLES:

sage: # needs sage.libs.pari
sage: t = matrix(QQ, 3, [3, 0, -2, 0, -2, 0, 0, 0, 0]); t
[ 3  0 -2]
[ 0 -2  0]
[ 0  0  0]
sage: t.fcp('X')                # factored charpoly
(X - 3) * X * (X + 2)
sage: v = kernel(t*(t+2)); v    # an invariant subspace
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[0 1 0]
[0 0 1]
sage: D = t.decomposition_of_subspace(v); D
[ (Vector space of degree 3 and dimension 1 over Rational Field
    Basis matrix: [0 0 1],
   True),
  (Vector space of degree 3 and dimension 1 over Rational Field
    Basis matrix: [0 1 0],
   True) ]
sage: t.restrict(D[0][0])
[0]
sage: t.restrict(D[1][0])
[-2]
>>> from sage.all import *
>>> # needs sage.libs.pari
>>> t = matrix(QQ, Integer(3), [Integer(3), Integer(0), -Integer(2), Integer(0), -Integer(2), Integer(0), Integer(0), Integer(0), Integer(0)]); t
[ 3  0 -2]
[ 0 -2  0]
[ 0  0  0]
>>> t.fcp('X')                # factored charpoly
(X - 3) * X * (X + 2)
>>> v = kernel(t*(t+Integer(2))); v    # an invariant subspace
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[0 1 0]
[0 0 1]
>>> D = t.decomposition_of_subspace(v); D
[ (Vector space of degree 3 and dimension 1 over Rational Field
    Basis matrix: [0 0 1],
   True),
  (Vector space of degree 3 and dimension 1 over Rational Field
    Basis matrix: [0 1 0],
   True) ]
>>> t.restrict(D[Integer(0)][Integer(0)])
[0]
>>> t.restrict(D[Integer(1)][Integer(0)])
[-2]

We do a decomposition over ZZ:

sage: a = matrix(ZZ, 6, [0, 0, -2, 0, 2, 0,
....:                    2, -4, -2, 0, 2, 0,
....:                    0, 0, -2, -2, 0, 0,
....:                    2, 0, -2, -4, 2, -2,
....:                    0, 2, 0, -2, -2, 0,
....:                    0, 2, 0, -2, 0, 0])
sage: a.decomposition_of_subspace(ZZ^6)                                     # needs sage.libs.pari
[ (Free module of degree 6 and rank 2 over Integer Ring
    Echelon basis matrix:
    [ 1  0  1 -1  1 -1]
    [ 0  1  0 -1  2 -1],
   False),
  (Free module of degree 6 and rank 4 over Integer Ring
    Echelon basis matrix:
    [ 1  0 -1  0  1  0]
    [ 0  1  0  0  0  0]
    [ 0  0  0  1  0  0]
    [ 0  0  0  0  0  1],
   False) ]
>>> from sage.all import *
>>> a = matrix(ZZ, Integer(6), [Integer(0), Integer(0), -Integer(2), Integer(0), Integer(2), Integer(0),
...                    Integer(2), -Integer(4), -Integer(2), Integer(0), Integer(2), Integer(0),
...                    Integer(0), Integer(0), -Integer(2), -Integer(2), Integer(0), Integer(0),
...                    Integer(2), Integer(0), -Integer(2), -Integer(4), Integer(2), -Integer(2),
...                    Integer(0), Integer(2), Integer(0), -Integer(2), -Integer(2), Integer(0),
...                    Integer(0), Integer(2), Integer(0), -Integer(2), Integer(0), Integer(0)])
>>> a.decomposition_of_subspace(ZZ**Integer(6))                                     # needs sage.libs.pari
[ (Free module of degree 6 and rank 2 over Integer Ring
    Echelon basis matrix:
    [ 1  0  1 -1  1 -1]
    [ 0  1  0 -1  2 -1],
   False),
  (Free module of degree 6 and rank 4 over Integer Ring
    Echelon basis matrix:
    [ 1  0 -1  0  1  0]
    [ 0  1  0  0  0  0]
    [ 0  0  0  1  0  0]
    [ 0  0  0  0  0  1],
   False) ]
denominator()[source]

Return the least common multiple of the denominators of the elements of self.

If there is no denominator function for the base field, or no LCM function for the denominators, raise a TypeError.

EXAMPLES:

sage: A = MatrixSpace(QQ, 2)([1/2, 1/3, 1/5, 1/7])
sage: A.denominator()
210
>>> from sage.all import *
>>> A = MatrixSpace(QQ, Integer(2))([Integer(1)/Integer(2), Integer(1)/Integer(3), Integer(1)/Integer(5), Integer(1)/Integer(7)])
>>> A.denominator()
210

A trivial example:

sage: A = matrix(QQ, 0,2)
sage: A.denominator()
1
>>> from sage.all import *
>>> A = matrix(QQ, Integer(0),Integer(2))
>>> A.denominator()
1

Denominators are not defined for real numbers:

sage: A = MatrixSpace(RealField(),2)([1,2,3,4])
sage: A.denominator()
Traceback (most recent call last):
...
TypeError: denominator not defined for elements of the base ring
>>> from sage.all import *
>>> A = MatrixSpace(RealField(),Integer(2))([Integer(1),Integer(2),Integer(3),Integer(4)])
>>> A.denominator()
Traceback (most recent call last):
...
TypeError: denominator not defined for elements of the base ring

We can even compute the denominator of matrix over the fraction field of \(\ZZ[x]\).

sage: K.<x> = Frac(ZZ['x'])
sage: A = MatrixSpace(K,2)([1/x, 2/(x+1), 1, 5/(x^3)])
sage: A.denominator()
x^4 + x^3
>>> from sage.all import *
>>> K = Frac(ZZ['x'], names=('x',)); (x,) = K._first_ngens(1)
>>> A = MatrixSpace(K,Integer(2))([Integer(1)/x, Integer(2)/(x+Integer(1)), Integer(1), Integer(5)/(x**Integer(3))])
>>> A.denominator()
x^4 + x^3

Here’s an example involving a cyclotomic field:

sage: # needs sage.rings.number_field
sage: K.<z> = CyclotomicField(3)
sage: M = MatrixSpace(K, 3, sparse=True)
sage: A = M([(1+z)/3, (2+z)/3, z/3, 1, 1+z, -2, 1, 5, -1+z])
sage: print(A)
[1/3*z + 1/3 1/3*z + 2/3       1/3*z]
[          1       z + 1          -2]
[          1           5       z - 1]
sage: print(A.denominator())
3
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> K = CyclotomicField(Integer(3), names=('z',)); (z,) = K._first_ngens(1)
>>> M = MatrixSpace(K, Integer(3), sparse=True)
>>> A = M([(Integer(1)+z)/Integer(3), (Integer(2)+z)/Integer(3), z/Integer(3), Integer(1), Integer(1)+z, -Integer(2), Integer(1), Integer(5), -Integer(1)+z])
>>> print(A)
[1/3*z + 1/3 1/3*z + 2/3       1/3*z]
[          1       z + 1          -2]
[          1           5       z - 1]
>>> print(A.denominator())
3
density()[source]

Return the density of the matrix.

By density we understand the ratio of the number of nonzero positions and the self.nrows() * self.ncols(), i.e. the number of possible nonzero positions.

EXAMPLES:

First, note that the density parameter does not ensure the density of a matrix, it is only an upper bound.

sage: A = random_matrix(GF(127), 200, 200, density=0.3)
sage: A.density() <= 0.3
True
>>> from sage.all import *
>>> A = random_matrix(GF(Integer(127)), Integer(200), Integer(200), density=RealNumber('0.3'))
>>> A.density() <= RealNumber('0.3')
True

sage: A = matrix(QQ, 3,3, [0,1,2,3,0,0,6,7,8])
sage: A.density()
2/3
>>> from sage.all import *
>>> A = matrix(QQ, Integer(3),Integer(3), [Integer(0),Integer(1),Integer(2),Integer(3),Integer(0),Integer(0),Integer(6),Integer(7),Integer(8)])
>>> A.density()
2/3

sage: a = matrix([[],[],[],[]])
sage: a.density()
0
>>> from sage.all import *
>>> a = matrix([[],[],[],[]])
>>> a.density()
0
derivative(*args)[source]

Derivative with respect to variables supplied in args.

Multiple variables and iteration counts may be supplied; see documentation for the global derivative() function for more details.

EXAMPLES:

sage: # needs sage.symbolic
sage: v = vector([1,x,x^2])
sage: v.derivative(x)
(0, 1, 2*x)
sage: type(v.derivative(x)) == type(v)
True
sage: v = vector([1,x,x^2], sparse=True)
sage: v.derivative(x)
(0, 1, 2*x)
sage: type(v.derivative(x)) == type(v)
True
sage: v.derivative(x,x)
(0, 0, 2)
>>> from sage.all import *
>>> # needs sage.symbolic
>>> v = vector([Integer(1),x,x**Integer(2)])
>>> v.derivative(x)
(0, 1, 2*x)
>>> type(v.derivative(x)) == type(v)
True
>>> v = vector([Integer(1),x,x**Integer(2)], sparse=True)
>>> v.derivative(x)
(0, 1, 2*x)
>>> type(v.derivative(x)) == type(v)
True
>>> v.derivative(x,x)
(0, 0, 2)
det(*args, **kwds)[source]

Synonym for self.determinant(…).

EXAMPLES:

sage: A = MatrixSpace(Integers(8), 3)([1,7,3, 1,1,1, 3,4,5])
sage: A.det()
6
>>> from sage.all import *
>>> A = MatrixSpace(Integers(Integer(8)), Integer(3))([Integer(1),Integer(7),Integer(3), Integer(1),Integer(1),Integer(1), Integer(3),Integer(4),Integer(5)])
>>> A.det()
6
determinant(algorithm=None)[source]

Return the determinant of self.

ALGORITHM:

If the base ring has a method _matrix_determinant(), we call it.

Otherwise, for small matrices (n less than 4), this is computed using the naive formula. In the specific case of matrices over the integers modulo a non-prime, the determinant of a lift is computed over the integers. In general, the characteristic polynomial is computed either using the Hessenberg form (specified by 'hessenberg') or the generic division-free algorithm (specified by 'df'). When the base ring is an exact field, the default choice is 'hessenberg', otherwise it is 'df'. Note that for matrices over most rings, more sophisticated algorithms can be used. (Type A.determinant? to see what is done for a specific matrix A.)

INPUT:

  • algorithm – string; one of
    • 'df' – generic O(n^4) division-free algorithm

    • 'hessenberg' – use the Hessenberg form of the matrix

EXAMPLES:

sage: A = MatrixSpace(Integers(8), 3)([1,7,3, 1,1,1, 3,4,5])
sage: A.determinant()
6
sage: A.determinant() is A.determinant()
True
sage: A[0,0] = 10
sage: A.determinant()
7
>>> from sage.all import *
>>> A = MatrixSpace(Integers(Integer(8)), Integer(3))([Integer(1),Integer(7),Integer(3), Integer(1),Integer(1),Integer(1), Integer(3),Integer(4),Integer(5)])
>>> A.determinant()
6
>>> A.determinant() is A.determinant()
True
>>> A[Integer(0),Integer(0)] = Integer(10)
>>> A.determinant()
7

We compute the determinant of the arbitrary 3x3 matrix:

sage: R = PolynomialRing(QQ, 9, 'x')
sage: A = matrix(R, 3, R.gens())
sage: A
[x0 x1 x2]
[x3 x4 x5]
[x6 x7 x8]
sage: A.determinant()
-x2*x4*x6 + x1*x5*x6 + x2*x3*x7 - x0*x5*x7 - x1*x3*x8 + x0*x4*x8
>>> from sage.all import *
>>> R = PolynomialRing(QQ, Integer(9), 'x')
>>> A = matrix(R, Integer(3), R.gens())
>>> A
[x0 x1 x2]
[x3 x4 x5]
[x6 x7 x8]
>>> A.determinant()
-x2*x4*x6 + x1*x5*x6 + x2*x3*x7 - x0*x5*x7 - x1*x3*x8 + x0*x4*x8

We create a matrix over \(\ZZ[x,y]\) and compute its determinant.

sage: R.<x,y> = PolynomialRing(IntegerRing(), 2)
sage: A = MatrixSpace(R,2)([x, y, x**2, y**2])
sage: A.determinant()
-x^2*y + x*y^2
>>> from sage.all import *
>>> R = PolynomialRing(IntegerRing(), Integer(2), names=('x', 'y',)); (x, y,) = R._first_ngens(2)
>>> A = MatrixSpace(R,Integer(2))([x, y, x**Integer(2), y**Integer(2)])
>>> A.determinant()
-x^2*y + x*y^2

A matrix over a non-domain:

sage: m = matrix(Integers(4), 2, [1,2,2,3])
sage: m.determinant()
3
>>> from sage.all import *
>>> m = matrix(Integers(Integer(4)), Integer(2), [Integer(1),Integer(2),Integer(2),Integer(3)])
>>> m.determinant()
3
diagonal()[source]

Return the diagonal entries of self.

OUTPUT:

A list containing the entries of the matrix that have equal row and column indices, in order of the indices. Behavior is not limited to square matrices.

EXAMPLES:

sage: A = matrix([[2,5], [3,7]]); A
[2 5]
[3 7]
sage: A.diagonal()
[2, 7]
>>> from sage.all import *
>>> A = matrix([[Integer(2),Integer(5)], [Integer(3),Integer(7)]]); A
[2 5]
[3 7]
>>> A.diagonal()
[2, 7]

Two rectangular matrices.

sage: B = matrix(3, 7, range(21)); B
[ 0  1  2  3  4  5  6]
[ 7  8  9 10 11 12 13]
[14 15 16 17 18 19 20]
sage: B.diagonal()
[0, 8, 16]

sage: C = matrix(3, 2, range(6)); C
[0 1]
[2 3]
[4 5]
sage: C.diagonal()
[0, 3]
>>> from sage.all import *
>>> B = matrix(Integer(3), Integer(7), range(Integer(21))); B
[ 0  1  2  3  4  5  6]
[ 7  8  9 10 11 12 13]
[14 15 16 17 18 19 20]
>>> B.diagonal()
[0, 8, 16]

>>> C = matrix(Integer(3), Integer(2), range(Integer(6))); C
[0 1]
[2 3]
[4 5]
>>> C.diagonal()
[0, 3]

Empty matrices behave properly.

sage: E = matrix(0, 5, []); E
[]
sage: E.diagonal()
[]
>>> from sage.all import *
>>> E = matrix(Integer(0), Integer(5), []); E
[]
>>> E.diagonal()
[]
diagonalization(base_field=None)[source]

Return a diagonal matrix similar to self along with the transformation matrix.

INPUT:

  • base_field – if given, self is regarded as a matrix over it

OUTPUT: a diagonal matrix \(D\) and an invertible matrix \(P\) such that \(P^{-1}AP=D\), if self is a diagonalizable matrix \(A\).

EXAMPLES:

sage: # needs sage.libs.pari
sage: A = matrix(QQ, 4, [-4, 6, 3, 3, -3, 5, 3, 3, 3, -6, -4, -3, -3, 6, 3, 2])
sage: A
[-4  6  3  3]
[-3  5  3  3]
[ 3 -6 -4 -3]
[-3  6  3  2]
sage: A.is_diagonalizable()
True
sage: A.diagonalization()
(
[ 2  0  0  0]  [ 1  1  0  0]
[ 0 -1  0  0]  [ 1  0  1  0]
[ 0  0 -1  0]  [-1  0  0  1]
[ 0  0  0 -1], [ 1  1 -2 -1]
)
sage: D, P = A.diagonalization()
sage: P^-1*A*P == D
True

sage: # needs sage.libs.pari
sage: A = matrix(QQ, 2, [0, 2, 1, 0])
sage: A.is_diagonalizable()
False
sage: A.is_diagonalizable(QQbar)                                            # needs sage.rings.number_field
True
sage: D, P = A.diagonalization(QQbar)                                       # needs sage.rings.number_field
sage: P^-1*A*P == D                                                         # needs sage.rings.number_field
True
>>> from sage.all import *
>>> # needs sage.libs.pari
>>> A = matrix(QQ, Integer(4), [-Integer(4), Integer(6), Integer(3), Integer(3), -Integer(3), Integer(5), Integer(3), Integer(3), Integer(3), -Integer(6), -Integer(4), -Integer(3), -Integer(3), Integer(6), Integer(3), Integer(2)])
>>> A
[-4  6  3  3]
[-3  5  3  3]
[ 3 -6 -4 -3]
[-3  6  3  2]
>>> A.is_diagonalizable()
True
>>> A.diagonalization()
(
[ 2  0  0  0]  [ 1  1  0  0]
[ 0 -1  0  0]  [ 1  0  1  0]
[ 0  0 -1  0]  [-1  0  0  1]
[ 0  0  0 -1], [ 1  1 -2 -1]
)
>>> D, P = A.diagonalization()
>>> P**-Integer(1)*A*P == D
True

>>> # needs sage.libs.pari
>>> A = matrix(QQ, Integer(2), [Integer(0), Integer(2), Integer(1), Integer(0)])
>>> A.is_diagonalizable()
False
>>> A.is_diagonalizable(QQbar)                                            # needs sage.rings.number_field
True
>>> D, P = A.diagonalization(QQbar)                                       # needs sage.rings.number_field
>>> P**-Integer(1)*A*P == D                                                         # needs sage.rings.number_field
True

Matrices may fail to be diagonalizable for various reasons:

sage: A = matrix(QQ, 2, [1,2,3, 4,5,6]); A
[1 2 3]
[4 5 6]
sage: A.diagonalization()
Traceback (most recent call last):
...
TypeError: not a square matrix

sage: B = matrix(ZZ, 2, [1, 2, 3, 4]); B
[1 2]
[3 4]
sage: B.diagonalization()
Traceback (most recent call last):
...
ValueError: matrix entries must be from a field

sage: C = matrix(RR, 2, [1., 2., 3., 4.]); C
[1.00000000000000 2.00000000000000]
[3.00000000000000 4.00000000000000]
sage: C.diagonalization()
Traceback (most recent call last):
...
ValueError: base field must be exact,
but Real Field with 53 bits of precision is not

sage: D = matrix(QQ, 2, [0, 2, 1, 0]); D
[0 2]
[1 0]
sage: D.diagonalization()                                                   # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: not diagonalizable over Rational Field

sage: E = matrix(QQ, 2, [3, 1, 0, 3]); E
[3 1]
[0 3]
sage: E.diagonalization()                                                   # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: not diagonalizable
sage: E.jordan_form()                                                       # needs sage.combinat sage.libs.pari
[3 1]
[0 3]
>>> 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.diagonalization()
Traceback (most recent call last):
...
TypeError: not a square matrix

>>> B = matrix(ZZ, Integer(2), [Integer(1), Integer(2), Integer(3), Integer(4)]); B
[1 2]
[3 4]
>>> B.diagonalization()
Traceback (most recent call last):
...
ValueError: matrix entries must be from a field

>>> C = matrix(RR, Integer(2), [RealNumber('1.'), RealNumber('2.'), RealNumber('3.'), RealNumber('4.')]); C
[1.00000000000000 2.00000000000000]
[3.00000000000000 4.00000000000000]
>>> C.diagonalization()
Traceback (most recent call last):
...
ValueError: base field must be exact,
but Real Field with 53 bits of precision is not

>>> D = matrix(QQ, Integer(2), [Integer(0), Integer(2), Integer(1), Integer(0)]); D
[0 2]
[1 0]
>>> D.diagonalization()                                                   # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: not diagonalizable over Rational Field

>>> E = matrix(QQ, Integer(2), [Integer(3), Integer(1), Integer(0), Integer(3)]); E
[3 1]
[0 3]
>>> E.diagonalization()                                                   # needs sage.libs.pari
Traceback (most recent call last):
...
ValueError: not diagonalizable
>>> E.jordan_form()                                                       # needs sage.combinat sage.libs.pari
[3 1]
[0 3]
echelon_form(algorithm='default', cutoff=0, **kwds)[source]

Return the echelon form of self.

Note

This row reduction does not use division if the matrix is not over a field (e.g., if the matrix is over the integers). If you want to calculate the echelon form using division, then use rref(), which assumes that the matrix entries are in a field (specifically, the field of fractions of the base ring of the matrix).

INPUT:

  • algorithm – string. Which algorithm to use. Choices are

    • 'default': Let Sage choose an algorithm (default).

    • 'classical': Gauss elimination.

    • 'partial_pivoting': Gauss elimination, using partial pivoting (if base ring has absolute value)

    • 'scaled_partial_pivoting' – Gauss elimination, using scaled partial pivoting (if base ring has absolute value)

    • 'scaled_partial_pivoting_valuation': Gauss elimination, using scaled partial pivoting (if base ring has valuation)

    • 'strassen': use a Strassen divide and conquer algorithm (if available)

  • cutoff – integer; only used if the Strassen algorithm is selected

  • transformation – boolean; whether to also return the transformation matrix. Some matrix backends do not provide this information, in which case this option is ignored.

OUTPUT:

The reduced row echelon form of self, as an immutable matrix. Note that self is not changed by this command. Use echelonize() to change self in place.

If the optional parameter transformation=True is specified, the output consists of a pair \((E,T)\) of matrices where \(E\) is the echelon form of self and \(T\) is the transformation matrix.

EXAMPLES:

sage: MS = MatrixSpace(GF(19), 2, 3)
sage: C = MS.matrix([1,2,3,4,5,6])
sage: C.rank()
2
sage: C.nullity()
0
sage: C.echelon_form()
[ 1  0 18]
[ 0  1  2]
>>> from sage.all import *
>>> MS = MatrixSpace(GF(Integer(19)), Integer(2), Integer(3))
>>> C = MS.matrix([Integer(1),Integer(2),Integer(3),Integer(4),Integer(5),Integer(6)])
>>> C.rank()
2
>>> C.nullity()
0
>>> C.echelon_form()
[ 1  0 18]
[ 0  1  2]

The matrix library used for \(\ZZ/p\)-matrices does not return the transformation matrix, so the transformation option is ignored:

sage: C.echelon_form(transformation=True)
[ 1  0 18]
[ 0  1  2]

sage: D = matrix(ZZ, 2, 3, [1,2,3,4,5,6])
sage: D.echelon_form(transformation=True)
(
[1 2 3]  [ 1  0]
[0 3 6], [ 4 -1]
)
sage: E, T = D.echelon_form(transformation=True)
sage: T*D == E
True
>>> from sage.all import *
>>> C.echelon_form(transformation=True)
[ 1  0 18]
[ 0  1  2]

>>> D = matrix(ZZ, Integer(2), Integer(3), [Integer(1),Integer(2),Integer(3),Integer(4),Integer(5),Integer(6)])
>>> D.echelon_form(transformation=True)
(
[1 2 3]  [ 1  0]
[0 3 6], [ 4 -1]
)
>>> E, T = D.echelon_form(transformation=True)
>>> T*D == E
True
echelonize(algorithm='default', cutoff=0, **kwds)[source]

Transform self into a matrix in echelon form over the same base ring as self.

Note

This row reduction does not use division if the matrix is not over a field (e.g., if the matrix is over the integers). If you want to calculate the echelon form using division, then use rref(), which assumes that the matrix entries are in a field (specifically, the field of fractions of the base ring of the matrix).

INPUT:

  • algorithm – string. Which algorithm to use. Choices are

    • 'default': Let Sage choose an algorithm (default).

    • 'classical': Gauss elimination.

    • 'partial_pivoting': Gauss elimination, using partial pivoting (if base ring has absolute value)

    • 'scaled_partial_pivoting': Gauss elimination, using scaled partial pivoting (if base ring has absolute value)

    • 'scaled_partial_pivoting_valuation': Gauss elimination, using scaled partial pivoting (if base ring has valuation)

    • 'strassen': use a Strassen divide and conquer algorithm (if available)

  • cutoff – integer; only used if the Strassen algorithm is selected

  • transformation – boolean; whether to also return the transformation matrix. Some matrix backends do not provide this information, in which case this option is ignored.

OUTPUT:

The matrix self is put into echelon form. Nothing is returned unless the keyword option transformation=True is specified, in which case the transformation matrix is returned.

EXAMPLES:

sage: a = matrix(QQ, 3,3, range(9)); a
[0 1 2]
[3 4 5]
[6 7 8]
sage: a.echelonize()
sage: a
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]
>>> from sage.all import *
>>> a = matrix(QQ, Integer(3),Integer(3), range(Integer(9))); a
[0 1 2]
[3 4 5]
[6 7 8]
>>> a.echelonize()
>>> a
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]

An immutable matrix cannot be transformed into echelon form. Use self.echelon_form() instead:

sage: a = matrix(QQ, 3,3, range(9)); a.set_immutable()
sage: a.echelonize()
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).
sage: a.echelon_form()
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]
>>> from sage.all import *
>>> a = matrix(QQ, Integer(3),Integer(3), range(Integer(9))); a.set_immutable()
>>> a.echelonize()
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.echelon_form()
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]

Echelon form over the integers is what is also classically often known as Hermite normal form:

sage: a = matrix(ZZ, 3,3, range(9))
sage: a.echelonize(); a
[ 3  0 -3]
[ 0  1  2]
[ 0  0  0]
>>> from sage.all import *
>>> a = matrix(ZZ, Integer(3),Integer(3), range(Integer(9)))
>>> a.echelonize(); a
[ 3  0 -3]
[ 0  1  2]
[ 0  0  0]

We compute an echelon form both over a domain and fraction field:

sage: R.<x,y> = QQ[]
sage: a = matrix(R, 2, [x,y, x,y])
sage: a.echelon_form()               # not very useful? -- why two copies of the same row?                  # needs sage.rings.function_field
[x y]
[x y]
>>> from sage.all import *
>>> R = QQ['x, y']; (x, y,) = R._first_ngens(2)
>>> a = matrix(R, Integer(2), [x,y, x,y])
>>> a.echelon_form()               # not very useful? -- why two copies of the same row?                  # needs sage.rings.function_field
[x y]
[x y]

sage: b = a.change_ring(R.fraction_field())
sage: b.echelon_form()               # potentially useful
[  1 y/x]
[  0   0]
>>> from sage.all import *
>>> b = a.change_ring(R.fraction_field())
>>> b.echelon_form()               # potentially useful
[  1 y/x]
[  0   0]

We check that the echelon form works for matrices over \(p\)-adics. See Issue #17272:

sage: # needs sage.rings.padics
sage: R = ZpCA(5,5,print_mode='val-unit')
sage: A = matrix(R, 3,3, [250,2369,1147, 106,927,362, 90,398,2483])
sage: A
[5^3 * 2 + O(5^5)    2369 + O(5^5)    1147 + O(5^5)]
[    106 + O(5^5)     927 + O(5^5)     362 + O(5^5)]
[ 5 * 18 + O(5^5)     398 + O(5^5)    2483 + O(5^5)]
sage: K = R.fraction_field()
sage: A.change_ring(K).augment(identity_matrix(K,3)).echelon_form()
[      1 + O(5^5)           O(5^5)           O(5^5) 5 * 212 + O(5^5)    3031 + O(5^5)    2201 + O(5^5)]
[          O(5^5)       1 + O(5^5)           O(5^5)    1348 + O(5^5) 5 * 306 + O(5^5)    2648 + O(5^5)]
[          O(5^5)           O(5^5)       1 + O(5^5)    1987 + O(5^5) 5 * 263 + O(5^5)     154 + O(5^5)]
>>> from sage.all import *
>>> # needs sage.rings.padics
>>> R = ZpCA(Integer(5),Integer(5),print_mode='val-unit')
>>> A = matrix(R, Integer(3),Integer(3), [Integer(250),Integer(2369),Integer(1147), Integer(106),Integer(927),Integer(362), Integer(90),Integer(398),Integer(2483)])
>>> A
[5^3 * 2 + O(5^5)    2369 + O(5^5)    1147 + O(5^5)]
[    106 + O(5^5)     927 + O(5^5)     362 + O(5^5)]
[ 5 * 18 + O(5^5)     398 + O(5^5)    2483 + O(5^5)]
>>> K = R.fraction_field()
>>> A.change_ring(K).augment(identity_matrix(K,Integer(3))).echelon_form()
[      1 + O(5^5)           O(5^5)           O(5^5) 5 * 212 + O(5^5)    3031 + O(5^5)    2201 + O(5^5)]
[          O(5^5)       1 + O(5^5)           O(5^5)    1348 + O(5^5) 5 * 306 + O(5^5)    2648 + O(5^5)]
[          O(5^5)           O(5^5)       1 + O(5^5)    1987 + O(5^5) 5 * 263 + O(5^5)     154 + O(5^5)]

Echelon form is not defined over arbitrary rings:

sage: a = matrix(Integers(9), 3,3, range(9))
sage: a.echelon_form()
Traceback (most recent call last):
...
NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 9'.
>>> from sage.all import *
>>> a = matrix(Integers(Integer(9)), Integer(3),Integer(3), range(Integer(9)))
>>> a.echelon_form()
Traceback (most recent call last):
...
NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 9'.

Involving a sparse matrix:

sage: m = matrix(3,[1, 1, 1, 1, 0, 2, 1, 2, 0], sparse=True); m
[1 1 1]
[1 0 2]
[1 2 0]
sage: m.echelon_form()
[ 1  0  2]
[ 0  1 -1]
[ 0  0  0]
sage: m.echelonize(); m
[ 1  0  2]
[ 0  1 -1]
[ 0  0  0]
>>> from sage.all import *
>>> m = matrix(Integer(3),[Integer(1), Integer(1), Integer(1), Integer(1), Integer(0), Integer(2), Integer(1), Integer(2), Integer(0)], sparse=True); m
[1 1 1]
[1 0 2]
[1 2 0]
>>> m.echelon_form()
[ 1  0  2]
[ 0  1 -1]
[ 0  0  0]
>>> m.echelonize(); m
[ 1  0  2]
[ 0  1 -1]
[ 0  0  0]

The transformation matrix is optionally returned:

sage: m_original = m
sage: transformation_matrix = m.echelonize(transformation=True)
sage: m == transformation_matrix * m_original
True
>>> from sage.all import *
>>> m_original = m
>>> transformation_matrix = m.echelonize(transformation=True)
>>> m == transformation_matrix * m_original
True
eigenmatrix_left(other=None)[source]

Return matrices \(D\) and \(P\), where \(D\) is a diagonal matrix of eigenvalues and the rows of \(P\) are corresponding eigenvectors (or zero vectors).

INPUT:

  • other – a square matrix \(B\) (default: None) in a generalized eigenvalue problem; if None, an ordinary eigenvalue problem is solved

OUTPUT:

If self is a square matrix \(A\), then the output is a diagonal matrix \(D\) and a matrix \(P\) such that

\[P A = D P,\]

where the rows of \(P\) are eigenvectors of \(A\) and the diagonal entries of \(D\) are the corresponding eigenvalues.

If a matrix \(B\) is passed as optional argument, the output is a solution to the generalized eigenvalue problem such that

\[P A = D P B.\]

The ordinary eigenvalue problem is equivalent to the generalized one if \(B\) is the identity matrix.

The generalized eigenvector decomposition is currently only implemented for matrices over RDF and CDF.

EXAMPLES:

sage: # needs sage.rings.number_field
sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: D, P = A.eigenmatrix_left()
sage: D
[                  0                   0                   0]
[                  0 -1.348469228349535?                   0]
[                  0                   0  13.34846922834954?]
sage: P
[                   1                   -2                    1]
[                   1  0.3101020514433644? -0.3797958971132713?]
[                   1   1.289897948556636?   1.579795897113272?]
sage: P*A == D*P
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> D, P = A.eigenmatrix_left()
>>> D
[                  0                   0                   0]
[                  0 -1.348469228349535?                   0]
[                  0                   0  13.34846922834954?]
>>> P
[                   1                   -2                    1]
[                   1  0.3101020514433644? -0.3797958971132713?]
[                   1   1.289897948556636?   1.579795897113272?]
>>> P*A == D*P
True

Because \(P\) is invertible, \(A\) is diagonalizable.

sage: A == (~P)*D*P                                                         # needs sage.rings.number_field
True
>>> from sage.all import *
>>> A == (~P)*D*P                                                         # needs sage.rings.number_field
True

The matrix \(P\) may contain zero rows corresponding to eigenvalues for which the algebraic multiplicity is greater than the geometric multiplicity. In these cases, the matrix is not diagonalizable.

sage: # needs sage.rings.number_field
sage: A = jordan_block(2, 3); A
[2 1 0]
[0 2 1]
[0 0 2]
sage: D, P = A.eigenmatrix_left()
sage: D
[2 0 0]
[0 2 0]
[0 0 2]
sage: P
[0 0 1]
[0 0 0]
[0 0 0]
sage: P*A == D*P
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = jordan_block(Integer(2), Integer(3)); A
[2 1 0]
[0 2 1]
[0 0 2]
>>> D, P = A.eigenmatrix_left()
>>> D
[2 0 0]
[0 2 0]
[0 0 2]
>>> P
[0 0 1]
[0 0 0]
[0 0 0]
>>> P*A == D*P
True

A generalized eigenvector decomposition:

sage: # needs scipy
sage: A = matrix(RDF, [[1, -2], [3, 4]])
sage: B = matrix(RDF, [[0, 7], [2, -3]])
sage: D, P = A.eigenmatrix_left(B)
sage: (P * A - D * P * B).norm() < 1e-14
True
>>> from sage.all import *
>>> # needs scipy
>>> A = matrix(RDF, [[Integer(1), -Integer(2)], [Integer(3), Integer(4)]])
>>> B = matrix(RDF, [[Integer(0), Integer(7)], [Integer(2), -Integer(3)]])
>>> D, P = A.eigenmatrix_left(B)
>>> (P * A - D * P * B).norm() < RealNumber('1e-14')
True

The matrix \(B\) in a generalized eigenvalue problem may be singular:

sage: # needs scipy sage.rings.complex_double sage.symbolic
sage: A = matrix.identity(CDF, 2)
sage: B = matrix(CDF, [[2, 1+I], [4, 2+2*I]])
sage: D, P = A.eigenmatrix_left(B)
sage: D.diagonal()  # tol 1e-14
[0.2 - 0.1*I, +infinity]
>>> from sage.all import *
>>> # needs scipy sage.rings.complex_double sage.symbolic
>>> A = matrix.identity(CDF, Integer(2))
>>> B = matrix(CDF, [[Integer(2), Integer(1)+I], [Integer(4), Integer(2)+Integer(2)*I]])
>>> D, P = A.eigenmatrix_left(B)
>>> D.diagonal()  # tol 1e-14
[0.2 - 0.1*I, +infinity]

In this case, we can still verify the eigenvector equation for the first eigenvalue and first eigenvector:

sage: # needs scipy sage.rings.complex_double sage.symbolic
sage: l = D[0, 0]
sage: v = P[0, :]
sage: (v * A - l * v * B).norm() < 1e-14
True
>>> from sage.all import *
>>> # needs scipy sage.rings.complex_double sage.symbolic
>>> l = D[Integer(0), Integer(0)]
>>> v = P[Integer(0), :]
>>> (v * A - l * v * B).norm() < RealNumber('1e-14')
True

The second eigenvector is contained in the left kernel of \(B\):

sage: (P[1, :] * B).norm() < 1e-14                                          # needs scipy sage.rings.complex_double sage.symbolic
True
>>> from sage.all import *
>>> (P[Integer(1), :] * B).norm() < RealNumber('1e-14')                                          # needs scipy sage.rings.complex_double sage.symbolic
True
eigenmatrix_right(other=None)[source]

Return matrices \(D\) and \(P\), where \(D\) is a diagonal matrix of eigenvalues and the columns of \(P\) are corresponding eigenvectors (or zero vectors).

INPUT:

  • other – a square matrix \(B\) (default: None) in a generalized eigenvalue problem; if None, an ordinary eigenvalue problem is solved

OUTPUT:

If self is a square matrix \(A\), then the output is a diagonal matrix \(D\) and a matrix \(P\) such that

\[A P = P D,\]

where the columns of \(P\) are eigenvectors of \(A\) and the diagonal entries of \(D\) are the corresponding eigenvalues.

If a matrix \(B\) is passed as optional argument, the output is a solution to the generalized eigenvalue problem such that

\[A P = B P D.\]

The ordinary eigenvalue problem is equivalent to the generalized one if \(B\) is the identity matrix.

The generalized eigenvector decomposition is currently only implemented for matrices over RDF and CDF.

EXAMPLES:

sage: # needs sage.rings.number_field
sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: D, P = A.eigenmatrix_right()
sage: D
[                  0                   0                   0]
[                  0 -1.348469228349535?                   0]
[                  0                   0  13.34846922834954?]
sage: P
[                   1                    1                    1]
[                  -2  0.1303061543300932?   3.069693845669907?]
[                   1 -0.7393876913398137?   5.139387691339814?]
sage: A*P == P*D
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> D, P = A.eigenmatrix_right()
>>> D
[                  0                   0                   0]
[                  0 -1.348469228349535?                   0]
[                  0                   0  13.34846922834954?]
>>> P
[                   1                    1                    1]
[                  -2  0.1303061543300932?   3.069693845669907?]
[                   1 -0.7393876913398137?   5.139387691339814?]
>>> A*P == P*D
True

Because \(P\) is invertible, \(A\) is diagonalizable.

sage: A == P*D*(~P)                                                         # needs sage.rings.number_field
True
>>> from sage.all import *
>>> A == P*D*(~P)                                                         # needs sage.rings.number_field
True

The matrix \(P\) may contain zero columns corresponding to eigenvalues for which the algebraic multiplicity is greater than the geometric multiplicity. In these cases, the matrix is not diagonalizable.

sage: # needs sage.rings.number_field
sage: A = jordan_block(2, 3); A
[2 1 0]
[0 2 1]
[0 0 2]
sage: D, P = A.eigenmatrix_right()
sage: D
[2 0 0]
[0 2 0]
[0 0 2]
sage: P
[1 0 0]
[0 0 0]
[0 0 0]
sage: A*P == P*D
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = jordan_block(Integer(2), Integer(3)); A
[2 1 0]
[0 2 1]
[0 0 2]
>>> D, P = A.eigenmatrix_right()
>>> D
[2 0 0]
[0 2 0]
[0 0 2]
>>> P
[1 0 0]
[0 0 0]
[0 0 0]
>>> A*P == P*D
True

A generalized eigenvector decomposition:

sage: # needs scipy
sage: A = matrix(RDF, [[1, -2], [3, 4]])
sage: B = matrix(RDF, [[0, 7], [2, -3]])
sage: D, P = A.eigenmatrix_right(B)
sage: (A * P - B * P * D).norm() < 1e-14
True
>>> from sage.all import *
>>> # needs scipy
>>> A = matrix(RDF, [[Integer(1), -Integer(2)], [Integer(3), Integer(4)]])
>>> B = matrix(RDF, [[Integer(0), Integer(7)], [Integer(2), -Integer(3)]])
>>> D, P = A.eigenmatrix_right(B)
>>> (A * P - B * P * D).norm() < RealNumber('1e-14')
True

The matrix \(B\) in a generalized eigenvalue problem may be singular:

sage: # needs scipy
sage: A = matrix.identity(RDF, 2)
sage: B = matrix(RDF, [[3, 5], [6, 10]])
sage: D, P = A.eigenmatrix_right(B); D   # tol 1e-14
[0.07692307692307694                 0.0]
[                0.0           +infinity]
>>> from sage.all import *
>>> # needs scipy
>>> A = matrix.identity(RDF, Integer(2))
>>> B = matrix(RDF, [[Integer(3), Integer(5)], [Integer(6), Integer(10)]])
>>> D, P = A.eigenmatrix_right(B); D   # tol 1e-14
[0.07692307692307694                 0.0]
[                0.0           +infinity]

In this case, we can still verify the eigenvector equation for the first eigenvalue and first eigenvector:

sage: # needs scipy
sage: l = D[0, 0]
sage: v = P[:, 0]
sage: (A * v  - B * v * l).norm() < 1e-14
True
>>> from sage.all import *
>>> # needs scipy
>>> l = D[Integer(0), Integer(0)]
>>> v = P[:, Integer(0)]
>>> (A * v  - B * v * l).norm() < RealNumber('1e-14')
True

The second eigenvector is contained in the right kernel of \(B\):

sage: (B * P[:, 1]).norm() < 1e-14                                          # needs scipy
True
>>> from sage.all import *
>>> (B * P[:, Integer(1)]).norm() < RealNumber('1e-14')                                          # needs scipy
True
eigenspaces_left(format='all', var='a', algebraic_multiplicity=False)[source]

Compute the left eigenspaces of a matrix.

Note that eigenspaces_left() and left_eigenspaces() are identical methods. Here “left” refers to the eigenvectors being placed to the left of the matrix.

INPUT:

  • self – a square matrix over an exact field. For inexact matrices consult the numerical or symbolic matrix classes

  • format – one of:

    • 'all' – attempts to create every eigenspace. This will always be possible for matrices with rational entries

    • 'galois' – for each irreducible factor of the characteristic polynomial, a single eigenspace will be output for a single root/eigenvalue for the irreducible factor

    • None – default; uses the 'all' format if the base ring is contained in an algebraically closed field which is implemented. Otherwise, uses the 'galois' format.

  • var – string (default: 'a'); variable name used to represent elements of the root field of each irreducible factor of the characteristic polynomial. If var='a', then the root fields will be in terms of a0, a1, a2, ..., where the numbering runs across all the irreducible factors of the characteristic polynomial, even for linear factors.

  • algebraic_multiplicity – boolean (default: False); whether to include the algebraic multiplicity of each eigenvalue in the output. See the discussion below.

OUTPUT:

If algebraic_multiplicity=False, return a list of pairs \((e, V)\) where \(e\) is an eigenvalue of the matrix, and \(V\) is the corresponding left eigenspace. For Galois conjugates of eigenvalues, there may be just one representative eigenspace, depending on the format keyword.

If algebraic_multiplicity=True, return a list of triples \((e, V, n)\) where \(e\) and \(V\) are as above and \(n\) is the algebraic multiplicity of the eigenvalue.

Warning

Uses a somewhat naive algorithm (simply factors the characteristic polynomial and computes kernels directly over the extension field).

EXAMPLES:

We compute the left eigenspaces of a \(3\times 3\) rational matrix. First, we request 'all' of the eigenvalues, so the results are in the field of algebraic numbers, QQbar. Then we request just one eigenspace per irreducible factor of the characteristic polynomial with format='galois'.

sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: es = A.eigenspaces_left(format='all'); es                             # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (-1.348469228349535?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                   1  0.3101020514433644? -0.3797958971132713?]),
  (13.34846922834954?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                 1 1.289897948556636? 1.579795897113272?]) ]

sage: # needs sage.rings.number_field
sage: es = A.eigenspaces_left(format='galois'); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5]) ]
sage: es = A.eigenspaces_left(format='galois',
....:                         algebraic_multiplicity=True); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1],
   1),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5],
   1) ]
sage: e, v, n = es[0]; v = v.basis()[0]
sage: delta = e*v - v*A
sage: abs(abs(delta)) < 1e-10
True
>>> from sage.all import *
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> es = A.eigenspaces_left(format='all'); es                             # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (-1.348469228349535?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                   1  0.3101020514433644? -0.3797958971132713?]),
  (13.34846922834954?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                 1 1.289897948556636? 1.579795897113272?]) ]

>>> # needs sage.rings.number_field
>>> es = A.eigenspaces_left(format='galois'); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5]) ]
>>> es = A.eigenspaces_left(format='galois',
...                         algebraic_multiplicity=True); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1],
   1),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5],
   1) ]
>>> e, v, n = es[Integer(0)]; v = v.basis()[Integer(0)]
>>> delta = e*v - v*A
>>> abs(abs(delta)) < RealNumber('1e-10')
True

The same computation, but with implicit base change to a field.

sage: A = matrix(ZZ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: A.eigenspaces_left(format='galois')                                   # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5]) ]
>>> from sage.all import *
>>> A = matrix(ZZ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> A.eigenspaces_left(format='galois')                                   # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [            1 1/15*a1 + 2/5 2/15*a1 - 1/5]) ]

We compute the left eigenspaces of the matrix of the Hecke operator \(T_2\) on level 43 modular symbols, both with all eigenvalues (the default) and with one subspace per factor.

sage: # needs sage.modular
sage: A = ModularSymbols(43).T(2).matrix(); A
[ 3  0  0  0  0  0 -1]
[ 0 -2  1  0  0  0  0]
[ 0 -1  1  1  0 -1  0]
[ 0 -1  0 -1  2 -1  1]
[ 0 -1  0  1  1 -1  1]
[ 0  0 -2  0  2 -2  1]
[ 0  0 -1  0  1  0 -1]
sage: A.base_ring()
Rational Field
sage: f = A.charpoly(); f
x^7 + x^6 - 12*x^5 - 16*x^4 + 36*x^3 + 52*x^2 - 32*x - 48
sage: factor(f)
(x - 3) * (x + 2)^2 * (x^2 - 2)^2
sage: A.eigenspaces_left(algebraic_multiplicity=True)
[ (3,
   Vector space of degree 7 and dimension 1 over Rational Field
     User basis matrix:
     [   1    0  1/7    0 -1/7    0 -2/7],
   1),
  (-2,
   Vector space of degree 7 and dimension 2 over Rational Field
     User basis matrix:
     [ 0  1  0  1 -1  1 -1]
     [ 0  0  1  0 -1  2 -1],
   2),
  (-1.414213562373095?,
   Vector space of degree 7 and dimension 2 over Algebraic Field
     User basis matrix:
     [                  0                   1                   0                  -1 0.4142135623730951?                   1                  -1]
     [                  0                   0                   1                   0                  -1                   0  2.414213562373095?],
   2),
  (1.414213562373095?,
   Vector space of degree 7 and dimension 2 over Algebraic Field
     User basis matrix:
     [                   0                    1                    0                   -1  -2.414213562373095?                    1                   -1]
     [                   0                    0                    1                    0                   -1                    0 -0.4142135623730951?],
   2) ]
sage: A.eigenspaces_left(format='galois', algebraic_multiplicity=True)
[ (3,
   Vector space of degree 7 and dimension 1 over Rational Field
     User basis matrix:
     [   1    0  1/7    0 -1/7    0 -2/7],
   1),
  (-2,
   Vector space of degree 7 and dimension 2 over Rational Field
     User basis matrix:
     [ 0  1  0  1 -1  1 -1]
     [ 0  0  1  0 -1  2 -1],
   2),
  (a2,
   Vector space of degree 7 and dimension 2
    over Number Field in a2 with defining polynomial x^2 - 2
     User basis matrix:
     [      0       1       0      -1 -a2 - 1       1      -1]
     [      0       0       1       0      -1       0 -a2 + 1],
   2) ]
>>> from sage.all import *
>>> # needs sage.modular
>>> A = ModularSymbols(Integer(43)).T(Integer(2)).matrix(); A
[ 3  0  0  0  0  0 -1]
[ 0 -2  1  0  0  0  0]
[ 0 -1  1  1  0 -1  0]
[ 0 -1  0 -1  2 -1  1]
[ 0 -1  0  1  1 -1  1]
[ 0  0 -2  0  2 -2  1]
[ 0  0 -1  0  1  0 -1]
>>> A.base_ring()
Rational Field
>>> f = A.charpoly(); f
x^7 + x^6 - 12*x^5 - 16*x^4 + 36*x^3 + 52*x^2 - 32*x - 48
>>> factor(f)
(x - 3) * (x + 2)^2 * (x^2 - 2)^2
>>> A.eigenspaces_left(algebraic_multiplicity=True)
[ (3,
   Vector space of degree 7 and dimension 1 over Rational Field
     User basis matrix:
     [   1    0  1/7    0 -1/7    0 -2/7],
   1),
  (-2,
   Vector space of degree 7 and dimension 2 over Rational Field
     User basis matrix:
     [ 0  1  0  1 -1  1 -1]
     [ 0  0  1  0 -1  2 -1],
   2),
  (-1.414213562373095?,
   Vector space of degree 7 and dimension 2 over Algebraic Field
     User basis matrix:
     [                  0                   1                   0                  -1 0.4142135623730951?                   1                  -1]
     [                  0                   0                   1                   0                  -1                   0  2.414213562373095?],
   2),
  (1.414213562373095?,
   Vector space of degree 7 and dimension 2 over Algebraic Field
     User basis matrix:
     [                   0                    1                    0                   -1  -2.414213562373095?                    1                   -1]
     [                   0                    0                    1                    0                   -1                    0 -0.4142135623730951?],
   2) ]
>>> A.eigenspaces_left(format='galois', algebraic_multiplicity=True)
[ (3,
   Vector space of degree 7 and dimension 1 over Rational Field
     User basis matrix:
     [   1    0  1/7    0 -1/7    0 -2/7],
   1),
  (-2,
   Vector space of degree 7 and dimension 2 over Rational Field
     User basis matrix:
     [ 0  1  0  1 -1  1 -1]
     [ 0  0  1  0 -1  2 -1],
   2),
  (a2,
   Vector space of degree 7 and dimension 2
    over Number Field in a2 with defining polynomial x^2 - 2
     User basis matrix:
     [      0       1       0      -1 -a2 - 1       1      -1]
     [      0       0       1       0      -1       0 -a2 + 1],
   2) ]

Next we compute the left eigenspaces over the finite field of order 11.

sage: # needs sage.modular sage.rings.finite_rings
sage: A = ModularSymbols(43, base_ring=GF(11), sign=1).T(2).matrix(); A
[ 3  0  9  0]
[ 0  9  0 10]
[ 0  0 10  1]
[ 0  0  1  1]
sage: A.base_ring()
Finite Field of size 11
sage: A.charpoly()
x^4 + 10*x^3 + 3*x^2 + 2*x + 1
sage: A.eigenspaces_left(format='galois', var='beta')
[ (9,
   Vector space of degree 4 and dimension 1 over Finite Field of size 11
     User basis matrix: [0 1 5 6]),
 (3, Vector space of degree 4 and dimension 1 over Finite Field of size 11
     User basis matrix: [1 0 1 6]),
 (beta2, Vector space of degree 4 and dimension 1
          over Univariate Quotient Polynomial Ring in beta2
           over Finite Field of size 11 with modulus x^2 + 9
         User basis matrix: [        0         0         1 beta2 + 1])
]
>>> from sage.all import *
>>> # needs sage.modular sage.rings.finite_rings
>>> A = ModularSymbols(Integer(43), base_ring=GF(Integer(11)), sign=Integer(1)).T(Integer(2)).matrix(); A
[ 3  0  9  0]
[ 0  9  0 10]
[ 0  0 10  1]
[ 0  0  1  1]
>>> A.base_ring()
Finite Field of size 11
>>> A.charpoly()
x^4 + 10*x^3 + 3*x^2 + 2*x + 1
>>> A.eigenspaces_left(format='galois', var='beta')
[ (9,
   Vector space of degree 4 and dimension 1 over Finite Field of size 11
     User basis matrix: [0 1 5 6]),
 (3, Vector space of degree 4 and dimension 1 over Finite Field of size 11
     User basis matrix: [1 0 1 6]),
 (beta2, Vector space of degree 4 and dimension 1
          over Univariate Quotient Polynomial Ring in beta2
           over Finite Field of size 11 with modulus x^2 + 9
         User basis matrix: [        0         0         1 beta2 + 1])
]

This method is only applicable to exact matrices. The “eigenmatrix” routines for matrices with double-precision floating-point entries (RDF, CDF) are the best alternative. (Since some platforms return eigenvectors that are the negatives of those given here, this one example is not tested here.) There are also “eigenmatrix” routines for matrices with symbolic entries.

sage: A = matrix(QQ, 3, 3, range(9))
sage: A.change_ring(RR).eigenspaces_left()
Traceback (most recent call last):
...
NotImplementedError: eigenspaces cannot be computed reliably
for inexact rings such as Real Field with 53 bits of precision,
consult numerical or symbolic matrix classes for other options

sage: # needs scipy
sage: em = A.change_ring(RDF).eigenmatrix_left()
sage: eigenvalues = em[0]; eigenvalues.dense_matrix()  # abs tol 1e-13
[13.348469228349522                0.0                 0.0]
[               0.0 -1.348469228349534                 0.0]
[               0.0                0.0                 0.0]
sage: eigenvectors = em[1]; eigenvectors # not tested
[ 0.440242867...  0.567868371...  0.695493875...]
[ 0.897878732...  0.278434036... -0.341010658...]
[ 0.408248290... -0.816496580...  0.408248290...]

sage: # needs sage.symbolic
sage: x, y = var('x y')
sage: S = matrix([[x, y], [y, 3*x^2]])
sage: em = S.eigenmatrix_left()
sage: eigenvalues = em[0]; eigenvalues
[3/2*x^2 + 1/2*x - 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)                                                       0]
[                                                      0 3/2*x^2 + 1/2*x + 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)]
sage: eigenvectors = em[1]; eigenvectors
[                                                    1 1/2*(3*x^2 - x - sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]
[                                                    1 1/2*(3*x^2 - x + sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]
>>> from sage.all import *
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9)))
>>> A.change_ring(RR).eigenspaces_left()
Traceback (most recent call last):
...
NotImplementedError: eigenspaces cannot be computed reliably
for inexact rings such as Real Field with 53 bits of precision,
consult numerical or symbolic matrix classes for other options

>>> # needs scipy
>>> em = A.change_ring(RDF).eigenmatrix_left()
>>> eigenvalues = em[Integer(0)]; eigenvalues.dense_matrix()  # abs tol 1e-13
[13.348469228349522                0.0                 0.0]
[               0.0 -1.348469228349534                 0.0]
[               0.0                0.0                 0.0]
>>> eigenvectors = em[Integer(1)]; eigenvectors # not tested
[ 0.440242867...  0.567868371...  0.695493875...]
[ 0.897878732...  0.278434036... -0.341010658...]
[ 0.408248290... -0.816496580...  0.408248290...]

>>> # needs sage.symbolic
>>> x, y = var('x y')
>>> S = matrix([[x, y], [y, Integer(3)*x**Integer(2)]])
>>> em = S.eigenmatrix_left()
>>> eigenvalues = em[Integer(0)]; eigenvalues
[3/2*x^2 + 1/2*x - 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)                                                       0]
[                                                      0 3/2*x^2 + 1/2*x + 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)]
>>> eigenvectors = em[Integer(1)]; eigenvectors
[                                                    1 1/2*(3*x^2 - x - sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]
[                                                    1 1/2*(3*x^2 - x + sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]

A request for 'all' the eigenvalues, when it is not possible, will raise an error. Using the 'galois' format option is more likely to be successful.

sage: # needs sage.rings.finite_rings
sage: F.<b> = FiniteField(11^2)
sage: A = matrix(F, [[b + 1, b + 1], [10*b + 4, 5*b + 4]])
sage: A.eigenspaces_left(format='all')                                      # needs sage.rings.number_field
Traceback (most recent call last):
...
NotImplementedError: unable to construct eigenspaces for eigenvalues outside the base field,
try the keyword option: format='galois'
sage: A.eigenspaces_left(format='galois')                                   # needs sage.rings.number_field
[ (a0,
   Vector space of degree 2 and dimension 1 over
    Univariate Quotient Polynomial Ring in a0 over
     Finite Field in b of size 11^2
    with modulus x^2 + (5*b + 6)*x + 8*b + 10
     User basis matrix:
     [               1 6*b*a0 + 3*b + 1]) ]
>>> from sage.all import *
>>> # needs sage.rings.finite_rings
>>> F = FiniteField(Integer(11)**Integer(2), names=('b',)); (b,) = F._first_ngens(1)
>>> A = matrix(F, [[b + Integer(1), b + Integer(1)], [Integer(10)*b + Integer(4), Integer(5)*b + Integer(4)]])
>>> A.eigenspaces_left(format='all')                                      # needs sage.rings.number_field
Traceback (most recent call last):
...
NotImplementedError: unable to construct eigenspaces for eigenvalues outside the base field,
try the keyword option: format='galois'
>>> A.eigenspaces_left(format='galois')                                   # needs sage.rings.number_field
[ (a0,
   Vector space of degree 2 and dimension 1 over
    Univariate Quotient Polynomial Ring in a0 over
     Finite Field in b of size 11^2
    with modulus x^2 + (5*b + 6)*x + 8*b + 10
     User basis matrix:
     [               1 6*b*a0 + 3*b + 1]) ]
eigenspaces_right(format='all', var='a', algebraic_multiplicity=False)[source]

Compute the right eigenspaces of a matrix.

Note that eigenspaces_right() and right_eigenspaces() are identical methods. Here “right” refers to the eigenvectors being placed to the right of the matrix.

INPUT:

  • self – a square matrix over an exact field. For inexact matrices consult the numerical or symbolic matrix classes

  • format – (default: None)

    • 'all' – attempts to create every eigenspace. This will always be possible for matrices with rational entries

    • 'galois' – for each irreducible factor of the characteristic polynomial, a single eigenspace will be output for a single root/eigenvalue for the irreducible factor

    • None – uses the ‘all’ format if the base ring is contained in an algebraically closed field which is implemented. Otherwise, uses the ‘galois’ format.

  • var – (default: 'a') variable name used to represent elements of the root field of each irreducible factor of the characteristic polynomial. If var=’a’, then the root fields will be in terms of a0, a1, a2, …., where the numbering runs across all the irreducible factors of the characteristic polynomial, even for linear factors.

  • algebraic_multiplicity – (default: False) whether or not to include the algebraic multiplicity of each eigenvalue in the output. See the discussion below.

OUTPUT:

If algebraic_multiplicity=False, return a list of pairs (e, V) where e is an eigenvalue of the matrix, and V is the corresponding left eigenspace. For Galois conjugates of eigenvalues, there may be just one representative eigenspace, depending on the format keyword.

If algebraic_multiplicity=True, return a list of triples (e, V, n) where e and V are as above and n is the algebraic multiplicity of the eigenvalue.

Warning

Uses a somewhat naive algorithm (simply factors the characteristic polynomial and computes kernels directly over the extension field).

EXAMPLES:

Right eigenspaces are computed from the left eigenspaces of the transpose of the matrix. As such, there is a greater collection of illustrative examples at the eigenspaces_left().

We compute the right eigenspaces of a \(3\times 3\) rational matrix.

sage: # needs sage.rings.number_field
sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: A.eigenspaces_right()
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (-1.348469228349535?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                   1  0.1303061543300932? -0.7393876913398137?]),
  (13.34846922834954?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                 1 3.069693845669907? 5.139387691339814?]) ]
sage: es = A.eigenspaces_right(format='galois'); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5]) ]
sage: es = A.eigenspaces_right(format='galois',
....:                          algebraic_multiplicity=True); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1],
   1),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5],
   1) ]
sage: e, v, n = es[0]; v = v.basis()[0]
sage: delta = v*e - A*v
sage: abs(abs(delta)) < 1e-10
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> A.eigenspaces_right()
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (-1.348469228349535?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                   1  0.1303061543300932? -0.7393876913398137?]),
  (13.34846922834954?,
   Vector space of degree 3 and dimension 1 over Algebraic Field
     User basis matrix:
     [                 1 3.069693845669907? 5.139387691339814?]) ]
>>> es = A.eigenspaces_right(format='galois'); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5]) ]
>>> es = A.eigenspaces_right(format='galois',
...                          algebraic_multiplicity=True); es
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1],
   1),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5],
   1) ]
>>> e, v, n = es[Integer(0)]; v = v.basis()[Integer(0)]
>>> delta = v*e - A*v
>>> abs(abs(delta)) < RealNumber('1e-10')
True

The same computation, but with implicit base change to a field:

sage: A = matrix(ZZ, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: A.eigenspaces_right(format='galois')                                  # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5]) ]
>>> from sage.all import *
>>> A = matrix(ZZ, Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> A.eigenspaces_right(format='galois')                                  # needs sage.rings.number_field
[ (0,
   Vector space of degree 3 and dimension 1 over Rational Field
     User basis matrix:
     [ 1 -2  1]),
  (a1,
   Vector space of degree 3 and dimension 1 over
    Number Field in a1 with defining polynomial x^2 - 12*x - 18
     User basis matrix:
     [           1 1/5*a1 + 2/5 2/5*a1 - 1/5]) ]

This method is only applicable to exact matrices. The “eigenmatrix” routines for matrices with double-precision floating-point entries (RDF, CDF) are the best alternative. (Since some platforms return eigenvectors that are the negatives of those given here, this one example is not tested here.) There are also “eigenmatrix” routines for matrices with symbolic entries.

sage: B = matrix(RR, 3, 3, range(9))
sage: B.eigenspaces_right()
Traceback (most recent call last):
...
NotImplementedError: eigenspaces cannot be computed reliably
for inexact rings such as Real Field with 53 bits of precision,
consult numerical or symbolic matrix classes for other options

sage: # needs scipy
sage: em = B.change_ring(RDF).eigenmatrix_right()
sage: eigenvalues = em[0]; eigenvalues.dense_matrix()  # abs tol 1e-13
[13.348469228349522                0.0                0.0]
[               0.0 -1.348469228349534                0.0]
[               0.0                0.0                0.0]
sage: eigenvectors = em[1]; eigenvectors # not tested
[ 0.164763817...  0.799699663...  0.408248290...]
[ 0.505774475...  0.104205787... -0.816496580...]
[ 0.846785134... -0.591288087...  0.408248290...]

sage: # needs sage.symbolic
sage: x, y = var('x y')
sage: S = matrix([[x, y], [y, 3*x^2]])
sage: em = S.eigenmatrix_right()
sage: eigenvalues = em[0]; eigenvalues
[3/2*x^2 + 1/2*x - 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)                                                       0]
[                                                      0 3/2*x^2 + 1/2*x + 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)]
sage: eigenvectors = em[1]; eigenvectors
[                                                    1                                                     1]
[1/2*(3*x^2 - x - sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y 1/2*(3*x^2 - x + sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]
>>> from sage.all import *
>>> B = matrix(RR, Integer(3), Integer(3), range(Integer(9)))
>>> B.eigenspaces_right()
Traceback (most recent call last):
...
NotImplementedError: eigenspaces cannot be computed reliably
for inexact rings such as Real Field with 53 bits of precision,
consult numerical or symbolic matrix classes for other options

>>> # needs scipy
>>> em = B.change_ring(RDF).eigenmatrix_right()
>>> eigenvalues = em[Integer(0)]; eigenvalues.dense_matrix()  # abs tol 1e-13
[13.348469228349522                0.0                0.0]
[               0.0 -1.348469228349534                0.0]
[               0.0                0.0                0.0]
>>> eigenvectors = em[Integer(1)]; eigenvectors # not tested
[ 0.164763817...  0.799699663...  0.408248290...]
[ 0.505774475...  0.104205787... -0.816496580...]
[ 0.846785134... -0.591288087...  0.408248290...]

>>> # needs sage.symbolic
>>> x, y = var('x y')
>>> S = matrix([[x, y], [y, Integer(3)*x**Integer(2)]])
>>> em = S.eigenmatrix_right()
>>> eigenvalues = em[Integer(0)]; eigenvalues
[3/2*x^2 + 1/2*x - 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)                                                       0]
[                                                      0 3/2*x^2 + 1/2*x + 1/2*sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2)]
>>> eigenvectors = em[Integer(1)]; eigenvectors
[                                                    1                                                     1]
[1/2*(3*x^2 - x - sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y 1/2*(3*x^2 - x + sqrt(9*x^4 - 6*x^3 + x^2 + 4*y^2))/y]
eigenvalue_multiplicity(s)[source]

Return the multiplicity of s as a generalized eigenvalue of the matrix.

EXAMPLES:

sage: M = Matrix(QQ, [[0,1],[0,0]])
sage: M.eigenvalue_multiplicity(0)
2
sage: M.eigenvalue_multiplicity(1)
0

sage: M = posets.DiamondPoset(5).coxeter_transformation()                   # needs sage.graphs sage.libs.flint
sage: [M.eigenvalue_multiplicity(x) for x in [-1, 1]]                       # needs sage.graphs sage.libs.flint
[3, 2]
>>> from sage.all import *
>>> M = Matrix(QQ, [[Integer(0),Integer(1)],[Integer(0),Integer(0)]])
>>> M.eigenvalue_multiplicity(Integer(0))
2
>>> M.eigenvalue_multiplicity(Integer(1))
0

>>> M = posets.DiamondPoset(Integer(5)).coxeter_transformation()                   # needs sage.graphs sage.libs.flint
>>> [M.eigenvalue_multiplicity(x) for x in [-Integer(1), Integer(1)]]                       # needs sage.graphs sage.libs.flint
[3, 2]
eigenvalues(extend=True)[source]

Return a sequence of the eigenvalues of a matrix, with multiplicity. If the eigenvalues are roots of polynomials in QQ, then QQbar elements are returned that represent each separate root.

If the option extend is set to False, only eigenvalues in the base ring are considered.

EXAMPLES:

sage: a = matrix(ZZ, 4, range(16)); a
[ 0  1  2  3]
[ 4  5  6  7]
[ 8  9 10 11]
[12 13 14 15]
sage: sorted(a.eigenvalues(), reverse=True)                                 # needs sage.rings.number_field
[32.46424919657298?, 0, 0, -2.464249196572981?]
>>> from sage.all import *
>>> a = matrix(ZZ, Integer(4), range(Integer(16))); a
[ 0  1  2  3]
[ 4  5  6  7]
[ 8  9 10 11]
[12 13 14 15]
>>> sorted(a.eigenvalues(), reverse=True)                                 # needs sage.rings.number_field
[32.46424919657298?, 0, 0, -2.464249196572981?]

sage: a = matrix([(1, 9, -1, -1),
....:             (-2, 0, -10, 2),
....:             (-1, 0, 15, -2),
....:             (0, 1, 0, -1)])
sage: a.eigenvalues()                                                       # needs sage.rings.number_field
[-0.9386318578049146?,
 15.50655435353258?,
 0.2160387521361705? - 4.713151979747493?*I,
 0.2160387521361705? + 4.713151979747493?*I]
>>> from sage.all import *
>>> a = matrix([(Integer(1), Integer(9), -Integer(1), -Integer(1)),
...             (-Integer(2), Integer(0), -Integer(10), Integer(2)),
...             (-Integer(1), Integer(0), Integer(15), -Integer(2)),
...             (Integer(0), Integer(1), Integer(0), -Integer(1))])
>>> a.eigenvalues()                                                       # needs sage.rings.number_field
[-0.9386318578049146?,
 15.50655435353258?,
 0.2160387521361705? - 4.713151979747493?*I,
 0.2160387521361705? + 4.713151979747493?*I]

A symmetric matrix a + a.transpose() should have real eigenvalues

sage: b = a + a.transpose()
sage: ev = b.eigenvalues(); ev                                              # needs sage.rings.number_field
[-8.35066086057957?, -1.107247901349379?,
 5.718651326708515?, 33.73925743522043?]
>>> from sage.all import *
>>> b = a + a.transpose()
>>> ev = b.eigenvalues(); ev                                              # needs sage.rings.number_field
[-8.35066086057957?, -1.107247901349379?,
 5.718651326708515?, 33.73925743522043?]

The eigenvalues are elements of QQbar, so they really represent exact roots of polynomials, not just approximations.

sage: e = ev[0]; e                                                          # needs sage.rings.number_field
-8.35066086057957?
sage: p = e.minpoly(); p                                                    # needs sage.rings.number_field
x^4 - 30*x^3 - 171*x^2 + 1460*x + 1784
sage: p(e) == 0                                                             # needs sage.rings.number_field
True
>>> from sage.all import *
>>> e = ev[Integer(0)]; e                                                          # needs sage.rings.number_field
-8.35066086057957?
>>> p = e.minpoly(); p                                                    # needs sage.rings.number_field
x^4 - 30*x^3 - 171*x^2 + 1460*x + 1784
>>> p(e) == Integer(0)                                                             # needs sage.rings.number_field
True

To perform computations on the eigenvalue as an element of a number field, you can always convert back to a number field element.

sage: e.as_number_field_element()                                           # needs sage.rings.number_field
(Number Field in a
  with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264,
 a + 7,
 Ring morphism:
   From: Number Field in a with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264
   To:   Algebraic Real Field
   Defn: a |--> -15.35066086057957?)
>>> from sage.all import *
>>> e.as_number_field_element()                                           # needs sage.rings.number_field
(Number Field in a
  with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264,
 a + 7,
 Ring morphism:
   From: Number Field in a with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264
   To:   Algebraic Real Field
   Defn: a |--> -15.35066086057957?)

Notice the effect of the extend option.

sage: M = matrix(QQ, [[0,-1,0], [1,0,0], [0,0,2]])
sage: M.eigenvalues()                                                       # needs sage.rings.number_field
[2, -1*I, 1*I]
sage: M.eigenvalues(extend=False)                                           # needs sage.libs.pari
[2]
>>> from sage.all import *
>>> M = matrix(QQ, [[Integer(0),-Integer(1),Integer(0)], [Integer(1),Integer(0),Integer(0)], [Integer(0),Integer(0),Integer(2)]])
>>> M.eigenvalues()                                                       # needs sage.rings.number_field
[2, -1*I, 1*I]
>>> M.eigenvalues(extend=False)                                           # needs sage.libs.pari
[2]

The method also works for matrices over finite fields:

sage: M = matrix(GF(3), [[0,1,1], [1,2,0], [2,0,1]])
sage: ev = sorted(M.eigenvalues()); ev                                      # needs sage.rings.finite_rings
[2*z3, 2*z3 + 1, 2*z3 + 2]
>>> from sage.all import *
>>> M = matrix(GF(Integer(3)), [[Integer(0),Integer(1),Integer(1)], [Integer(1),Integer(2),Integer(0)], [Integer(2),Integer(0),Integer(1)]])
>>> ev = sorted(M.eigenvalues()); ev                                      # needs sage.rings.finite_rings
[2*z3, 2*z3 + 1, 2*z3 + 2]

Similarly as in the case of QQbar, the eigenvalues belong to some algebraic closure but they can be converted to elements of a finite field:

sage: e = ev[0]                                                             # needs sage.rings.finite_rings
sage: e.parent()                                                            # needs sage.rings.finite_rings
Algebraic closure of Finite Field of size 3
sage: e.as_finite_field_element()                                           # needs sage.rings.finite_rings
(Finite Field in z3 of size 3^3,
 2*z3,
 Ring morphism:
   From: Finite Field in z3 of size 3^3
   To:   Algebraic closure of Finite Field of size 3
   Defn: z3 |--> z3)
>>> from sage.all import *
>>> e = ev[Integer(0)]                                                             # needs sage.rings.finite_rings
>>> e.parent()                                                            # needs sage.rings.finite_rings
Algebraic closure of Finite Field of size 3
>>> e.as_finite_field_element()                                           # needs sage.rings.finite_rings
(Finite Field in z3 of size 3^3,
 2*z3,
 Ring morphism:
   From: Finite Field in z3 of size 3^3
   To:   Algebraic closure of Finite Field of size 3
   Defn: z3 |--> z3)
eigenvectors_left(other=None, extend=True)[source]

Compute the left eigenvectors of a matrix.

INPUT:

  • other – a square matrix \(B\) (default: None) in a generalized eigenvalue problem; if None, an ordinary eigenvalue problem is solved (currently supported only if the base ring of self is RDF or CDF)

  • extend – boolean (default: True)

OUTPUT:

For each distinct eigenvalue, returns a list of the form (e,V,n) where e is the eigenvalue, V is a list of eigenvectors forming a basis for the corresponding left eigenspace, and n is the algebraic multiplicity of the eigenvalue.

If the option extend is set to False, then only the eigenvalues that live in the base ring are considered.

EXAMPLES:

We compute the left eigenvectors of a \(3\times 3\) rational matrix.

sage: # needs sage.rings.number_field
sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: es = A.eigenvectors_left(); es
[(0, [ (1, -2, 1) ], 1),
 (-1.348469228349535?, [(1, 0.3101020514433644?, -0.3797958971132713?)], 1),
 (13.34846922834954?, [(1, 1.289897948556636?, 1.579795897113272?)], 1)]
sage: eval, [evec], mult = es[0]
sage: delta = eval*evec - evec*A
sage: abs(abs(delta)) < 1e-10
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> es = A.eigenvectors_left(); es
[(0, [ (1, -2, 1) ], 1),
 (-1.348469228349535?, [(1, 0.3101020514433644?, -0.3797958971132713?)], 1),
 (13.34846922834954?, [(1, 1.289897948556636?, 1.579795897113272?)], 1)]
>>> eval, [evec], mult = es[Integer(0)]
>>> delta = eval*evec - evec*A
>>> abs(abs(delta)) < RealNumber('1e-10')
True

Notice the difference between considering ring extensions or not.

sage: M = matrix(QQ, [[0,-1,0], [1,0,0], [0,0,2]])
sage: M.eigenvectors_left()                                                 # needs sage.rings.number_field
[(2,    [ (0, 0, 1) ],  1),
 (-1*I, [(1, -1*I, 0)], 1),
 (1*I,  [(1, 1*I, 0)],  1)]
sage: M.eigenvectors_left(extend=False)                                     # needs sage.rings.number_field
[(2,    [ (0, 0, 1) ],  1)]
>>> from sage.all import *
>>> M = matrix(QQ, [[Integer(0),-Integer(1),Integer(0)], [Integer(1),Integer(0),Integer(0)], [Integer(0),Integer(0),Integer(2)]])
>>> M.eigenvectors_left()                                                 # needs sage.rings.number_field
[(2,    [ (0, 0, 1) ],  1),
 (-1*I, [(1, -1*I, 0)], 1),
 (1*I,  [(1, 1*I, 0)],  1)]
>>> M.eigenvectors_left(extend=False)                                     # needs sage.rings.number_field
[(2,    [ (0, 0, 1) ],  1)]
eigenvectors_right(other=None, extend=True)[source]

Compute the right eigenvectors of a matrix.

INPUT:

  • other – a square matrix \(B\) (default: None) in a generalized eigenvalue problem; if None, an ordinary eigenvalue problem is solved (currently supported only if the base ring of self is RDF or CDF)

  • extend – boolean (default: True)

OUTPUT:

For each distinct eigenvalue, returns a list of the form (e,V,n) where e is the eigenvalue, V is a list of eigenvectors forming a basis for the corresponding right eigenspace, and n is the algebraic multiplicity of the eigenvalue. If extend = True (the default), this will return eigenspaces over the algebraic closure of the base field where this is implemented; otherwise it will restrict to eigenvalues in the base field.

EXAMPLES:

We compute the right eigenvectors of a \(3\times 3\) rational matrix.

sage: # needs sage.rings.number_field
sage: A = matrix(QQ, 3, 3, range(9)); A
[0 1 2]
[3 4 5]
[6 7 8]
sage: es = A.eigenvectors_right(); es
[(0, [ (1, -2, 1) ], 1),
 (-1.348469228349535?, [(1, 0.1303061543300932?, -0.7393876913398137?)], 1),
 (13.34846922834954?, [(1, 3.069693845669907?, 5.139387691339814?)], 1)]
sage: A.eigenvectors_right(extend=False)
[(0, [ (1, -2, 1) ], 1)]
sage: eval, [evec], mult = es[0]
sage: delta = eval*evec - A*evec
sage: abs(abs(delta)) < 1e-10
True
>>> from sage.all import *
>>> # needs sage.rings.number_field
>>> A = matrix(QQ, Integer(3), Integer(3), range(Integer(9))); A
[0 1 2]
[3 4 5]
[6 7 8]
>>> es = A.eigenvectors_right(); es
[(0, [ (1, -2, 1) ], 1),
 (-1.348469228349535?, [(1, 0.1303061543300932?, -0.7393876913398137?)], 1),
 (13.34846922834954?, [(1, 3.069693845669907?, 5.139387691339814?)], 1)]
>>> A.eigenvectors_right(extend=False)
[(0, [ (1, -2, 1) ], 1)]
>>> eval, [evec], mult = es[Integer(0)]
>>> delta = eval*evec - A*evec
>>> abs(abs(delta)) < RealNumber('1e-10')
True
elementary_divisors(algorithm=None)[source]

If self is a matrix over a principal ideal domain \(R\), return elements \(d_i\) for \(1 \le i \le k = \min(r,s)\) where \(r\) and \(s\) are the number of rows and columns of self, such that the cokernel of self is isomorphic to

\[R/(d_1) \oplus R/(d_2) \oplus R/(d_k)\]

with \(d_i \mid d_{i+1}\) for all \(i\). These are the diagonal entries of the Smith form of self (see smith_form()).

INPUT:

  • algorithm – ignored

EXAMPLES:

sage: x = polygen(ZZ, 'x')
sage: OE.<w> = EquationOrder(x^2 - x + 2)                                   # needs sage.rings.number_field
sage: m = Matrix([[1, w], [w, 7]])                                          # needs sage.rings.number_field
sage: m.elementary_divisors()                                               # needs sage.rings.number_field
[1, -w + 9]
>>> from sage.all import *
>>> x = polygen(ZZ, 'x')
>>> OE = EquationOrder(x**Integer(2) - x + Integer(2), names=('w',)); (w,) = OE._first_ngens(1)# needs sage.rings.number_field
>>> m = Matrix([[Integer(1), w], [w, Integer(7)]])                                          # needs sage.rings.number_field
>>> m.elementary_divisors()                                               # needs sage.rings.number_field
[1, -w + 9]

See also

smith_form()

elementwise_product(right)[source]

Return the elementwise product of two matrices of the same size (also known as the Hadamard product).

INPUT:

  • right – the right operand of the product. A matrix of the same size as self such that multiplication of elements of the base rings of self and right is defined, once Sage’s coercion model is applied. If the matrices have different sizes, or if multiplication of individual entries cannot be achieved, a TypeError will result.

OUTPUT:

A matrix of the same size as self and right. The entry in location \((i,j)\) of the output is the product of the two entries in location \((i,j)\) of self and right (in that order).

The parent of the result is determined by Sage’s coercion model. If the base rings are identical, then the result is dense or sparse according to this property for the left operand. If the base rings must be adjusted for one, or both, matrices then the result will be sparse only if both operands are sparse. No subdivisions are present in the result.

If the type of the result is not to your liking, or the ring could be “tighter,” adjust the operands with change_ring(). Adjust sparse versus dense inputs with the methods sparse_matrix() and dense_matrix().

EXAMPLES:

sage: A = matrix(ZZ, 2, 3, range(6))
sage: B = matrix(QQ, 2, 3, [5, 1/3, 2/7, 11/2, -3/2, 8])
sage: C = A.elementwise_product(B)
sage: C
[   0  1/3  4/7]
[33/2   -6   40]
sage: C.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field
>>> from sage.all import *
>>> A = matrix(ZZ, Integer(2), Integer(3), range(Integer(6)))
>>> B = matrix(QQ, Integer(2), Integer(3), [Integer(5), Integer(1)/Integer(3), Integer(2)/Integer(7), Integer(11)/Integer(2), -Integer(3)/Integer(2), Integer(8)])
>>> C = A.elementwise_product(B)
>>> C
[   0  1/3  4/7]
[33/2   -6   40]
>>> C.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field

Notice the base ring of the results in the next two examples.

sage: x = polygen(ZZ, 'x')
sage: D = matrix(ZZ['x'],2,[1+x^2,2,3,4-x])
sage: E = matrix(QQ,2,[1,2,3,4])
sage: F = D.elementwise_product(E)
sage: F
[  x^2 + 1         4]
[        9 -4*x + 16]
sage: F.parent()
Full MatrixSpace of 2 by 2 dense matrices
 over Univariate Polynomial Ring in x over Rational Field
>>> from sage.all import *
>>> x = polygen(ZZ, 'x')
>>> D = matrix(ZZ['x'],Integer(2),[Integer(1)+x**Integer(2),Integer(2),Integer(3),Integer(4)-x])
>>> E = matrix(QQ,Integer(2),[Integer(1),Integer(2),Integer(3),Integer(4)])
>>> F = D.elementwise_product(E)
>>> F
[  x^2 + 1         4]
[        9 -4*x + 16]
>>> F.parent()
Full MatrixSpace of 2 by 2 dense matrices
 over Univariate Polynomial Ring in x over Rational Field

sage: G = matrix(GF(3), 2, [0, 1, 2, 2])
sage: H = matrix(ZZ, 2, [1, 2, 3, 4])
sage: J = G.elementwise_product(H)
sage: J
[0 2]
[0 2]
sage: J.parent()
Full MatrixSpace of 2 by 2 dense matrices
 over Finite Field of size 3
>>> from sage.all import *
>>> G = matrix(GF(Integer(3)), Integer(2), [Integer(0), Integer(1), Integer(2), Integer(2)])
>>> H = matrix(ZZ, Integer(2), [Integer(1), Integer(2), Integer(3), Integer(4)])
>>> J = G.elementwise_product(H)
>>> J
[0 2]
[0 2]
>>> J.parent()
Full MatrixSpace of 2 by 2 dense matrices
 over Finite Field of size 3

Non-commutative rings behave as expected. These are the usual quaternions.

sage: # needs sage.combinat
sage: R.<i,j,k> = QuaternionAlgebra(-1, -1)
sage: A = matrix(R, 2, [1,i,j,k])
sage: B = matrix(R, 2, [i,i,i,i])
sage: A.elementwise_product(B)
[ i -1]
[-k  j]
sage: B.elementwise_product(A)
[ i -1]
[ k -j]
>>> from sage.all import *
>>> # needs sage.combinat
>>> R = QuaternionAlgebra(-Integer(1), -Integer(1), names=('i', 'j', 'k',)); (i, j, k,) = R._first_ngens(3)
>>> A = matrix(R, Integer(2), [Integer(1),i,j,k])
>>> B = matrix(R, Integer(2), [i,i,i,i])
>>> A.elementwise_product(B)
[ i -1]
[-k  j]
>>> B.elementwise_product(A)
[ i -1]
[ k -j]

Input that is not a matrix will raise an error.

sage: A = random_matrix(ZZ, 5, 10, x=20)
sage: A.elementwise_product(vector(ZZ, [1,2,3,4]))
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
'Full MatrixSpace of 5 by 10 dense matrices over Integer Ring' and
'Ambient free module of rank 4 over the principal ideal domain Integer Ring'

sage: A = matrix(2, 2, range(4))
sage: A.elementwise_product(polygen(parent(A)))
Traceback (most recent call last):
...
TypeError: elementwise_product() argument should be a matrix
or coercible to a matrix
>>> from sage.all import *
>>> A = random_matrix(ZZ, Integer(5), Integer(10), x=Integer(20))
>>> A.elementwise_product(vector(ZZ, [Integer(1),Integer(2),Integer(3),Integer(4)]))
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
'Full MatrixSpace of 5 by 10 dense matrices over Integer Ring' and
'Ambient free module of rank 4 over the principal ideal domain Integer Ring'

>>> A = matrix(Integer(2), Integer(2), range(Integer(4)))
>>> A.elementwise_product(polygen(parent(A)))
Traceback (most recent call last):
...
TypeError: elementwise_product() argument should be a matrix
or coercible to a matrix

Matrices of different sizes for operands will raise an error.

sage: A = random_matrix(ZZ, 5, 10, x=20)
sage: B = random_matrix(ZZ, 10, 5, x=40)
sage: A.elementwise_product(B)
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
'Full MatrixSpace of 5 by 10 dense matrices over Integer Ring' and
'Full MatrixSpace of 10 by 5 dense matrices over Integer Ring'
>>> from sage.all import *
>>> A = random_matrix(ZZ, Integer(5), Integer(10), x=Integer(20))
>>> B = random_matrix(ZZ, Integer(10), Integer(5), x=Integer(40))
>>> A.elementwise_product(B)
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
'Full MatrixSpace of 5 by 10 dense matrices over Integer Ring' and
'Full MatrixSpace of 10 by 5 dense matrices over Integer Ring'

Some pairs of rings do not have a common parent where multiplication makes sense. This will raise an error.

sage: A = matrix(QQ, 3, 2, range(6))
sage: B = matrix(GF(3), 3, [2]*6)
sage: A.elementwise_product(B)
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
 'Full MatrixSpace of 3 by 2 dense matrices over Rational Field' and
 'Full MatrixSpace of 3 by 2 dense matrices over Finite Field of size 3'
>>> from sage.all import *
>>> A = matrix(QQ, Integer(3), Integer(2), range(Integer(6)))
>>> B = matrix(GF(Integer(3)), Integer(3), [Integer(2)]*Integer(6))
>>> A.elementwise_product(B)
Traceback (most recent call last):
...
TypeError: no common canonical parent for objects with parents:
 'Full MatrixSpace of 3 by 2 dense matrices over Rational Field' and
 'Full MatrixSpace of 3 by 2 dense matrices over Finite Field of size 3'

We illustrate various combinations of sparse and dense matrices. The usual coercion rules apply:

sage: A = matrix(ZZ, 5, 6, range(30), sparse=False)
sage: B = matrix(ZZ, 5, 6, range(30), sparse=True)
sage: C = matrix(QQ, 5, 6, range(30), sparse=True)
sage: A.elementwise_product(C).is_sparse()
True
sage: B.elementwise_product(C).is_sparse()
True
sage: A.elementwise_product(B).is_dense()
True
sage: B.elementwise_product(A).is_dense()
True
>>> from sage.all import *
>>> A = matrix(ZZ, Integer(5), Integer(6), range(Integer(30)), sparse=False)
>>> B = matrix(ZZ, Integer(5), Integer(6), range(Integer(30)), sparse=True)
>>> C = matrix(QQ, Integer(5), Integer(6), range(Integer(30)), sparse=True)
>>> A.elementwise_product(C).is_sparse()
True
>>> B.elementwise_product(C).is_sparse()
True
>>> A.elementwise_product(B).is_dense()
True
>>> B.elementwise_product(A).is_dense()
True
exp()[source]

Calculate the exponential of this matrix X, which is the matrix.

\[e^X = \sum_{k=0}^{\infty} \frac{X^k}{k!}.\]

This function depends on maxima’s matrix exponentiation function, which does not deal well with floating point numbers. If the matrix has floating point numbers, they will be rounded automatically to rational numbers during the computation. If you want approximations to the exponential that are calculated numerically, you may get better results by first converting your matrix to RDF or CDF, as shown in the last example.

EXAMPLES:

sage: # needs sage.symbolic
sage: a = matrix([[1,2], [3,4]])
sage: a.exp()
[-1/22*((sqrt(33) - 11)*e^sqrt(33) - sqrt(33) - 11)*e^(-1/2*sqrt(33) + 5/2)              2/33*(sqrt(33)*e^sqrt(33) - sqrt(33))*e^(-1/2*sqrt(33) + 5/2)]
[             1/11*(sqrt(33)*e^sqrt(33) - sqrt(33))*e^(-1/2*sqrt(33) + 5/2)  1/22*((sqrt(33) + 11)*e^sqrt(33) - sqrt(33) + 11)*e^(-1/2*sqrt(33) + 5/2)]

sage: type(a.exp())                                                         # needs sage.symbolic
<class 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>

sage: a = matrix([[1/2,2/3], [3/4,4/5]])
sage: a.exp()                                                               # needs sage.symbolic
[-1/418*((3*sqrt(209) - 209)*e^(1/10*sqrt(209)) - 3*sqrt(209) - 209)*e^(-1/20*sqrt(209) + 13/20)                   20/627*(sqrt(209)*e^(1/10*sqrt(209)) - sqrt(209))*e^(-1/20*sqrt(209) + 13/20)]
[                  15/418*(sqrt(209)*e^(1/10*sqrt(209)) - sqrt(209))*e^(-1/20*sqrt(209) + 13/20)  1/418*((3*sqrt(209) + 209)*e^(1/10*sqrt(209)) - 3*sqrt(209) + 209)*e^(-1/20*sqrt(209) + 13/20)]

sage: a = matrix(RR, [[1,pi.n()], [1e2,1e-2]])                              # needs sage.symbolic
sage: a.exp()                                                               # needs sage.symbolic
[ 1/11882424341266*((11*sqrt(227345670387496707609) + 5941212170633)*e^(3/1275529100*sqrt(227345670387496707609)) - 11*sqrt(227345670387496707609) + 5941212170633)*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)                            445243650/75781890129165569203*(sqrt(227345670387496707609)*e^(3/1275529100*sqrt(227345670387496707609)) - sqrt(227345670387496707609))*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)]
[                                     10000/53470909535697*(sqrt(227345670387496707609)*e^(3/1275529100*sqrt(227345670387496707609)) - sqrt(227345670387496707609))*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200) -1/11882424341266*((11*sqrt(227345670387496707609) - 5941212170633)*e^(3/1275529100*sqrt(227345670387496707609)) - 11*sqrt(227345670387496707609) - 5941212170633)*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)]
sage: a.change_ring(RDF).exp()  # rel tol 6e-14                             # needs sage.symbolic
[42748127.31532951 7368259.244159399]
[234538976.1381042 40426191.45156228]
>>> from sage.all import *
>>> # needs sage.symbolic
>>> a = matrix([[Integer(1),Integer(2)], [Integer(3),Integer(4)]])
>>> a.exp()
[-1/22*((sqrt(33) - 11)*e^sqrt(33) - sqrt(33) - 11)*e^(-1/2*sqrt(33) + 5/2)              2/33*(sqrt(33)*e^sqrt(33) - sqrt(33))*e^(-1/2*sqrt(33) + 5/2)]
[             1/11*(sqrt(33)*e^sqrt(33) - sqrt(33))*e^(-1/2*sqrt(33) + 5/2)  1/22*((sqrt(33) + 11)*e^sqrt(33) - sqrt(33) + 11)*e^(-1/2*sqrt(33) + 5/2)]

>>> type(a.exp())                                                         # needs sage.symbolic
<class 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>

>>> a = matrix([[Integer(1)/Integer(2),Integer(2)/Integer(3)], [Integer(3)/Integer(4),Integer(4)/Integer(5)]])
>>> a.exp()                                                               # needs sage.symbolic
[-1/418*((3*sqrt(209) - 209)*e^(1/10*sqrt(209)) - 3*sqrt(209) - 209)*e^(-1/20*sqrt(209) + 13/20)                   20/627*(sqrt(209)*e^(1/10*sqrt(209)) - sqrt(209))*e^(-1/20*sqrt(209) + 13/20)]
[                  15/418*(sqrt(209)*e^(1/10*sqrt(209)) - sqrt(209))*e^(-1/20*sqrt(209) + 13/20)  1/418*((3*sqrt(209) + 209)*e^(1/10*sqrt(209)) - 3*sqrt(209) + 209)*e^(-1/20*sqrt(209) + 13/20)]

>>> a = matrix(RR, [[Integer(1),pi.n()], [RealNumber('1e2'),RealNumber('1e-2')]])                              # needs sage.symbolic
>>> a.exp()                                                               # needs sage.symbolic
[ 1/11882424341266*((11*sqrt(227345670387496707609) + 5941212170633)*e^(3/1275529100*sqrt(227345670387496707609)) - 11*sqrt(227345670387496707609) + 5941212170633)*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)                            445243650/75781890129165569203*(sqrt(227345670387496707609)*e^(3/1275529100*sqrt(227345670387496707609)) - sqrt(227345670387496707609))*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)]
[                                     10000/53470909535697*(sqrt(227345670387496707609)*e^(3/1275529100*sqrt(227345670387496707609)) - sqrt(227345670387496707609))*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200) -1/11882424341266*((11*sqrt(227345670387496707609) - 5941212170633)*e^(3/1275529100*sqrt(227345670387496707609)) - 11*sqrt(227345670387496707609) - 5941212170633)*e^(-3/2551058200*sqrt(227345670387496707609) + 101/200)]
>>> a.change_ring(RDF).exp()  # rel tol 6e-14                             # needs sage.symbolic
[42748127.31532951 7368259.244159399]
[234538976.1381042 40426191.45156228]
extended_echelon_form(subdivide=False, **kwds)[source]

Return the echelon form of self augmented with an identity matrix.

INPUT:

  • subdivide – boolean (default: False); whether to subdivide the returned matrix. See the description of the (output) below for details.

  • kwds – additional keywords that can be passed to the method that computes the echelon form

OUTPUT:

If \(A\) is an \(m\times n\) matrix, add the \(m\) columns of an \(m\times m\) identity matrix to the right of self. Then row-reduce this \(m\times(n+m)\) matrix. This matrix is returned as an immutable matrix.

If subdivide is True then the returned matrix has a single division among the columns and a single division among the rows. The column subdivision has \(n\) columns to the left and \(m\) columns to the right. The row division separates the nonzero rows from the zero rows, when restricted to the first \(n\) columns.

For a nonsingular matrix the final \(m\) columns of the extended echelon form are the inverse of self. For a matrix of any size, the final \(m\) columns provide a matrix that transforms self to echelon form when it multiplies self from the left. When the base ring is a field, the uniqueness of reduced row-echelon form implies that this transformation matrix can be taken as the coefficients giving a canonical set of linear combinations of the rows of self that yield reduced row-echelon form.

When subdivided as described above, and again over a field, the parts of the subdivision in the upper-left corner and lower-right corner satisfy several interesting relationships with the row space, column space, left kernel and right kernel of self. See the examples below.

Note

This method returns an echelon form. If the base ring is not a field, no attempt is made to move to the fraction field. See an example below where the base ring is changed manually.

EXAMPLES:

The four relationships at the end of this example hold in general.

sage: A = matrix(QQ, [[2, -1, 7, -1, 0, -3],
....:                 [-1, 1, -5, 3, 4, 4],
....:                 [2, -1, 7, 0, 2, -2],
....:                 [2, 0, 4, 3, 6, 1],
....:                 [2, -1, 7, 0, 2, -2]])
sage: E = A.extended_echelon_form(subdivide=True); E
[   1    0    2    0    0   -1|   0   -1    0    1   -1]
[   0    1   -3    0   -2    0|   0   -2    0    2   -3]
[   0    0    0    1    2    1|   0  2/3    0 -1/3  2/3]
[-----------------------------+------------------------]
[   0    0    0    0    0    0|   1  2/3    0 -1/3 -1/3]
[   0    0    0    0    0    0|   0    0    1    0   -1]
sage: J = E.matrix_from_columns(range(6,11)); J
[   0   -1    0    1   -1]
[   0   -2    0    2   -3]
[   0  2/3    0 -1/3  2/3]
[   1  2/3    0 -1/3 -1/3]
[   0    0    1    0   -1]
sage: J*A == A.rref()
True
sage: C = E.subdivision(0,0); C
[ 1  0  2  0  0 -1]
[ 0  1 -3  0 -2  0]
[ 0  0  0  1  2  1]
sage: L = E.subdivision(1,1); L
[   1  2/3    0 -1/3 -1/3]
[   0    0    1    0   -1]
sage: A.right_kernel() == C.right_kernel()
True
sage: A.row_space() == C.row_space()
True
sage: A.column_space() == L.right_kernel()
True
sage: A.left_kernel() == L.row_space()
True
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(2), -Integer(1), Integer(7), -Integer(1), Integer(0), -Integer(3)],
...                 [-Integer(1), Integer(1), -Integer(5), Integer(3), Integer(4), Integer(4)],
...                 [Integer(2), -Integer(1), Integer(7), Integer(0), Integer(2), -Integer(2)],
...                 [Integer(2), Integer(0), Integer(4), Integer(3), Integer(6), Integer(1)],
...                 [Integer(2), -Integer(1), Integer(7), Integer(0), Integer(2), -Integer(2)]])
>>> E = A.extended_echelon_form(subdivide=True); E
[   1    0    2    0    0   -1|   0   -1    0    1   -1]
[   0    1   -3    0   -2    0|   0   -2    0    2   -3]
[   0    0    0    1    2    1|   0  2/3    0 -1/3  2/3]
[-----------------------------+------------------------]
[   0    0    0    0    0    0|   1  2/3    0 -1/3 -1/3]
[   0    0    0    0    0    0|   0    0    1    0   -1]
>>> J = E.matrix_from_columns(range(Integer(6),Integer(11))); J
[   0   -1    0    1   -1]
[   0   -2    0    2   -3]
[   0  2/3    0 -1/3  2/3]
[   1  2/3    0 -1/3 -1/3]
[   0    0    1    0   -1]
>>> J*A == A.rref()
True
>>> C = E.subdivision(Integer(0),Integer(0)); C
[ 1  0  2  0  0 -1]
[ 0  1 -3  0 -2  0]
[ 0  0  0  1  2  1]
>>> L = E.subdivision(Integer(1),Integer(1)); L
[   1  2/3    0 -1/3 -1/3]
[   0    0    1    0   -1]
>>> A.right_kernel() == C.right_kernel()
True
>>> A.row_space() == C.row_space()
True
>>> A.column_space() == L.right_kernel()
True
>>> A.left_kernel() == L.row_space()
True

For a nonsingular matrix, the right half of the extended echelon form is the inverse matrix.

sage: B = matrix(QQ, [[1,3,4], [1,4,4], [0,-2,-1]])
sage: E = B.extended_echelon_form()
sage: J = E.matrix_from_columns(range(3,6)); J
[-4  5  4]
[-1  1  0]
[ 2 -2 -1]
sage: J == B.inverse()
True
>>> from sage.all import *
>>> B = matrix(QQ, [[Integer(1),Integer(3),Integer(4)], [Integer(1),Integer(4),Integer(4)], [Integer(0),-Integer(2),-Integer(1)]])
>>> E = B.extended_echelon_form()
>>> J = E.matrix_from_columns(range(Integer(3),Integer(6))); J
[-4  5  4]
[-1  1  0]
[ 2 -2 -1]
>>> J == B.inverse()
True

The result is in echelon form, so if the base ring is not a field, the leading entry of each row may not be 1. But you can easily change to the fraction field if necessary.

sage: A = matrix(ZZ, [[16, 20, 4, 5, -4, 13, 5],
....:                 [10, 13, 3, -3, 7, 11, 6],
....:                 [-12, -15, -3, -3, 2, -10, -4],
....:                 [10, 13, 3, 3, -1, 9, 4],
....:                 [4, 5, 1, 8, -10, 1, -1]])
sage: E = A.extended_echelon_form(subdivide=True); E
[ 2  0 -2  2 -9 -3 -4| 0  4 -3 -9  4]
[ 0  1  1  2  0  1  1| 0  1  2  1  1]
[ 0  0  0  3 -4 -1 -1| 0  3  1 -3  3]
[--------------------+--------------]
[ 0  0  0  0  0  0  0| 1  6  3 -6  5]
[ 0  0  0  0  0  0  0| 0  7  2 -7  6]
sage: J = E.matrix_from_columns(range(7,12)); J
[ 0  4 -3 -9  4]
[ 0  1  2  1  1]
[ 0  3  1 -3  3]
[ 1  6  3 -6  5]
[ 0  7  2 -7  6]
sage: J*A == A.echelon_form()
True
sage: B = A.change_ring(QQ)
sage: B.extended_echelon_form(subdivide=True)
[     1      0     -1      0  -19/6   -7/6   -5/3|     0      0 -89/42   -5/2    1/7]
[     0      1      1      0    8/3    5/3    5/3|     0      0  34/21      2   -1/7]
[     0      0      0      1   -4/3   -1/3   -1/3|     0      0   1/21      0    1/7]
[------------------------------------------------+----------------------------------]
[     0      0      0      0      0      0      0|     1      0    9/7      0   -1/7]
[     0      0      0      0      0      0      0|     0      1    2/7     -1    6/7]
>>> from sage.all import *
>>> A = matrix(ZZ, [[Integer(16), Integer(20), Integer(4), Integer(5), -Integer(4), Integer(13), Integer(5)],
...                 [Integer(10), Integer(13), Integer(3), -Integer(3), Integer(7), Integer(11), Integer(6)],
...                 [-Integer(12), -Integer(15), -Integer(3), -Integer(3), Integer(2), -Integer(10), -Integer(4)],
...                 [Integer(10), Integer(13