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_vectormethodsRobert Bradshaw (2007-06-14): added
subdividemethodJaap Spies (2007-11-14): implemented
_binomial,_chooseauxiliary functionsWilliam Stein (2007-11-18): added
_gram_schmidt_noscalemethodDavid Loeffler (2008-12-05): added
smith_formmethodDavid Loeffler (2009-06-01): added
_echelon_form_PIDmethodSebastian Pancratz (2009-06-25): implemented
adjointandcharpolymethods; fixedadjointreflecting the change that_adjointis now implemented inMatrix; used the division-free algorithm forcharpolyRob Beezer (2009-07-13): added
elementwise_productmethodMiguel Marco (2010-06-19): modified eigenvalues and eigenvectors functions to allow the option
extend=FalseThierry Monteil (2010-10-05): bugfix for Issue #10063, so that the determinant is computed even for rings for which the
is_fieldmethod is not implemented.Rob Beezer (2010-12-13): added
conjugate_transposemethodRob Beezer (2011-02-05): refactored all of the matrix kernel routines; added
extended_echelon_form,right_kernel_matrix,QR,_gram_schmidt_noscale,is_similarmethodsMoritz Minzlaff (2011-03-17): corrected
_echelon_form_PIDmethod for matrices of one row, fixed in Issue #9053Rob Beezer (2011-06-09): added
is_normal,is_diagonalizable,LU,cyclic_subspace,zigzag_form,rational_formmethodsRob Beezer (2012-05-27): added
indefinite_factorization,is_positive_definite,choleskymethodsDarij Grinberg (2013-10-01): added first (slow) pfaffian implementation
Mario Pernici (2014-07-01): modified
rook_vectormethodRob Beezer (2015-05-25): modified
is_similarmethodSamuel Lelièvre (2020-09-18): improved method
LLL_grambased on a patch by William Stein posted at Issue #5178, moving the method from its initial location insage.matrix.integer_matrix_denseMichael Jung (2020-10-02): added Bär-Faddeev-LeVerrier algorithm for the Pfaffian
Moritz Firsching(2020-10-05): added
quantum_determinantDima Pasechnik (2022-11-08): fixed
echelonizefor inexact matrices
- class sage.matrix.matrix2.Matrix[source]¶
Bases:
MatrixBase 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
selfis 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 formflag– an optional flag passed toqflllgramAccording to pari:qflllgram’s documentation the options are:0– (default), assume thatselfhas either exact (integral or rational) or real floating point entries. The matrix is rescaled, converted to integers and the behavior is then as inflag=1.1– assume that G is integral. Computations involving Gram-Schmidt vectors are approximate, with precision varying as needed.
OUTPUT:
A dense matrix
Uover the integers with determinant 1 such thatU.T * M * Uis 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,
RDFor 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
RDFandCDF). 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
selfas a unitary matrix and an upper-triangular matrix.INPUT:
full– (default:True) ifTruethen the returned matrices have dimensions as described below. IfFalsetheRmatrix has no zero rows and the columns ofQare a basis for the column space ofself.
OUTPUT:
If
selfis an \(m\times n\) matrix andfull=Truethen 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=Falsethen \(Q\) has \(m\) rows and the columns form an orthonormal basis for the column space ofself. 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=Trueare 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
fullkeyword.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 deprecatedadjoint()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._adjugateis 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 deprecatedadjoint()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._adjugateis 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); ``Trueto make the output a sparse matrixphi– arbitrary Python function or callable objectR– (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
selfis 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, sophiis callable andphi.domain()andphi.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
Brepresenting the matrix uniquely.Vertices are labeled 1 to
nrowson the left andnrows+ 1 tonrows+ncolson 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
bistochastic_as_sum_of_permutations– for more information on this method.
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
selfunder permutations of rows and columns as a list of pairs ofPermutationGroupElementobjects.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
Mto 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=Trueat 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. AValueErroris 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].
See also
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
Sin 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
ValueErrorresults.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 classicalblock_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
copyfunction 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
selfafter 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
selfand 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 fieldv– a vector with a degree equal to the size of the matrix and entries compatible with the entries of the matrixvar– (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
varkeyword, 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 matrixalgorithm– string (default:'spin');'spin': involves iterating the action ofselfon 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 isTrueexactly when the charpoly ofselfon V is irreducible.(optional) list – list of pairs (W,t), where W is a vector space and t is a boolean, and t is
Trueexactly when the charpoly of the transpose ofselfon 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
selfon M leaves M invariant. Return the decomposition of M as a list of pairs (W, is_irred) where is_irred isTrueif the charpoly ofselfacting on the factor W is irreducible.Additional inputs besides M are passed onto the decomposition command.
INPUT:
M– a subspace of the free moduleselfacts oncheck_restrict– boolean (default:True); call restrict with or without checkkwds– keywords that will be forwarded todecomposition()
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. (TypeA.determinant?to see what is done for a specific matrixA.)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
selfalong with the transformation matrix.INPUT:
base_field– if given,selfis regarded as a matrix over it
OUTPUT: a diagonal matrix \(D\) and an invertible matrix \(P\) such that \(P^{-1}AP=D\), if
selfis 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 selectedtransformation– 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 thatselfis not changed by this command. Useechelonize()to changeselfin place.If the optional parameter
transformation=Trueis specified, the output consists of a pair \((E,T)\) of matrices where \(E\) is the echelon form ofselfand \(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
transformationoption 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
selfinto a matrix in echelon form over the same base ring asself.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 selectedtransformation– 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
selfis put into echelon form. Nothing is returned unless the keyword optiontransformation=Trueis 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; ifNone, an ordinary eigenvalue problem is solved
OUTPUT:
If
selfis 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
RDFandCDF.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; ifNone, an ordinary eigenvalue problem is solved
OUTPUT:
If
selfis 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
RDFandCDF.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()andleft_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 classesformat– 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 factorNone– 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. Ifvar='a', then the root fields will be in terms ofa0, 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 theformatkeyword.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 withformat='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()andright_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 classesformat– (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 factorNone– 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
formatkeyword.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
sas 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, thenQQbarelements are returned that represent each separate root.If the option
extendis set toFalse, 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 eigenvaluessage: 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
extendoption.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; ifNone, an ordinary eigenvalue problem is solved (currently supported only if the base ring ofselfisRDForCDF)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; ifNone, an ordinary eigenvalue problem is solved (currently supported only if the base ring ofselfisRDForCDF)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
selfis 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 ofselfis 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(seesmith_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
- 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 asselfsuch that multiplication of elements of the base rings ofselfandrightis defined, once Sage’s coercion model is applied. If the matrices have different sizes, or if multiplication of individual entries cannot be achieved, aTypeErrorwill result.
OUTPUT:
A matrix of the same size as
selfandright. The entry in location \((i,j)\) of the output is the product of the two entries in location \((i,j)\) ofselfandright(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 methodssparse_matrix()anddense_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
selfaugmented 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
subdivideisTruethen 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 transformsselfto echelon form when it multipliesselffrom 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 ofselfthat 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