Base class for matrices, part 2#
For design documentation see matrix/docs.py.
AUTHORS:
William Stein: initial version
Jaap Spies (2006-02-24): added
prod_of_row_sums
,permanent
,permanental_minor
,rook_vector
methodsRobert Bradshaw (2007-06-14): added
subdivide
methodJaap Spies (2007-11-14): implemented
_binomial
,_choose
auxiliary functionsWilliam Stein (2007-11-18): added
_gram_schmidt_noscale
methodDavid Loeffler (2008-12-05): added
smith_form
methodDavid Loeffler (2009-06-01): added
_echelon_form_PID
methodSebastian Pancratz (2009-06-25): implemented
adjoint
andcharpoly
methods; fixedadjoint
reflecting the change that_adjoint
is now implemented inMatrix
; used the division-free algorithm forcharpoly
Rob Beezer (2009-07-13): added
elementwise_product
methodMiguel Marco (2010-06-19): modified eigenvalues and eigenvectors functions to allow the option
extend=False
Thierry Monteil (2010-10-05): bugfix for github issue #10063, so that the determinant is computed even for rings for which the
is_field
method is not implemented.Rob Beezer (2010-12-13): added
conjugate_transpose
methodRob Beezer (2011-02-05): refactored all of the matrix kernel routines; added
extended_echelon_form
,right_kernel_matrix
,QR
,_gram_schmidt_noscale
,is_similar
methodsMoritz Minzlaff (2011-03-17): corrected
_echelon_form_PID
method for matrices of one row, fixed in github issue #9053Rob Beezer (2011-06-09): added
is_normal
,is_diagonalizable
,LU
,cyclic_subspace
,zigzag_form
,rational_form
methodsRob Beezer (2012-05-27): added
indefinite_factorization
,is_positive_definite
,cholesky
methodsDarij Grinberg (2013-10-01): added first (slow) pfaffian implementation
Mario Pernici (2014-07-01): modified
rook_vector
methodRob Beezer (2015-05-25): modified
is_similar
methodSamuel Lelièvre (2020-09-18): improved method
LLL_gram
based on a patch by William Stein posted at github issue #5178, moving the method from its initial location insage.matrix.integer_matrix_dense
Michael Jung (2020-10-02): added Bär-Faddeev-LeVerrier algorithm for the Pfaffian
Moritz Firsching(2020-10-05): added
quantum_determinant
Dima Pasechnik (2022-11-08): fixed
echelonize
for inexact matrices
- class sage.matrix.matrix2.Matrix#
Bases:
Matrix
Base class for matrices, part 2
- C#
Returns the conjugate matrix.
EXAMPLES:
sage: A = matrix(QQbar, [[ -3, 5 - 3*I, 7 - 4*I], ....: [7 + 3*I, -1 + 6*I, 3 + 5*I], ....: [3 + 3*I, -3 + 6*I, 5 + I]]) sage: A.C [ -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#
Returns the conjugate-transpose (Hermitian) matrix.
EXAMPLES:
sage: A = matrix(QQbar, [[ -3, 5 - 3*I, 7 - 4*I], ....: [7 + 3*I, -1 + 6*I, 3 + 5*I], ....: [3 + 3*I, -3 + 6*I, 5 + I]]) sage: A.H [ -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)#
Return the LLL transformation matrix for this Gram matrix.
That is, the transformation matrix U over ZZ of determinant 1 that transforms the lattice with this matrix as Gram matrix to a lattice that is LLL-reduced.
Always works when
self
is positive definite, might work in some semidefinite and indefinite cases.INPUT:
self
– the Gram matrix of a quadratic form or of a lattice equipped with a bilinear formflag
– an optional flag passed toqflllgram
. According to pari:qflllgram’s documentation the options are:0
– (default), assume thatself
has 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
U
over the integers with determinant 1 such thatU.T * M * U
is LLL-reduced.ALGORITHM:
Calls PARI’s pari:qflllgram.
EXAMPLES:
Create a Gram matrix and LLL-reduce it:
sage: M = Matrix(ZZ, 2, 2, [5, 3, 3, 2]) sage: U = M.LLL_gram() sage: MM = U.transpose() * M * U sage: M, U, MM ( [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() [ 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() [-3 -1] [ 1 0] sage: Matrix(ZZ, 2, 2, [1, 0, 0, -1]).LLL_gram() [ 0 -1] [ 1 0]
However, it might fail for others, either raising a
ValueError
:sage: Matrix(ZZ, 1, 1, [0]).LLL_gram() 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() 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 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() Traceback (most recent call last): ... ValueError: qflllgram failed, perhaps the matrix is not positive definite
- LU(pivot=None, format='plu')#
Finds a decomposition into a lower-triangular matrix and an upper-triangular matrix.
INPUT:
pivot
- pivoting strategy‘auto’ (default) - see if the matrix entries are ordered (i.e. if they have an absolute value method), and if so, use a the partial pivoting strategy. Otherwise, fall back to the nonzero strategy. This is the best choice for general routines that may call this for matrix entries of a variety of types.
‘partial’ - each column is examined for the element with the largest absolute value and the row containing this element is swapped into place.
‘nonzero’ - the first nonzero element in a column is located and the row with this element is used.
format
- contents of output, see more discussion below about output.‘plu’ (default) - a triple; matrices P, L and U such that A = P*L*U.
‘compact’ - a pair; row permutation as a tuple, and the matrices L and U combined into one matrix.
OUTPUT:
Suppose that \(A\) is an \(m\times n\) matrix, then an LU decomposition is a lower-triangular \(m\times m\) matrix \(L\) with every diagonal element equal to 1, and an upper-triangular \(m\times n\) matrix, \(U\) such that the product \(LU\), after a permutation of the rows, is then equal to \(A\). For the ‘plu’ format the permutation is returned as an \(m\times m\) permutation matrix \(P\) such that
\[A = PLU\]It is more common to place the permutation matrix just to the left of \(A\). If you desire this version, then use the inverse of \(P\) which is computed most efficiently as its transpose.
If the ‘partial’ pivoting strategy is used, then the non-diagonal entries of \(L\) will be less than or equal to 1 in absolute value. The ‘nonzero’ pivot strategy may be faster, but the growth of data structures for elements of the decomposition might counteract the advantage.
By necessity, returned matrices have a base ring equal to the fraction field of the base ring of the original matrix.
In the ‘compact’ format, the first returned value is a tuple that is a permutation of the rows of \(LU\) that yields \(A\). See the doctest for how you might employ this permutation. Then the matrices \(L\) and \(U\) are merged into one matrix – remove the diagonal of ones in \(L\) and the remaining nonzero entries can replace the entries of \(U\) beneath the diagonal.
The results are cached, only in the compact format, separately for each pivot strategy called. Repeated requests for the ‘plu’ format will require just a small amount of overhead in each call to bust out the compact format to the three matrices. Since only the compact format is cached, the components of the compact format are immutable, while the components of the ‘plu’ format are regenerated, and hence are mutable.
Notice that while \(U\) is similar to row-echelon form and the rows of \(U\) span the row space of \(A\), the rows of \(U\) are not generally linearly independent. Nor are the pivot columns (or rank) immediately obvious. However for rings without specialized echelon form routines, this method is about twice as fast as the generic echelon form routine since it only acts “below the diagonal”, as would be predicted from a theoretical analysis of the algorithms.
Note
This is an exact computation, so limited to exact rings. If you need numerical results, convert the base ring to the field of real double numbers,
RDF
or the field of complex double numbers,CDF
, which will use a faster routine that is careful about numerical subtleties.ALGORITHM:
“Gaussian Elimination with Partial Pivoting,” Algorithm 21.1 of [TB1997].
EXAMPLES:
Notice the difference in the \(L\) matrix as a result of different pivoting strategies. With partial pivoting, every entry of \(L\) has absolute value 1 or less.
sage: A = matrix(QQ, [[1, -1, 0, 2, 4, 7, -1], ....: [2, -1, 0, 6, 4, 8, -2], ....: [2, 0, 1, 4, 2, 6, 0], ....: [1, 0, -1, 8, -1, -1, -3], ....: [1, 1, 2, -2, -1, 1, 3]]) sage: P, L, U = A.LU(pivot='partial') sage: P [0 0 0 0 1] [1 0 0 0 0] [0 0 0 1 0] [0 0 1 0 0] [0 1 0 0 0] sage: L [ 1 0 0 0 0] [ 1/2 1 0 0 0] [ 1/2 1/3 1 0 0] [ 1 2/3 1/5 1 0] [ 1/2 -1/3 -2/5 0 1] sage: U [ 2 -1 0 6 4 8 -2] [ 0 3/2 2 -5 -3 -3 4] [ 0 0 -5/3 20/3 -2 -4 -10/3] [ 0 0 0 0 2/5 4/5 0] [ 0 0 0 0 1/5 2/5 0] sage: A == P*L*U True sage: P, L, U = A.LU(pivot='nonzero') sage: P [1 0 0 0 0] [0 1 0 0 0] [0 0 1 0 0] [0 0 0 1 0] [0 0 0 0 1] sage: L [ 1 0 0 0 0] [ 2 1 0 0 0] [ 2 2 1 0 0] [ 1 1 -1 1 0] [ 1 2 2 0 1] sage: U [ 1 -1 0 2 4 7 -1] [ 0 1 0 2 -4 -6 0] [ 0 0 1 -4 2 4 2] [ 0 0 0 0 1 2 0] [ 0 0 0 0 -1 -2 0] sage: A == P*L*U True
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]
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
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
The base ring of the matrix may be any field, or a ring which has a fraction field implemented in Sage. The ring needs to be exact (there is a numerical LU decomposition for matrices over
RDF
andCDF
). 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: 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 [ 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 [ 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() Finite Field in a of size 5^2 sage: C == P*L*U 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)
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)
- QR(full=True)#
Returns a factorization of
self
as a unitary matrix and an upper-triangular matrix.INPUT:
full
- default:True
- ifTrue
then the returned matrices have dimensions as described below. IfFalse
theR
matrix has no zero rows and the columns ofQ
are a basis for the column space ofself
.
OUTPUT:
If
self
is an \(m\times n\) matrix andfull=True
then this method returns a pair of matrices: \(Q\) is an \(m\times m\) unitary matrix (meaning its inverse is its conjugate-transpose) and \(R\) is an \(m\times n\) upper-triangular matrix with non-negative entries on the diagonal. For a matrix of full rank this factorization is unique (due to the restriction to positive entries on the diagonal).If
full=False
then \(Q\) has \(m\) rows and the columns form an orthonormal basis for the column space 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=True
are cached, hence \(Q\) and \(R\) are immutable matrices in this case.Note
This is an exact computation, so limited to exact rings. Also the base ring needs to have a fraction field implemented in Sage and this field must contain square roots. One example is the field of algebraic numbers,
QQbar
, as used in the examples below. If you need numerical results, convert the base ring to the field of complex double numbers,CDF
, which will use a faster routine that is careful about numerical subtleties.ALGORITHM:
“Modified Gram-Schmidt,” Algorithm 8.1 of [TB1997].
EXAMPLES:
For a nonsingular matrix, the QR decomposition is unique.
sage: 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
An example with complex numbers in
QQbar
, the field of algebraic numbers.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]
A rank-deficient rectangular matrix, with both values of the
full
keyword.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: 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]
Another rank-deficient rectangular matrix, with complex entries, as a reduced decomposition.
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]
Results of full decompositions are cached and thus returned immutable.
sage: A = random_matrix(QQbar, 2, 2) sage: Q, R = A.QR() sage: Q.is_mutable() False sage: R.is_mutable() False
Trivial cases return trivial results of the correct size, and we check \(Q\) itself in one case.
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]
- T#
Returns 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]
- adjoint(*args, **kwds)#
Deprecated: Use
adjugate()
instead. See github issue #10501 for details.
- adjoint_classical()#
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 [ 41/10 -1/28] [-33/13 5/3] sage: M * N [7363/1092 0] [ 0 7363/1092]
An alias is
adjoint_classical()
, which replaces the deprecatedadjoint()
method:sage: M.adjoint() ...: 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() [ 41/10 -1/28] [-33/13 5/3]
ALGORITHM:
Use PARI whenever the method
self._adjugate
is included to do so in an inheriting class. Otherwise, use a generic division-free algorithm that computes the adjugate matrix from the characteristic polynomial.The result is cached.
- adjugate()#
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 [ 41/10 -1/28] [-33/13 5/3] sage: M * N [7363/1092 0] [ 0 7363/1092]
An alias is
adjoint_classical()
, which replaces the deprecatedadjoint()
method:sage: M.adjoint() ...: 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() [ 41/10 -1/28] [-33/13 5/3]
ALGORITHM:
Use PARI whenever the method
self._adjugate
is included to do so in an inheriting class. Otherwise, use a generic division-free algorithm that computes the adjugate matrix from the characteristic polynomial.The result is cached.
- apply_map(phi, R=None, sparse=None)#
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
– True to make the output a sparse matrix; default Falsephi
- 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) sage: f = lambda x: k(x) sage: n = m.apply_map(f); n [0 1 2] [0 1 2] [0 1 2] sage: n.parent() 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 [0 1 2] [0 1 2] [0 1 2] sage: n.parent() Full MatrixSpace of 3 by 3 dense matrices over Finite Field in a of size 3^2
If self is subdivided, the result will be as well:
sage: m = matrix(2, 2, srange(4)) sage: m.subdivide(None, 1); m [0|1] [2|3] sage: m.apply_map(lambda x: x*x) [0|1] [4|9]
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
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
- apply_morphism(phi)#
Apply the morphism phi to the coefficients of this dense matrix.
The resulting matrix is over the codomain of phi.
INPUT:
phi
- a morphism, so phi is callable and phi.domain() and phi.codomain() are defined. The codomain must be a ring.
OUTPUT: a matrix over the codomain of phi
EXAMPLES:
sage: m = matrix(ZZ, 3, 3, range(9)) sage: phi = ZZ.hom(GF(5)) sage: m.apply_morphism(phi) [0 1 2] [3 4 0] [1 2 3] sage: parent(m.apply_morphism(phi)) Full MatrixSpace of 3 by 3 dense matrices over Finite Field of size 5
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]
- as_bipartite_graph()#
Construct a bipartite graph
B
representing the matrix uniquely.Vertices are labeled 1 to
nrows
on the left andnrows
+ 1 tonrows
+ncols
on the right, representing rows and columns correspondingly. Each row is connected to each column with an edge weighted by the value of the corresponding matrix entry.This graph is a helper for calculating automorphisms of a matrix under row and column permutations. See
automorphisms_of_rows_and_columns()
.OUTPUT:
A bipartite graph.
EXAMPLES:
sage: M = matrix(QQ, [[1/3, 7], [6, 1/4], [8, -5]]) sage: M [1/3 7] [ 6 1/4] [ 8 -5] sage: B = M.as_bipartite_graph() sage: 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
- as_sum_of_permutations()#
Returns 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: decomp = sage.combinat.permutation.bistochastic_as_sum_of_permutations(M) sage: print(decomp) 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 = sage.combinat.permutation.bistochastic_as_sum_of_permutations(M) Traceback (most recent call last): ... ValueError: The matrix is not bistochastic
- automorphisms_of_rows_and_columns()#
Return the automorphisms of
self
under permutations of rows and columns as a list of pairs ofPermutationGroupElement
objects.EXAMPLES:
sage: M = matrix(ZZ,[[1,0],[1,0],[0,1]]) sage: M [1 0] [1 0] [0 1] sage: A = M.automorphisms_of_rows_and_columns() sage: 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
One can now apply these automorphisms to
M
to show that it leaves it invariant:sage: all(M.with_permuted_rows_and_columns(*i) == M for i in A) True
Check that github 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() [((), ()), ((1,3)(2,5), (1,3)(2,5))]
- block_ldlt(classical=False)#
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
– (default:False
) whether or not to attempt a classical non-block \(LDL^{T}\) factorization with no row/column swaps.
Warning
Not all matrices have a classical \(LDL^{T}\) factorization. Set
classical=True
at your own risk, preferably after verifying that your matrix is positive-definite and (over inexact rings) not ill-conditioned.OUTPUT:
If the input matrix is not Hermitian, the output from this function is undefined. Otherwise, we return a triple \((P,L,D)\) such that \(A = PLDL^{*}P^{T}\) and
\(P\) is a permutation matrix,
\(L\) is unit lower-triangular,
\(D\) is a block-diagonal matrix whose blocks are of size one or two.
With
classical=True
, the permutation matrix \(P\) is always an identity matrix and the diagonal blocks are always one-by-one. AValueError
is raised if the matrix has no classical \(LDL^{T}\) factorization.ALGORITHM:
We essentially follow “Algorithm A” in the paper by Bunch and Kaufman [BK1977] that describes the stable pivoting strategy. The same scheme is described by Higham [Hig2002].
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
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] )
The same is true of the following complex Hermitian matrix:
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] )
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() ( [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: 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
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] )
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] )
An indefinite Hermitian matrix that happens to have a classical factorization:
sage: F.<I> = QuadraticField(-1) sage: A = matrix(F, [[ 2, 4 - 2*I, 2 + 2*I], ....: [4 + 2*I, 8, 10*I], ....: [2 - 2*I, -10*I, -3]]) sage: A.block_ldlt(classical=True)[1:] ( [ 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)#
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') T^2 - 5*T - 2
- charpoly(var='x', algorithm=None)#
Returns 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:"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]
An example over \(\QQ\):
sage: A = MatrixSpace(QQ,3)(range(9)) sage: A.charpoly('x') x^3 - 12*x^2 - 18*x sage: A.trace() 12 sage: 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
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
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
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
Here is an example over a 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
Over integers modulo \(n\) with composite \(n\):
sage: A = Mat(Integers(6),3,3)(range(9)) sage: A.charpoly() x^3
Here is an example over a general commutative ring, that is to say, as of version 4.0.2, Sage does not even positively determine that
S
in the following example is an integral domain. But the computation of the characteristic polynomial succeeds as follows:sage: 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]]) sage: 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
- cholesky()#
Returns the Cholesky decomposition of a Hermitian matrix.
INPUT:
A positive-definite matrix. Generally, the base ring for the entries of the matrix needs to be a subfield of the algebraic numbers (
QQbar
). Examples include the rational numbers (QQ
), some number fields, and real algebraic numbers and the algebraic numbers themselves. Symbolic matrices can also occasionally be factored.OUTPUT:
For a matrix \(A\) the routine returns a lower triangular matrix \(L\) such that,
\[A = LL^\ast\]where \(L^\ast\) is the conjugate-transpose. If the matrix is not positive-definite (for example, if it is not Hermitian) then a
ValueError
results.If possible, the output matrix will be over the fraction field of the base ring of the input matrix. If that fraction field is missing the requisite square roots but if no imaginaries are encountered, then the algebraic-reals will be used. Otherwise, the algebraic closure of the fraction field (typically
QQbar
) will be used.ALGORITHM:
First we ensure that the matrix \(A\)
is_hermitian()
. Afterwards, we attempt to compute a 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() sage: 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
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() sage: L [ 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() Full MatrixSpace of 4 by 4 dense matrices over Algebraic Real Field sage: L*L.transpose() == A 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: 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() sage: 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
The field of algebraic numbers is an ideal setting for this computation:
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
Results are cached, hence immutable. Use the
copy
function if you need to make a change:sage: A = matrix(QQ, [[ 4, -2, 4, 2], ....: [-2, 10, -2, -7], ....: [ 4, -2, 8, 4], ....: [ 2, -7, 4, 7]]) sage: L = A.cholesky() sage: L.is_immutable() True sage: from copy import copy sage: LC = copy(L) sage: LC[0,0] = 1000 sage: LC [1000 0 0 0] [ -1 3 0 0] [ 2 0 2 0] [ 1 -2 1 1]
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 True
Even symbolic matrices can sometimes be factored:
sage: A = matrix(SR, [[pi,0], [0,pi]]) # optional - sage.symbolic sage: A.cholesky() # optional - 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'
The matrix may not be Hermitian:
sage: F.<a> = FiniteField(5^4) sage: A = matrix(F, [[2+a^3, 3], [3, 3]]) sage: A.cholesky() Traceback (most recent call last): ... ValueError: matrix is not Hermitian
The matrix may not be positive-definite:
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
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
- column_module()#
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]
- column_space()#
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: 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]
- conjugate()#
Return the conjugate of self, i.e. the matrix whose entries are the conjugates of the entries of self.
EXAMPLES:
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]
A matrix over a not-totally-real number field:
sage: x = polygen(ZZ, 'x') sage: K.<j> = NumberField(x^2 + 5) # optional - sage.rings.number_field sage: M = matrix(K, [[1+j,1], [0,2*j]]) # optional - sage.rings.number_field sage: M.conjugate() # optional - sage.rings.number_field [-j + 1 1] [ 0 -2*j]
There is a shortcut for the conjugate:
sage: M.C [-j + 1 1] [ 0 -2*j]
There is also a shortcut for the conjugate transpose, or “Hermitian transpose”:
sage: M.H [-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
- conjugate_transpose()#
Return the transpose of
self
after each entry has been converted to its complex conjugate.Note
This function is sometimes known as the “adjoint” of a matrix, though there is substantial variation and some confusion with the use of that term.
OUTPUT:
A matrix formed by taking the complex conjugate of every entry of
self
and then transposing the resulting matrix.Complex conjugation is implemented for many subfields of the complex numbers. See the examples below, or more at
conjugate()
.EXAMPLES:
sage: M = matrix(SR, 2, 2, [[2-I, 3+4*I], [9-6*I, 5*I]]) # optional - sage.symbolic sage: M.base_ring() # optional - sage.symbolic Symbolic Ring sage: M.conjugate_transpose() # optional - sage.symbolic [ I + 2 6*I + 9] [-4*I + 3 -5*I] 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]
There is also a shortcut for the conjugate transpose, or “Hermitian transpose”:
sage: M.H [ 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
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: 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]
Furthermore, this method can be applied to matrices over quadratic extensions of finite fields:
sage: F.<a> = GF(9,'a') sage: N = matrix(F, 2, [0,a,-a,1]); N [ 0 a] [2*a 1] sage: N.conjugate_transpose() [ 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'
- cyclic_subspace(v, var=None, basis='echelon')#
Create a cyclic subspace for a vector, and optionally, a minimal polynomial for the iterated powers.
These subspaces are also known as Krylov subspaces. They are spanned by the vectors
\[\{v, Av, A^2v, A^3v, \dots \}\]INPUT:
self
- a square matrix with entries from a field.v
- a vector with a degree equal to the size of the matrix and entries compatible with the entries of the matrix.var
- default:None
- if specified as a string or a generator of a polynomial ring, then this will be used to construct a polynomial reflecting a relation of linear dependence on the powers \(A^iv\) and this will cause the polynomial to be returned along with the subspace. A generator must create polynomials with coefficients from the same field as the matrix entries.basis
- default:echelon
- the basis for the subspace is “echelonized” by default, but the keyword ‘iterates’ will return a subspace with a user basis equal to the largest linearly independent set \(\{v, Av, A^2v, A^3v, \dots, A^{k-1}v \}\).
OUTPUT:
Suppose \(k\) is the smallest power such that \(\{v, Av, A^2v, A^3v, \dots, A^{k}v \}\) is linearly dependent. Then the subspace returned will have dimension \(k\) and be spanned by the powers \(0\) through \(k-1\).
If a polynomial is requested through the use of the
var
keyword, then a pair is returned, with the polynomial first and the subspace second. The polynomial is the unique monic polynomial whose coefficients provide a relation of linear dependence on the first \(k\) powers.For less convenient, but more flexible output, see the helper method “_cyclic_subspace” in this module.
EXAMPLES:
sage: A = matrix(QQ, [[5,4,2,1],[0,1,-1,-1],[-1,-1,3,0],[1,1,-1,2]]) sage: v = vector(QQ, [0,1,0,0]) sage: E = A.cyclic_subspace(v); E Vector space of degree 4 and dimension 3 over Rational Field Basis matrix: [ 1 0 0 0] [ 0 1 0 0] [ 0 0 1 -1] sage: F = A.cyclic_subspace(v, basis='iterates'); F Vector space of degree 4 and dimension 3 over Rational Field User basis matrix: [ 0 1 0 0] [ 4 1 -1 1] [23 1 -8 8] sage: E == F True sage: p, S = A.cyclic_subspace(v, var='T'); p T^3 - 9*T^2 + 24*T - 16 sage: p.degree() == E.dimension() True
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
- decomposition(algorithm='spin', is_diagonalizable=False, dual=False)#
Returns 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
- ‘spin’ (default): algorithm involves iterating the action of self on a vector. ‘kernel’: naively just compute \(ker(f_i(A))\) for each factor \(f_i\).dual
- bool (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 bool, and t is True exactly when the charpoly of self on V is irreducible.(optional) list - list of pairs (W,t), where W is a vector space and t is a bool, and t is True exactly when the charpoly of the transpose of self on W is irreducible.
EXAMPLES:
sage: A = matrix(ZZ, 4, [3,4,5,6,7,3,8,10,14,5,6,7,2,2,10,9]) sage: B = matrix(QQ, 6, 6, range(36)) sage: B*11 [ 0 11 22 33 44 55] [ 66 77 88 99 110 121] [132 143 154 165 176 187] [198 209 220 231 242 253] [264 275 286 297 308 319] [330 341 352 363 374 385] sage: A.decomposition() [ (Ambient free module of rank 4 over the principal ideal domain Integer Ring, True) ] sage: B.decomposition() [ (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)#
Suppose the right action of self on M leaves M invariant. Return the decomposition of M as a list of pairs (W, is_irred) where is_irred is True if the charpoly of self acting on the factor W is irreducible.
Additional inputs besides M are passed onto the decomposition command.
INPUT:
\(M\) – A subspace of the free module
self
acts on.check_restrict
– A boolean (default:True
); Call restrictwith or without check.
kwds
– Keywords that will be forwarded todecomposition()
.
EXAMPLES:
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]
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) [ (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()#
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
A trivial example:
sage: A = matrix(QQ, 0,2) sage: 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
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
Here’s an example involving a cyclotomic 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
- density()#
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
sage: A = matrix(QQ,3,3,[0,1,2,3,0,0,6,7,8]) sage: A.density() 2/3
sage: a = matrix([[],[],[],[]]) sage: a.density() 0
- derivative(*args)#
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: 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)
- det(*args, **kwds)#
Synonym for self.determinant(…).
EXAMPLES:
sage: A = MatrixSpace(Integers(8),3)([1,7,3, 1,1,1, 3,4,5]) sage: A.det() 6
- determinant(algorithm=None)#
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:"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
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
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
A matrix over a non-domain:
sage: m = matrix(Integers(4), 2, [1,2,2,3]) sage: m.determinant() 3
- diagonal()#
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]
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]
Empty matrices behave properly.
sage: E = matrix(0, 5, []); E [] sage: E.diagonal() []
- diagonalization(base_field=None)#
Return a diagonal matrix similar to
self
along with the transformation matrix.INPUT:
base_field
– if given,self
is regarded as a matrix over it
OUTPUT: a diagonal matrix \(D\) and an invertible matrix \(P\) such that \(P^{-1}AP=D\), if
self
is a diagonalizable matrix \(A\).EXAMPLES:
sage: 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: A = matrix(QQ, 2, [0, 2, 1, 0]) sage: A.is_diagonalizable() False sage: A.is_diagonalizable(QQbar) True sage: D, P = A.diagonalization(QQbar) sage: P^-1*A*P == D True
Matrices may fail to be diagonalizable for various reasons:
sage: A = matrix(QQ, 2, [1,2,3,4,5,6]) sage: 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]) sage: 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.]) sage: 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]) sage: D [0 2] [1 0] sage: D.diagonalization() Traceback (most recent call last): ... ValueError: not diagonalizable over Rational Field sage: E = matrix(QQ, 2, [3, 1, 0, 3]) sage: E [3 1] [0 3] sage: E.diagonalization() Traceback (most recent call last): ... ValueError: not diagonalizable sage: E.jordan_form() [3 1] [0 3]
- echelon_form(algorithm='default', cutoff=0, **kwds)#
Return the echelon form of
self
.Note
This row reduction does not use division if the matrix is not over a field (e.g., if the matrix is over the integers). If you want to calculate the echelon form using division, then use
rref()
, which assumes that the matrix entries are in a field (specifically, the field of fractions of the base ring of the matrix).INPUT:
algorithm
– string. Which algorithm to use. Choices are'default'
: Let Sage choose an algorithm (default).'classical'
: Gauss elimination.'partial_pivoting'
: Gauss elimination, using partial pivoting (if base ring has absolute value)'scaled_partial_pivoting'
: Gauss elimination, using scaled partial pivoting (if base ring has absolute value)'scaled_partial_pivoting_valuation'
: Gauss elimination, using scaled partial pivoting (if base ring has valuation)'strassen'
: use a Strassen divide and conquer algorithm (if available)
cutoff
– integer. Only used if the Strassen algorithm is selected.transformation
– boolean. Whether to also return the transformation matrix. Some matrix backends do not provide this information, in which case this option is ignored.
OUTPUT:
The reduced row echelon form of
self
, as an immutable matrix. Note that self is not changed by this command. Useechelonize()
to changeself
in place.If the optional parameter
transformation=True
is specified, the output consists of a pair \((E,T)\) of matrices where \(E\) is the echelon form ofself
and \(T\) is the transformation matrix.EXAMPLES:
sage: MS = MatrixSpace(GF(19),2,3) sage: C = MS.matrix([1,2,3,4,5,6]) sage: C.rank() 2 sage: C.nullity() 0 sage: C.echelon_form() [ 1 0 18] [ 0 1 2]
The matrix library used for \(\ZZ/p\)-matrices does not return the transformation matrix, so the
transformation
option is ignored:sage: C.echelon_form(transformation=True) [ 1 0 18] [ 0 1 2] sage: D = matrix(ZZ, 2, 3, [1,2,3,4,5,6]) sage: D.echelon_form(transformation=True) ( [1 2 3] [ 1 0] [0 3 6], [ 4 -1] ) sage: E, T = D.echelon_form(transformation=True) sage: T*D == E True
- echelonize(algorithm='default', cutoff=0, **kwds)#
Transform
self
into a matrix in echelon form over the same base ring as self.Note
This row reduction does not use division if the matrix is not over a field (e.g., if the matrix is over the integers). If you want to calculate the echelon form using division, then use
rref()
, which assumes that the matrix entries are in a field (specifically, the field of fractions of the base ring of the matrix).INPUT:
algorithm
– string. Which algorithm to use. Choices are'default'
: Let Sage choose an algorithm (default).'classical'
: Gauss elimination.'partial_pivoting'
: Gauss elimination, using partial pivoting (if base ring has absolute value)'scaled_partial_pivoting'
: Gauss elimination, using scaled partial pivoting (if base ring has absolute value)'scaled_partial_pivoting_valuation'
: Gauss elimination, using scaled partial pivoting (if base ring has valuation)'strassen'
: use a Strassen divide and conquer algorithm (if available)
cutoff
– integer. Only used if the Strassen algorithm is selected.transformation
– boolean. Whether to also return the transformation matrix. Some matrix backends do not provide this information, in which case this option is ignored.
OUTPUT:
The matrix
self
is put into echelon form. Nothing is returned unless the keyword optiontransformation=True
is specified, in which case the transformation matrix is returned.EXAMPLES:
sage: a = matrix(QQ,3,3,range(9)); a [0 1 2] [3 4 5] [6 7 8] sage: a.echelonize() sage: a [ 1 0 -1] [ 0 1 2] [ 0 0 0]
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]
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]
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? [x y] [x y]
sage: b = a.change_ring(R.fraction_field()) sage: b.echelon_form() # potentially useful [ 1 y/x] [ 0 0]
We check that the echelon form works for matrices over p-adics. See github issue #17272:
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)]
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'.
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]
The transformation matrix is optionally returned:
sage: m_original = m sage: transformation_matrix = m.echelonize(transformation=True) sage: m == transformation_matrix * m_original True
- eigenmatrix_left(other=None)#
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
self
is a square matrix \(A\), then the output is a diagonal matrix \(D\) and a matrix \(P\) such that\[P A = D P,\]where the rows of \(P\) are eigenvectors of \(A\) and the diagonal entries of \(D\) are the corresponding eigenvalues.
If a matrix \(B\) is passed as optional argument, the output is a solution to the generalized eigenvalue problem such that
\[P A = D P B.\]The ordinary eigenvalue problem is equivalent to the generalized one if \(B\) is the identity matrix.
The generalized eigenvector decomposition is currently only implemented for matrices over
RDF
andCDF
.EXAMPLES:
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
Because \(P\) is invertible, \(A\) is diagonalizable.
sage: A == (~P)*D*P 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: 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
A generalized eigenvector decomposition:
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
The matrix \(B\) in a generalized eigenvalue problem may be singular:
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]
In this case, we can still verify the eigenvector equation for the first eigenvalue and first eigenvector:
sage: l = D[0, 0] sage: v = P[0, :] sage: (v * A - l * v * B).norm() < 1e-14 True
The second eigenvector is contained in the left kernel of \(B\):
sage: (P[1, :] * B).norm() < 1e-14 True
- eigenmatrix_right(other=None)#
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
self
is a square matrix \(A\), then the output is a diagonal matrix \(D\) and a matrix \(P\) such that\[A P = P D,\]where the columns of \(P\) are eigenvectors of \(A\) and the diagonal entries of \(D\) are the corresponding eigenvalues.
If a matrix \(B\) is passed as optional argument, the output is a solution to the generalized eigenvalue problem such that
\[A P = B P D.\]The ordinary eigenvalue problem is equivalent to the generalized one if \(B\) is the identity matrix.
The generalized eigenvector decomposition is currently only implemented for matrices over
RDF
andCDF
.EXAMPLES:
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
Because \(P\) is invertible, \(A\) is diagonalizable.
sage: A == P*D*(~P) 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: 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
A generalized eigenvector decomposition:
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
The matrix \(B\) in a generalized eigenvalue problem may be singular:
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]
In this case, we can still verify the eigenvector equation for the first eigenvalue and first eigenvector:
sage: l = D[0, 0] sage: v = P[:, 0] sage: (A * v - B * v * l).norm() < 1e-14 True
The second eigenvector is contained in the right kernel of \(B\):
sage: (B * P[:, 1]).norm() < 1e-14 True
- eigenspaces_left(format='all', var='a', algebraic_multiplicity=False)#
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 classes.format
- default:None
'all'
- attempts to create every eigenspace. This will always be possible for matrices with rational entries.'galois'
- for each irreducible factor of the characteristic polynomial, a single eigenspace will be output for a single root/eigenvalue for the irreducible factor.None
- Uses the ‘all’ format if the base ring is contained in an algebraically closed field which is implemented. Otherwise, uses the ‘galois’ format.
var
- default: ‘a’ - variable name used to represent elements of the root field of each irreducible factor of the characteristic polynomial. If var=’a’, then the root fields will be in terms of a0, a1, a2, …., where the numbering runs across all the irreducible factors of the characteristic polynomial, even for linear factors.algebraic_multiplicity
- default: False - whether or not to include the algebraic multiplicity of each eigenvalue in the output. See the discussion below.
OUTPUT:
If algebraic_multiplicity=False, return a list of pairs (e, V) where e is an eigenvalue of the matrix, and V is the corresponding left eigenspace. For Galois conjugates of eigenvalues, there may be just one representative eigenspace, depending on the
format
keyword.If algebraic_multiplicity=True, return a list of triples (e, V, n) where e and V are as above and n is the algebraic multiplicity of the eigenvalue.
Warning
Uses a somewhat naive algorithm (simply factors the characteristic polynomial and computes kernels directly over the extension field).
EXAMPLES:
We compute the left eigenspaces of a \(3\times 3\) rational matrix. First, we request \(all\) of the eigenvalues, so the results are in the field of algebraic numbers, \(QQbar\). Then we request just one eigenspace per irreducible factor of the characteristic polynomial with the \(galois\) keyword.
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 [ (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: 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
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') [ (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: 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) ]
Next we compute the left eigenspaces over the finite field of order 11.
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]) ]
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: 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: 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]
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: 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') 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') [ (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)#
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 classes.format
- default:None
'all'
- attempts to create every eigenspace. This will always be possible for matrices with rational entries.'galois'
- for each irreducible factor of the characteristic polynomial, a single eigenspace will be output for a single root/eigenvalue for the irreducible factor.None
- Uses the ‘all’ format if the base ring is contained in an algebraically closed field which is implemented. Otherwise, uses the ‘galois’ format.
var
- default: ‘a’ - variable name used to represent elements of the root field of each irreducible factor of the characteristic polynomial. If var=’a’, then the root fields will be in terms of a0, a1, a2, …., where the numbering runs across all the irreducible factors of the characteristic polynomial, even for linear factors.algebraic_multiplicity
- default: False - whether or not to include the algebraic multiplicity of each eigenvalue in the output. See the discussion below.
OUTPUT:
If algebraic_multiplicity=False, return a list of pairs (e, V) where e is an eigenvalue of the matrix, and V is the corresponding left eigenspace. For Galois conjugates of eigenvalues, there may be just one representative eigenspace, depending on the
format
keyword.If algebraic_multiplicity=True, return a list of triples (e, V, n) where e and V are as above and n is the algebraic multiplicity of the eigenvalue.
Warning
Uses a somewhat naive algorithm (simply factors the characteristic polynomial and computes kernels directly over the extension field).
EXAMPLES:
Right eigenspaces are computed from the left eigenspaces of the transpose of the matrix. As such, there is a greater collection of illustrative examples at the
eigenspaces_left()
.We compute the right eigenspaces of a \(3\times 3\) rational matrix.
sage: 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
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') [ (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: 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: 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]
- eigenvalue_multiplicity(s)#
Return the multiplicity of
s
as a generalized eigenvalue of the matrix.EXAMPLES:
sage: M = Matrix(QQ, [[0,1],[0,0]]) sage: M.eigenvalue_multiplicity(0) 2 sage: M.eigenvalue_multiplicity(1) 0 sage: M = posets.DiamondPoset(5).coxeter_transformation() sage: [M.eigenvalue_multiplicity(x) for x in [-1, 1]] [3, 2]
- eigenvalues(extend=True)#
Return a sequence of the eigenvalues of a matrix, with multiplicity. If the eigenvalues are roots of polynomials in QQ, then QQbar elements are returned that represent each separate root.
If the option extend is set to False, only eigenvalues in the base ring are considered.
EXAMPLES:
sage: a = matrix(ZZ, 4, range(16)); a [ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15] sage: sorted(a.eigenvalues(), reverse=True) [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() [-0.9386318578049146?, 15.50655435353258?, 0.2160387521361705? - 4.713151979747493?*I, 0.2160387521361705? + 4.713151979747493?*I]
A symmetric matrix a+a.transpose() should have real eigenvalues
sage: b=a+a.transpose() sage: ev = b.eigenvalues(); ev [-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 -8.35066086057957? sage: p = e.minpoly(); p x^4 - 30*x^3 - 171*x^2 + 1460*x + 1784 sage: p(e) == 0 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() (Number Field in a with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264, a + 7, Ring morphism: From: Number Field in a with defining polynomial y^4 - 2*y^3 - 507*y^2 - 3972*y - 4264 To: Algebraic Real Field Defn: a |--> -15.35066086057957?)
Notice the effect of the extend option.
sage: M=matrix(QQ,[[0,-1,0],[1,0,0],[0,0,2]]) sage: M.eigenvalues() [2, -1*I, 1*I] sage: M.eigenvalues(extend=False) [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 [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] sage: e.parent() Algebraic closure of Finite Field of size 3 sage: e.as_finite_field_element() (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)#
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 ofself
isRDF
orCDF
)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: 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
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() [(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) [(2, [ (0, 0, 1) ], 1)]
- eigenvectors_right(other=None, extend=True)#
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 ofself
isRDF
orCDF
)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: 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
- elementary_divisors()#
If self is a matrix over a principal ideal domain R, return elements \(d_i\) for \(1 \le i \le k = \min(r,s)\) where \(r\) and \(s\) are the number of rows and columns of self, such that the cokernel of self is isomorphic to
\[R/(d_1) \oplus R/(d_2) \oplus R/(d_k)\]with \(d_i \mid d_{i+1}\) for all \(i\). These are the diagonal entries of the Smith form of self (see
smith_form()
).EXAMPLES:
sage: x = polygen(ZZ, 'x') sage: OE.<w> = EquationOrder(x^2 - x + 2) # optional - sage.rings.number_field sage: m = Matrix([[1, w], [w, 7]]) # optional - sage.rings.number_field sage: m.elementary_divisors() # optional - sage.rings.number_field [1, -w + 9]
See also
- elementwise_product(right)#
Returns 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 asself
such that multiplication of elements of the base rings ofself
andright
is defined, once Sage’s coercion model is applied. If the matrices have different sizes, or if multiplication of individual entries cannot be achieved, aTypeError
will result.
OUTPUT:
A matrix of the same size as
self
andright
. The entry in location \((i,j)\) of the output is the product of the two entries in location \((i,j)\) ofself
andright
(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
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
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
Non-commutative rings behave as expected. These are the usual quaternions.
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]
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
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'
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'
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
- exp()#
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: 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()) <class 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'> sage: a=matrix([[1/2,2/3],[3/4,4/5]]) sage: a.exp() [-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]]) sage: a.exp() [ 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 1e-14 [42748127.31532951 7368259.244159399] [234538976.1381042 40426191.45156228]
- extended_echelon_form(subdivide=False, **kwds)#
Returns the echelon form of
self
augmented with an identity matrix.INPUT:
subdivide
- default:False
- determines if the returned matrix is subdivided. See the description of the (output) below for details.kwds
- additional keywords that can be passed to the method that computes the echelon form.
OUTPUT:
If \(A\) is an \(m\times n\) matrix, add the \(m\) columns of an \(m\times m\) identity matrix to the right of
self
. Then row-reduce this \(m\times(n+m)\) matrix. This matrix is returned as an immutable matrix.If
subdivide
isTrue
then the returned matrix has a single division among the columns and a single division among the rows. The column subdivision has \(n\) columns to the left and \(m\) columns to the right. The row division separates the non-zero 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 transformsself
to echelon form when it multipliesself
from the left. When the base ring is a field, the uniqueness of reduced row-echelon form implies that this transformation matrix can be taken as the coefficients giving a canonical set of linear combinations of the rows ofself
that yield reduced row-echelon form.When subdivided as described above, and again over a field, the parts of the subdivision in the upper-left corner and lower-right corner satisfy several interesting relationships with the row space, column space, left kernel and right kernel of
self
. See the examples below.Note
This method returns an echelon form. If the base ring is not a field, no attempt is made to move to the fraction field. See an example below where the base ring is changed manually.
EXAMPLES:
The four relationships at the end of this example hold in general.
sage: A = matrix(QQ, [[2, -1, 7, -1, 0, -3], ....: [-1, 1, -5, 3, 4, 4], ....: [2, -1, 7, 0, 2, -2], ....: [2, 0, 4, 3, 6, 1], ....: [2, -1, 7, 0, 2, -2]]) sage: E = A.extended_echelon_form(subdivide=True); E [ 1 0 2 0 0 -1| 0 -1 0 1 -1] [ 0 1 -3 0 -2 0| 0 -2 0 2 -3] [ 0 0 0 1 2 1| 0 2/3 0 -1/3 2/3] [-----------------------------+------------------------] [ 0 0 0 0 0 0| 1 2/3 0 -1/3 -1/3] [ 0 0 0 0 0 0| 0 0 1 0 -1] sage: J = E.matrix_from_columns(range(6,11)); J [ 0 -1 0 1 -1] [ 0 -2 0 2 -3] [ 0 2/3 0 -1/3 2/3] [ 1 2/3 0 -1/3 -1/3] [ 0 0 1 0 -1] sage: J*A == A.rref() True sage: C = E.subdivision(0,0); C [ 1 0 2 0 0 -1] [ 0 1 -3 0 -2 0] [ 0 0 0 1 2 1] sage: L = E.subdivision(1,1); L [ 1 2/3 0 -1/3 -1/3] [ 0 0 1 0 -1] sage: A.right_kernel() == C.right_kernel() True sage: A.row_space() == C.row_space() True sage: A.column_space() == L.right_kernel() True sage: A.left_kernel() == L.row_space() True
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
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]
Subdivided, or not, the result is immutable, so make a copy if you want to make changes.
sage: A = matrix(FiniteField(7), [[2,0,3], [5,5,3], [5,6,5]]) sage: E = A.extended_echelon_form() sage: E.is_mutable() False sage: F = A.extended_echelon_form(subdivide=True) sage: F [1 0 0|0 4 6] [0 1 0|4 2 2] [0 0 1|5 2 3] [-----+-----] sage: F.is_mutable() False sage: G = copy(F) sage: G.subdivide([],[]); G [1 0 0 0 4 6] [0 1 0 4 2 2] [0 0 1 5 2 3]
If you want to determine exactly which algorithm is used to compute the echelon form, you can add additional keywords to pass on to the
echelon_form()
routine employed on the augmented matrix. Sending the flaginclude_zero_rows
is a bit silly, since the extended echelon form will never have any zero rows.sage: A = matrix(ZZ, [[1,2], [5,0], [5,9]]) sage: E = A.extended_echelon_form(algorithm='padic', include_zero_rows=False) sage: E [ 1 0 36 1 -8] [ 0 1 5 0 -1] [ 0 0 45 1 -10]
- fcp(var='x')#
Return the factorization of the characteristic polynomial of self.
INPUT:
var
- (default: ‘x’) name of variable of charpoly
EXAMPLES:
sage: M = MatrixSpace(QQ,3,3) sage: A = M([1,9,-7,4/5,4,3,6,4,3]) sage: A.fcp() x^3 - 8*x^2 + 209/5*x - 286 sage: A = M([3, 0, -2, 0, -2, 0, 0, 0, 0]) sage: A.fcp('T') (T - 3) * T * (T + 2)
- find(f, indices=False)#
Find elements in this matrix satisfying the constraints in the function \(f\). The function is evaluated on each element of the matrix .
INPUT:
f
- a function that is evaluated on each element of this matrix.indices
- whether or not to return the indices and elements of this matrix that satisfy the function.
OUTPUT: If
indices
is not specified, return a matrix with 1 where \(f\) is satisfied and 0 where it is not. Ifindices
is specified, return a dictionary containing the elements of this matrix satisfying \(f\).EXAMPLES:
sage: M = matrix(4,3,[1, -1/2, -1, 1, -1, -1/2, -1, 0, 0, 2, 0, 1]) sage: M.find(lambda entry:entry==0) [0 0 0] [0 0 0] [0 1 1] [0 1 0]
sage: M.find(lambda u:u<0) [0 1 1] [0 1 1] [1 0 0] [0 0 0]
sage: M = matrix(4,3,[1, -1/2, -1, 1, -1, -1/2, -1, 0, 0, 2, 0, 1]) sage: len(M.find(lambda u:u<1 and u>-1,indices=True)) 5
sage: M.find(lambda u:u!=1/2) [1 1 1] [1 1 1] [1 1 1] [1 1 1]
sage: M.find(lambda u:u>1.2) [0 0 0] [0 0 0] [0 0 0] [1 0 0]
sage: sorted(M.find(lambda u:u!=0,indices=True).keys()) == M.nonzero_positions() True
- get_subdivisions()#
Returns the current subdivision of self.
EXAMPLES:
sage: M = matrix(5, 5, range(25)) sage: M.subdivisions() ([], []) sage: M.subdivide(2,3) sage: M.subdivisions() ([2], [3]) sage: N = M.parent()(1) sage: N.subdivide(M.subdivisions()); N [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]
- gram_schmidt(orthonormal=False)#
Performs Gram-Schmidt orthogonalization on the rows of the matrix, returning a new matrix and a matrix accomplishing the transformation.
INPUT:
self
- a matrix whose rows are to be orthogonalized.orthonormal
- default:False
- ifTrue
the returned orthogonal vectors are unit vectors. This keyword is ignored if the matrix is overRDF
orCDF
and the results are always orthonormal.
OUTPUT:
A pair of matrices,
G
andM
such that ifA
representsself
, where the parenthetical properties occur whenorthonormal = True
:A = M*G
The rows of
G
are an orthogonal (resp. orthonormal) set of vectors.G
times the conjugate-transpose ofG
is a diagonal (resp. identity) matrix.The row space of
G
equals the row space ofA
.M
is a full-rank matrix with zeros above the diagonal.
For exact rings, any zero vectors produced (when the original vectors are linearly dependent) are not output, thus the orthonormal set is linearly independent, and thus a basis for the row space of the original matrix.
Any notion of a Gram-Schmidt procedure requires that the base ring of the matrix has a fraction field implemented. In order to arrive at an orthonormal set, it must be possible to construct square roots of the elements of the base field. In Sage, your best option is the field of algebraic numbers,
QQbar
, which properly contains the rationals and number fields.If you have an approximate numerical matrix, then this routine requires that your base field be the real and complex double-precision floating point numbers,
RDF
andCDF
. In this case, the matrix is treated as having full rank, as no attempt is made to recognize linear dependence with approximate calculations.EXAMPLES:
Inexact Rings, Numerical Matrices:
First, the inexact rings,
CDF
andRDF
.sage: A = matrix(CDF, [[ 0.6454 + 0.7491*I, -0.8662 + 0.1489*I, 0.7656 - 0.00344*I], ....: [-0.2913 + 0.8057*I, 0.8321 + 0.8170*I, -0.6744 + 0.9248*I], ....: [ 0.2554 + 0.3517*I, -0.4454 - 0.1715*I, 0.8325 - 0.6282*I]]) sage: G, M = A.gram_schmidt() sage: G.round(6) # random signs [-0.422243 - 0.490087*I 0.566698 - 0.097416*I -0.500882 + 0.002251*I] [-0.057002 - 0.495035*I -0.35059 - 0.625323*I 0.255514 - 0.415284*I] [ 0.394105 - 0.421778*I -0.392266 - 0.039345*I -0.352905 + 0.62195*I] sage: M.round(6) # random [ -1.528503 0.0 0.0] [ 0.459974 - 0.40061*I -1.741233 0.0] [-0.934304 + 0.148868*I 0.54833 + 0.073202*I -0.550725] sage: (A - M*G).zero_at(10^-12) [0.0 0.0 0.0] [0.0 0.0 0.0] [0.0 0.0 0.0] sage: (G*G.conjugate_transpose()) # random [0.9999999999999999 0.0 0.0] [ 0.0 0.9999999999999997 0.0] [ 0.0 0.0 1.0]
A rectangular matrix. Note that the
orthonormal
keyword is ignored in these cases.sage: A = matrix(RDF, [[-0.978325, -0.751994, 0.925305, -0.200512, 0.420458], ....: [-0.474877, -0.983403, 0.089836, 0.132218, 0.672965]]) sage: G, M = A.gram_schmidt(orthonormal=False) sage: G.round(6).zero_at(10^-6) [-0.607223 -0.466745 0.574315 -0.124453 0.260968] [ 0.123203 -0.617909 -0.530578 0.289773 0.487368] sage: M.round(6).zero_at(10^-6) [1.611147 0.0] [0.958116 0.867778] sage: (A-M*G).zero_at(10^-12) [0.0 0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0 0.0] sage: (G*G.transpose()).round(6).zero_at(10^-6) [1.0 0.0] [0.0 1.0]
Even though a set of vectors may be linearly dependent, no effort is made to decide when a zero vector is really the result of a relation of linear dependence. So in this regard, input matrices are treated as being of full rank. Try one of the base rings that provide exact results if you need exact results.
sage: entries = [[1,1,2], [2,1,3], [3,1,4]] sage: A = matrix(QQ, entries) sage: A.rank() 2 sage: B = matrix(RDF, entries) sage: G, M = B.gram_schmidt() sage: G.round(6) # random signs [-0.408248 -0.408248 -0.816497] [ 0.707107 -0.707107 -0.0] [ -0.57735 -0.57735 0.57735] sage: M.round(10) # random [-2.4494897428 0.0 0.0] [-3.6742346142 0.7071067812 0.0] [-4.8989794856 1.4142135624 0.0] sage: (A - M*G).zero_at(1e-14) [0.0 0.0 0.0] [0.0 0.0 0.0] [0.0 0.0 0.0] sage: (G*G.transpose()) # abs tol 1e-14 [0.9999999999999997 0.0 0.0] [ 0.0 0.9999999999999998 0.0] [ 0.0 0.0 1.0]
Exact Rings, Orthonormalization:
To scale a vector to unit length requires taking a square root, which often takes us outside the base ring. For the integers and the rationals, the field of algebraic numbers (
QQbar
) is big enough to contain what we need, but the price is that the computations are very slow, hence mostly of value for small cases or instruction. Now we need to use theorthonormal
keyword.sage: A = matrix(QQbar, [[6, -8, 1], ....: [4, 1, 3], ....: [6, 3, 3], ....: [7, 1, -5], ....: [7, -3, 5]]) sage: G, M = A.gram_schmidt(orthonormal=True) sage: G [ 0.5970223141259934? -0.7960297521679913? 0.09950371902099891?] [ 0.6063218341690895? 0.5289635311888953? 0.5937772444966257?] [ 0.5252981913594170? 0.2941669871612735? -0.798453250866314?] sage: M [ 10.04987562112089? 0 0] [ 1.890570661398980? 4.735582601355131? 0] [ 1.492555785314984? 7.006153332071100? 1.638930357041381?] [ 2.885607851608969? 1.804330147889395? 7.963520581008761?] [ 7.064764050490923? 5.626248468100069? -1.197679876299471?] sage: M*G-A [0 0 0] [0 0 0] [0 0 0] [0 0 0] [0 0 0] sage: (G*G.transpose()-identity_matrix(3)).norm() < 10^-10 True sage: G.row_space() == A.row_space() True
After github issue #14047, the matrix can also be over the algebraic reals
AA
:sage: A = matrix(AA, [[6, -8, 1], ....: [4, 1, 3], ....: [6, 3, 3], ....: [7, 1, -5], ....: [7, -3, 5]]) sage: G, M = A.gram_schmidt(orthonormal=True) sage: G [ 0.5970223141259934? -0.7960297521679913? 0.09950371902099891?] [ 0.6063218341690895? 0.5289635311888953? 0.5937772444966257?] [ 0.5252981913594170? 0.2941669871612735? -0.798453250866314?] sage: M [ 10.04987562112089? 0 0] [ 1.890570661398980? 4.735582601355131? 0] [ 1.492555785314984? 7.006153332071100? 1.638930357041381?] [ 2.885607851608969? 1.804330147889395? 7.963520581008761?] [ 7.064764050490923? 5.626248468100069? -1.197679876299471?]
Starting with complex numbers with rational real and imaginary parts. Note the use of the conjugate-transpose when checking the orthonormality.
sage: A = matrix(QQbar, [[ -2, -I - 1, 4*I + 2, -1], ....: [-4*I, -2*I + 17, 0, 9*I + 1], ....: [ 1, -2*I - 6, -I + 11, -5*I + 1]]) sage: G, M = A.gram_schmidt(orthonormal=True) sage: (M*G-A).norm() < 10^-10 True sage: id3 = G*G.conjugate().transpose() sage: (id3 - identity_matrix(3)).norm() < 10^-10 True sage: G.row_space() == A.row_space() # long time True
A square matrix with small rank. The zero vectors produced as a result of linear dependence get eliminated, so the rows of
G
are a basis for the row space ofA
.sage: A = matrix(QQbar, [[2, -6, 3, 8], ....: [1, -3, 2, 5], ....: [0, 0, 2, 4], ....: [2, -6, 3, 8]]) sage: A.change_ring(QQ).rank() 2 sage: G, M = A.gram_schmidt(orthonormal=True) sage: G [ 0.1881441736767195? -0.5644325210301583? 0.2822162605150792? 0.7525766947068779?] [-0.2502818123591464? 0.750845437077439? 0.3688363550555841? 0.4873908977520218?] sage: M [10.630145812734649? 0] [ 6.208757731331742? 0.6718090752798139?] [ 3.574739299857670? 2.687236301119256?] [10.630145812734649? 0] sage: M*G-A [0 0 0 0] [0 0 0 0] [0 0 0 0] [0 0 0 0] sage: (G*G.transpose()-identity_matrix(2)).norm() < 10^-10 True sage: G.row_space() == A.row_space() True
Exact Rings, Orthogonalization:
If we forego scaling orthogonal vectors to unit vectors, we can apply Gram-Schmidt to a much greater variety of rings. Use the
orthonormal=False
keyword (or assume it as the default). Note that now the orthogonality check creates a diagonal matrix whose diagonal entries are the squares of the lengths of the vectors.First, in the rationals, without involving
QQbar
.sage: A = matrix(QQ, [[-1, 3, 2, 2], ....: [-1, 0, -1, 0], ....: [-1, -2, -3, -1], ....: [ 1, 1, 2, 0]]) sage: A.rank() 3 sage: G, M = A.gram_schmidt() sage: G [ -1 3 2 2] [-19/18 1/6 -8/9 1/9] [ 2/35 -4/35 -2/35 9/35] sage: M [ 1 0 0] [ -1/18 1 0] [-13/18 59/35 1] [ 1/3 -48/35 -2] sage: M*G-A [0 0 0 0] [0 0 0 0] [0 0 0 0] [0 0 0 0] sage: G*G.transpose() [ 18 0 0] [ 0 35/18 0] [ 0 0 3/35] sage: G.row_space() == A.row_space() True
A complex subfield of the complex numbers.
sage: C.<z> = CyclotomicField(5) sage: A = matrix(C, [[ -z^3 - 2*z, -z^3 - 1, 2*z^3 - 2*z^2 + 2*z, 1], ....: [ z^3 - 2*z^2 + 1, -z^3 + 2*z^2 - z - 1, -1, z^2 + z], ....: [-1/2*z^3 - 2*z^2 + z + 1, -z^3 + z - 2, -2*z^3 + 1/2*z^2, 2*z^2 - z + 2]]) sage: G, M = A.gram_schmidt(orthonormal=False) sage: G [ -z^3 - 2*z -z^3 - 1 2*z^3 - 2*z^2 + 2*z 1] [ 155/139*z^3 - 161/139*z^2 + 31/139*z + 13/139 -175/139*z^3 + 180/139*z^2 - 125/139*z - 142/139 230/139*z^3 + 124/139*z^2 + 6/139*z + 19/139 -14/139*z^3 + 92/139*z^2 - 6/139*z - 95/139] [-10359/19841*z^3 - 36739/39682*z^2 + 24961/39682*z - 11879/39682 -28209/39682*z^3 - 3671/19841*z^2 + 51549/39682*z - 38613/39682 -42769/39682*z^3 - 615/39682*z^2 - 1252/19841*z - 14392/19841 4895/19841*z^3 + 57885/39682*z^2 - 46094/19841*z + 65747/39682] sage: M [ 1 0 0] [ 14/139*z^3 + 47/139*z^2 + 145/139*z + 95/139 1 0] [ -7/278*z^3 + 199/278*z^2 + 183/139*z + 175/278 -3785/39682*z^3 + 3346/19841*z^2 - 3990/19841*z + 2039/19841 1] sage: M*G - A [0 0 0 0] [0 0 0 0] [0 0 0 0] sage: G*G.conjugate().transpose() [ 15*z^3 + 15*z^2 + 28 0 0] [ 0 463/139*z^3 + 463/139*z^2 + 1971/139 0] [ 0 0 230983/19841*z^3 + 230983/19841*z^2 + 1003433/39682] sage: G.row_space() == A.row_space() True
A slightly edited legacy example.
sage: A = matrix(ZZ, 3, [-1, 2, 5, -11, 1, 1, 1, -1, -3]); A [ -1 2 5] [-11 1 1] [ 1 -1 -3] sage: G, mu = A.gram_schmidt() sage: G [ -1 2 5] [ -52/5 -1/5 -2] [ 2/187 36/187 -14/187] sage: mu [ 1 0 0] [ 3/5 1 0] [ -3/5 -7/187 1] sage: G.row(0) * G.row(1) 0 sage: G.row(0) * G.row(2) 0 sage: G.row(1) * G.row(2) 0
The relation between mu and A is as follows.
sage: mu*G == A True
- hadamard_bound()#
Return an int n such that the absolute value of the determinant of this matrix is at most \(10^n\).
This is got using both the row norms and the column norms.
This function only makes sense when the base field can be coerced to the real double field RDF or the MPFR Real Field with 53-bits precision.
EXAMPLES:
sage: a = matrix(ZZ, 3, [1,2,5,7,-3,4,2,1,123]) sage: a.hadamard_bound() 4 sage: a.det() -2014 sage: 10^4 10000
In this example the Hadamard bound has to be computed (automatically) using MPFR instead of doubles, since doubles overflow:
sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292]) sage: a.hadamard_bound() 12215 sage: len(str(a.det())) 12215
- hermite_form(include_zero_rows=True, transformation=False)#
Return the Hermite form of self, if it is defined.
INPUT:
include_zero_rows
– bool (default: True); if False the zero rows in the output matrix are deleted.transformation
– bool (default: False) a matrix U such that U*self == H.
OUTPUT:
matrix H
(optional) transformation matrix U such that U*self == H, possibly with zero rows deleted…
EXAMPLES:
sage: M = FunctionField(GF(7),'x').maximal_order() sage: K.<x> = FunctionField(GF(7)); M = K.maximal_order() sage: A = matrix(M, 2, 3, [x, 1, 2*x, x, 1+x, 2]) sage: A.hermite_form() [ x 1 2*x] [ 0 x 5*x + 2] sage: A.hermite_form(transformation=True) ( [ x 1 2*x] [1 0] [ 0 x 5*x + 2], [6 1] ) sage: A = matrix(M, 2, 3, [x, 1, 2*x, 2*x, 2, 4*x]) sage: A.hermite_form(transformation=True, include_zero_rows=False) ([ x 1 2*x], [1 0]) sage: H, U = A.hermite_form(transformation=True, include_zero_rows=True); H, U ( [ x 1 2*x] [1 0] [ 0 0 0], [5 1] ) sage: U*A == H True sage: H, U = A.hermite_form(transformation=True, include_zero_rows=False) sage: U*A [ x 1 2*x] sage: U*A == H True
- hessenberg_form()#
Return Hessenberg form of self.
If the base ring is merely an integral domain (and not a field), the Hessenberg form will (in general) only be defined over the fraction field of the base ring.
EXAMPLES:
sage: A = matrix(ZZ,4,[2, 1, 1, -2, 2, 2, -1, -1, -1,1,2,3,4,5,6,7]) sage: h = A.hessenberg_form(); h [ 2 -7/2 -19/5 -2] [ 2 1/2 -17/5 -1] [ 0 25/4 15/2 5/2] [ 0 0 58/5 3] sage: parent(h) Full MatrixSpace of 4 by 4 dense matrices over Rational Field sage: A.hessenbergize() Traceback (most recent call last): ... TypeError: Hessenbergize only possible for matrices over a field
- hessenbergize()#
Transform self to Hessenberg form.
The hessenberg form of a matrix \(A\) is a matrix that is similar to \(A\), so has the same characteristic polynomial as \(A\), and is upper triangular except possible for entries right below the diagonal.
ALGORITHM: See Henri Cohen’s first book.
EXAMPLES:
sage: A = matrix(QQ,3, [2, 1, 1, -2, 2, 2, -1, -1, -1]) sage: A.hessenbergize(); A [ 2 3/2 1] [ -2 3 2] [ 0 -3 -2]
sage: A = matrix(QQ,4, [2, 1, 1, -2, 2, 2, -1, -1, -1,1,2,3,4,5,6,7]) sage: A.hessenbergize(); A [ 2 -7/2 -19/5 -2] [ 2 1/2 -17/5 -1] [ 0 25/4 15/2 5/2] [ 0 0 58/5 3]
You can’t Hessenbergize an immutable matrix:
sage: A = matrix(QQ, 3, [1..9]) sage: A.set_immutable() sage: A.hessenbergize() 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).
- image()#
Return the image of the homomorphism on rows defined by right multiplication by this matrix: that is, the row-space.
EXAMPLES:
sage: MS1 = MatrixSpace(ZZ,4) sage: MS2 = MatrixSpace(QQ,6) sage: A = MS1.matrix([3,4,5,6,7,3,8,10,14,5,6,7,2,2,10,9]) sage: B = MS2.random_element()
sage: image(A) Free module of degree 4 and rank 4 over Integer Ring Echelon basis matrix: [ 1 0 0 426] [ 0 1 0 518] [ 0 0 1 293] [ 0 0 0 687]
sage: image(B) == B.row_module() True sage: image(B) == B.transpose().column_module() True
See also
- indefinite_factorization(algorithm='symmetric', check=True)#
Decomposes a symmetric or Hermitian matrix into a lower triangular matrix and a diagonal matrix.
INPUT:
self
- a square matrix over a ring. The base ring must have an implemented fraction field.algorithm
- default:'symmetric'
. Either'symmetric'
or'hermitian'
, according to if the input matrix is symmetric or hermitian.check
- default:True
- ifTrue
then performs the check that the matrix is consistent with thealgorithm
keyword.
OUTPUT:
A lower triangular matrix \(L\) with each diagonal element equal to \(1\) and a vector of entries that form a diagonal matrix \(D\). The vector of diagonal entries can be easily used to form the matrix, as demonstrated below in the examples.
For a symmetric matrix, \(A\), these will be related by
\[A = LDL^T\]If \(A\) is Hermitian matrix, then the transpose of \(L\) should be replaced by the conjugate-transpose of \(L\).
If any leading principal submatrix (a square submatrix in the upper-left corner) is singular then this method will fail with a
ValueError
.ALGORITHM:
The algorithm employed only uses field operations, but the computation of each diagonal entry has the potential for division by zero. The number of operations is of order \(n^3/3\), which is half the count for an LU decomposition. This makes it an appropriate candidate for solving systems with symmetric (or Hermitian) coefficient matrices.
See also
EXAMPLES:
There is no requirement that a matrix be positive definite, as indicated by the negative entries in the resulting diagonal matrix. The default is that the input matrix is symmetric.
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.is_symmetric() True sage: L, d = A.indefinite_factorization() sage: D = diagonal_matrix(d) sage: L [ 1 0 0 0 0] [-2 1 0 0 0] [ 3 -2 1 0 0] [ 2 -1 0 1 0] [-3 1 -3 1 1] sage: D [ 3 0 0 0 0] [ 0 -1 0 0 0] [ 0 0 5 0 0] [ 0 0 0 -2 0] [ 0 0 0 0 -1] sage: A == L*D*L.transpose() True
Optionally, Hermitian matrices can be factored and the result has a similar property (but not identical). Here, the field is all complex numbers with rational real and imaginary parts. As theory predicts, the diagonal entries will be real numbers.
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_hermitian() True sage: L, d = B.indefinite_factorization(algorithm='hermitian') sage: D = diagonal_matrix(d) sage: L [ 1 0 0] [ I + 2 1 0] [ -I + 1 2*I + 1 1] sage: D [ 2 0 0] [ 0 -2 0] [ 0 0 3] sage: B == L*D*L.conjugate_transpose() True
If a leading principal submatrix has zero determinant, this algorithm will fail. This will never happen with a positive definite matrix.
sage: A = matrix(QQ, [[21, 15, 12, -2], ....: [15, 12, 9, 6], ....: [12, 9, 7, 3], ....: [-2, 6, 3, 8]]) sage: A.is_symmetric() True sage: A[0:3,0:3].det() == 0 True sage: A.indefinite_factorization() Traceback (most recent call last): ... ValueError: 3x3 leading principal submatrix is singular, so cannot create indefinite factorization
This algorithm only depends on field operations, so outside of the singular submatrix situation, any matrix may be factored. This provides a reasonable alternative to the Cholesky decomposition.
sage: F.<a> = FiniteField(5^3) sage: A = matrix(F, ....: [[ a^2 + 2*a, 4*a^2 + 3*a + 4, 3*a^2 + a, 2*a^2 + 2*a + 1], ....: [4*a^2 + 3*a + 4, 4*a^2 + 2, 3*a, 2*a^2 + 4*a + 2], ....: [ 3*a^2 + a, 3*a, 3*a^2 + 2, 3*a^2 + 2*a + 3], ....: [2*a^2 + 2*a + 1, 2*a^2 + 4*a + 2, 3*a^2 + 2*a + 3, 3*a^2 + 2*a + 4]]) sage: A.is_symmetric() True sage: L, d = A.indefinite_factorization() sage: D = diagonal_matrix(d) sage: L [ 1 0 0 0] [4*a^2 + 4*a + 3 1 0 0] [ 3 4*a^2 + a + 2 1 0] [ 4*a^2 + 4 2*a^2 + 3*a + 3 2*a^2 + 3*a + 1 1] sage: D [ a^2 + 2*a 0 0 0] [ 0 2*a^2 + 2*a + 4 0 0] [ 0 0 3*a^2 + 4*a + 3 0] [ 0 0 0 a^2 + 3*a] sage: A == L*D*L.transpose() True
This works correctly for the 0x0 matrix:
sage: Matrix(0).indefinite_factorization() ([], ())
- integer_kernel(ring='ZZ')#
Return the kernel of this matrix over the given ring (which should be either the base ring, or a PID whose fraction field is the base ring).
Assume that the base field of this matrix has a numerator and denominator functions for its elements, e.g., it is the rational numbers or a fraction field. This function computes a basis over the integers for the kernel of self.
If the matrix is not coercible into QQ, then the PID itself should be given as a second argument, as in the third example below.
EXAMPLES:
sage: A = MatrixSpace(QQ, 4)(range(16)) sage: A.integer_kernel() Free module of degree 4 and rank 2 over Integer Ring Echelon basis matrix: [ 1 0 -3 2] [ 0 1 -2 1]
The integer kernel even makes sense for matrices with fractional entries:
sage: A = MatrixSpace(QQ, 2)([1/2, 0, 0, 0]) sage: A.integer_kernel() Free module of degree 2 and rank 1 over Integer Ring Echelon basis matrix: [0 1]
An example over a bigger ring:
sage: x = polygen(ZZ, 'x') sage: L.<w> = NumberField(x^2 - x + 2) # optional - sage.rings.number_field sage: OL = L.ring_of_integers() # optional - sage.rings.number_field sage: A = matrix(L, 2, [1, w/2]) # optional - sage.rings.number_field sage: A.integer_kernel(OL) # optional - sage.rings.number_field Free module of degree 2 and rank 1 over Maximal Order in Number Field in w with defining polynomial x^2 - x + 2 Echelon basis matrix: [ -1 -w + 1]
- inverse()#
Returns the inverse of self, without changing self.
Note that one can use the Python inverse operator to obtain the inverse as well.
See also
EXAMPLES:
sage: m = matrix([[1,2],[3,4]]) sage: m^(-1) [ -2 1] [ 3/2 -1/2] sage: m.inverse() [ -2 1] [ 3/2 -1/2] sage: ~m [ -2 1] [ 3/2 -1/2]
sage: m = matrix([[1,2],[3,4]], sparse=True) sage: m^(-1) [ -2 1] [ 3/2 -1/2] sage: m.inverse() [ -2 1] [ 3/2 -1/2] sage: ~m [ -2 1] [ 3/2 -1/2]
- inverse_positive_definite()#
Compute the inverse of a positive-definite matrix.
In accord with
is_positive_definite()
, only Hermitian matrices are considered positive-definite. Positive-definite matrices have several factorizations (Cholesky, LDLT, et cetera) that allow them to be inverted in a fast, numerically-stable way. This method uses an appropriate factorization, and is akin to thecholinv
andchol2inv
functions available in R, Octave, and Stata.You should ensure that your matrix is positive-definite before using this method. When in doubt, use the generic
inverse()
method instead.OUTPUT:
If the given matrix is positive-definite, the return value is the same as that of the
inverse()
method. If the matrix is not positive-definite, the behavior of this function is undefined.EXAMPLES:
A simple two-by-two matrix with rational entries:
sage: A = matrix(QQ, [[ 2, -1], ....: [-1, 2]]) sage: A.is_positive_definite() True sage: A.inverse_positive_definite() [2/3 1/3] [1/3 2/3] sage: A.inverse_positive_definite() == A.inverse() True
A matrix containing real roots:
sage: A = matrix(AA, [ [1, 0, sqrt(2)], ....: [0, sqrt(3), 0 ], ....: [sqrt(2), 0, sqrt(5)] ]) sage: A.is_positive_definite() True sage: B = matrix(AA, [ [2*sqrt(5) + 5, 0, -sqrt(8*sqrt(5) + 18)], ....: [0, sqrt(1/3), 0], ....: [-sqrt(8*sqrt(5) + 18), 0, sqrt(5) + 2] ]) sage: A.inverse_positive_definite() == B True sage: A*B == A.matrix_space().identity_matrix() True
A Hermitian (but not symmetric) matrix with complex entries:
sage: A = matrix(QQbar, [ [ 1, 0, I ], ....: [ 0, sqrt(5), 0 ], ....: [-I, 0, 3 ] ]) sage: A.is_positive_definite() True sage: B = matrix(QQbar, [ [ 3/2, 0, -I/2 ], ....: [ 0, sqrt(1/5), 0 ], ....: [ I/2, 0, 1/2 ] ]) sage: A.inverse_positive_definite() == B True sage: A*B == A.matrix_space().identity_matrix() True
- is_Z_operator_on(K)#
Determine if this matrix is a Z-operator on a cone.
We say that a matrix \(L\) is a Z-operator on a closed convex cone \(K\) if the inner product of \(Lx\) and \(s\) is nonpositive for all pairs of orthogonal vectors \(x\) in \(K\) and \(s\) in the dual of \(K\). This property need only be checked for generators of \(K\) and its dual.
A matrix is a Z-operator on \(K\) if and only if its negation is a cross-positive operator on \(K\).
To reliably check whether or not this matrix is a Z operator, its base ring must be either exact (for example, the rationals) or the symbolic ring. An exact ring is more reliable, but in some cases a matrix whose entries contain symbolic constants like \(e\) and \(\pi\) will work.
INPUT:
K
– a polyhedral closed convex cone.
OUTPUT:
If the base ring of this matrix is exact, then
True
will be returned if and only if this matrix is a Z-operator onK
.If the base ring of this matrix is symbolic, then the situation is more complicated:
True
will be returned if it can be proven that this matrix is a Z-operator onK
.False
will be returned if it can be proven that this matrix is not a Z-operator onK
.False
will also be returned if we can’t decide; specifically if we arrive at a symbolic inequality that cannot be resolved.
REFERENCES:
A. Berman and R. J. Plemmons. Nonnegative Matrices in the Mathematical Sciences. SIAM, Philadelphia, 1994.
EXAMPLES:
Z-matrices are Z-operators on the nonnegative orthant:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [-1, -2, 0], ....: [ 0, 2, -7], ....: [-3, 0, 3] ]) sage: L.is_Z_operator_on(K) True
Symbolic entries also work in some easy cases:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [-1, -e, 0 ], ....: [ 0, 2, -pi], ....: [-sqrt(2), 0, 3 ] ]) sage: L.is_Z_operator_on(K) True
- is_bistochastic(normalized=True)#
Returns
True
if this matrix is bistochastic.A matrix is said to be bistochastic if both the sums of the entries of each row and the sum of the entries of each column are equal to 1 and all entries are nonnegative.
INPUT:
normalized
– if set toTrue
(default), checks that the sums are equal to 1. When set toFalse
, checks that the row sums and column sums are all equal to some constant possibly different from 1.
EXAMPLES:
The identity matrix is clearly bistochastic:
sage: Matrix(5,5,1).is_bistochastic() True
The same matrix, multiplied by 2, is not bistochastic anymore, though is verifies the constraints of
normalized == False
:sage: (2 * Matrix(5,5,1)).is_bistochastic() False sage: (2 * Matrix(5,5,1)).is_bistochastic(normalized = False) True
Here is a matrix whose row and column sums is 1, but not all entries are nonnegative:
sage: m = matrix([[-1,2],[2,-1]]) sage: m.is_bistochastic() False
- is_cross_positive_on(K)#
Determine if this matrix is cross-positive on a cone.
We say that a matrix \(L\) is cross-positive on a closed convex cone \(K\) if the inner product of \(Lx\) and \(s\) is nonnegative for all pairs of orthogonal vectors \(x\) in \(K\) and \(s\) in the dual of \(K\). This property need only be checked for generators of \(K\) and its dual.
To reliably check whether or not this matrix is cross-positive, its base ring must be either exact (for example, the rationals) or the symbolic ring. An exact ring is more reliable, but in some cases a matrix whose entries contain symbolic constants like \(e\) and \(\pi\) will work.
INPUT:
K
– a polyhedral closed convex cone.
OUTPUT:
If the base ring of this matrix is exact, then
True
will be returned if and only if this matrix is cross-positive onK
.If the base ring of this matrix is symbolic, then the situation is more complicated:
True
will be returned if it can be proven that this matrix is cross-positive onK
.False
will be returned if it can be proven that this matrix is not cross-positive onK
.False
will also be returned if we can’t decide; specifically if we arrive at a symbolic inequality that cannot be resolved.
REFERENCES:
H. Schneider and M. Vidyasagar. Cross-positive matrices. SIAM Journal on Numerical Analysis, 7:508-519, 1970.
EXAMPLES:
Negative Z-matrices are cross-positive operators on the nonnegative orthant:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [-1, 2, 0], ....: [ 0, 2, 7], ....: [ 3, 0, 3] ]) sage: L.is_cross_positive_on(K) True
Symbolic entries also work in some easy cases:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [-1, e, 0 ], ....: [ 0, 2, pi], ....: [ sqrt(2), 0, 3 ] ]) sage: L.is_cross_positive_on(K) True
- is_diagonal()#
Return
True
if this matrix is a diagonal matrix.OUTPUT:
boolean
EXAMPLES:
sage: m = matrix(QQ,2,2,range(4)) sage: m.is_diagonal() False sage: m = matrix(QQ,2,[5,0,0,5]) sage: m.is_diagonal() True sage: m = matrix(QQ,2,[1,0,0,1]) sage: m.is_diagonal() True sage: m = matrix(QQ,2,[1,1,1,1]) sage: m.is_diagonal() False
- is_diagonalizable(base_field=None)#
Determine if the matrix is similar to a diagonal matrix.
INPUT:
base_field
- a new field to use for entries of the matrix.
OUTPUT:
If
self
is the matrix \(A\), then it is diagonalizable if there is an invertible matrix \(P\) and a diagonal matrix \(D\) such that\[P^{-1}AP = D\]This routine returns
True
ifself
is diagonalizable. The diagonal entries of the matrix \(D\) are the eigenvalues of \(A\).A matrix not diagonalizable over the base field may become diagonalizable by extending the base field to contain all of the eigenvalues. Over the rationals, the field of algebraic numbers,
sage.rings.qqbar
is a good choice.To obtain the matrices \(D\) and \(P\), use the
diagonalization()
method.ALGORITHM:
For each eigenvalue, this routine checks that the algebraic multiplicity (number of occurrences as a root of the characteristic polynomial) is equal to the geometric multiplicity (dimension of the eigenspace), which is sufficient to ensure a basis of eigenvectors for the columns of \(P\).
EXAMPLES:
A matrix that is diagonalizable over the rationals:
sage: A = matrix(QQ, [[-7, 16, 12, 0, 6], ....: [-9, 15, 0, 12, -27], ....: [ 9, -8, 11, -12, 51], ....: [ 3, -4, 0, -1, 9], ....: [-1, 0, -4, 4, -12]]) sage: A.is_diagonalizable() True sage: A.diagonalization() ( [ 2 0 0 0 0] [ 1 1 0 1 0] [ 0 3 0 0 0] [ 1/2 0 1 0 1] [ 0 0 3 0 0] [ 1/6 1 -3/2 2/3 -14/9] [ 0 0 0 -1 0] [ -1/6 0 -1/4 0 -1/3] [ 0 0 0 0 -1], [ -1/6 -1/3 1/3 -1/3 4/9] )
This is a matrix not diagonalizable over the rationals, but you can still get its Jordan form.
sage: A = matrix(QQ, [[-3, -14, 2, -1, 15], ....: [4, 6, -2, 3, -8], ....: [-2, -14, 0, 0, 10], ....: [3, 13, -2, 0, -11], ....: [-1, 6, 1, -3, 1]]) sage: A.is_diagonalizable() False sage: A.jordan_form(subdivide=False) [-1 1 0 0 0] [ 0 -1 0 0 0] [ 0 0 2 1 0] [ 0 0 0 2 1] [ 0 0 0 0 2]
If any eigenvalue of a matrix is outside the base ring, then this routine raises an error. However, the ring can be extended to contain the eigenvalues.
sage: A = matrix(QQ, [[1, 0, 1, 1, -1], ....: [0, 1, 0, 4, 8], ....: [2, 1, 3, 5, 1], ....: [2, -1, 1, 0, -2], ....: [0, -1, -1, -5, -8]]) sage: [e in QQ for e in A.eigenvalues()] [False, False, False, False, False] sage: A.is_diagonalizable() False sage: A.diagonalization() Traceback (most recent call last): ... ValueError: not diagonalizable over Rational Field sage: [e in QQbar for e in A.eigenvalues()] [True, True, True, True, True] sage: A.is_diagonalizable(base_field=QQbar) True
Other exact fields may be employed, though it will not always be possible to extend their base fields to contain all the eigenvalues.
sage: F.<b> = FiniteField(5^2) sage: A = matrix(F, [[ 4, 3*b + 2, 3*b + 1, 3*b + 4], ....: [2*b + 1, 4*b, 0, 2], ....: [ 4*b, b + 2, 2*b + 3, 3], ....: [ 2*b, 3*b, 4*b + 4, 3*b + 3]]) sage: A.is_diagonalizable() False sage: A.jordan_form() [ 4 1| 0 0] [ 0 4| 0 0] [---------------+---------------] [ 0 0|2*b + 1 1] [ 0 0| 0 2*b + 1] sage: F.<c> = QuadraticField(-7) sage: A = matrix(F, [[ c + 3, 2*c - 2, -2*c + 2, c - 1], ....: [2*c + 10, 13*c + 15, -13*c - 17, 11*c + 31], ....: [2*c + 10, 14*c + 10, -14*c - 12, 12*c + 30], ....: [ 0, 2*c - 2, -2*c + 2, 2*c + 2]]) sage: A.is_diagonalizable() True sage: A.diagonalization() ( [ 4 0 0 0] [ 1 0 1 0] [ 0 -2 0 0] [ 4 1 0 1] [ 0 0 c + 3 0] [ 5 1 -2/9 10/9] [ 0 0 0 c + 3], [ 1 0 -4/9 2/9] )
A trivial matrix is diagonalizable, trivially.
sage: A = matrix(QQ, 0, 0) sage: A.is_diagonalizable() True
A matrix must be square to be diagonalizable.
sage: A = matrix(QQ, 3, 4) sage: A.is_diagonalizable() Traceback (most recent call last): ... TypeError: not a square matrix
The matrix must have entries from a field, and it must be an exact field.
sage: A = matrix(ZZ, 4, range(16)) sage: A.is_diagonalizable() Traceback (most recent call last): ... ValueError: matrix entries must be from a field sage: A = matrix(RDF, 4, range(16)) sage: A.is_diagonalizable() Traceback (most recent call last): ... ValueError: base field must be exact, but Real Double Field is not
- is_lyapunov_like_on(K)#
Determine if this matrix is Lyapunov-like on a cone.
We say that a matrix \(L\) is Lyapunov-like on a closed convex cone \(K\) if the inner product of \(Lx\) and \(s\) is zero for all pairs of orthogonal vectors \(x\) in \(K\) and \(s\) in the dual of \(K\). This property need only be checked for generators of \(K\) and its dual.
An operator is Lyapunov-like on \(K\) if and only if both the operator itself and its negation are cross-positive on \(K\).
To reliably check whether or not this matrix is Lyapunov-like, its base ring must be either exact (for example, the rationals) or the symbolic ring. An exact ring is more reliable, but in some cases a matrix whose entries contain symbolic constants like \(e\) and \(\pi\) will work.
INPUT:
K
– a polyhedral closed convex cone.
OUTPUT:
If the base ring of this matrix is exact, then
True
will be returned if and only if this matrix is Lyapunov-like onK
.If the base ring of this matrix is symbolic, then the situation is more complicated:
True
will be returned if it can be proven that this matrix is Lyapunov-like onK
.False
will be returned if it can be proven that this matrix is not Lyapunov-like onK
.False
will also be returned if we can’t decide; specifically if we arrive at a symbolic inequality that cannot be resolved.
REFERENCES:
EXAMPLES:
Diagonal matrices are Lyapunov-like operators on the nonnegative orthant:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = diagonal_matrix(random_vector(QQ,3)) sage: L.is_lyapunov_like_on(K) True
Symbolic entries also work in some easy cases:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [e, 0, 0 ], ....: [0, pi, 0 ], ....: [0, 0, sqrt(2)] ]) sage: L.is_lyapunov_like_on(K) True
- is_normal()#
Returns
True
if the matrix commutes with its conjugate-transpose.OUTPUT:
True
if the matrix is square and commutes with its conjugate-transpose, andFalse
otherwise.Normal matrices are precisely those that can be diagonalized by a unitary matrix.
This routine is for matrices over exact rings and so may not work properly for matrices over
RR
orCC
. For matrices with approximate entries, the rings of double-precision floating-point numbers,RDF
andCDF
, are a better choice since thesage.matrix.matrix_double_dense.Matrix_double_dense.is_normal()
method has a tolerance parameter. This provides control over allowing for minor discrepancies between entries when checking equality.The result is cached.
EXAMPLES:
Hermitian matrices are normal.
sage: A = matrix(QQ, 5, 5, range(25)) + I*matrix(QQ, 5, 5, range(0, 50, 2)) sage: B = A*A.conjugate_transpose() sage: B.is_hermitian() True sage: B.is_normal() True
Circulant matrices are normal.
sage: G = graphs.CirculantGraph(20, [3, 7]) sage: D = digraphs.Circuit(20) sage: A = 3*D.adjacency_matrix() - 5*G.adjacency_matrix() sage: A.is_normal() True
Skew-symmetric matrices are normal.
sage: A = matrix(QQ, 5, 5, range(25)) sage: B = A - A.transpose() sage: B.is_skew_symmetric() True sage: B.is_normal() True
A small matrix that does not fit into any of the usual categories of normal matrices.
sage: A = matrix(ZZ, [[1, -1], ....: [1, 1]]) sage: A.is_normal() True sage: not A.is_hermitian() and not A.is_skew_symmetric() True
Sage has several fields besides the entire complex numbers where conjugation is non-trivial.
sage: F.<b> = QuadraticField(-7) sage: C = matrix(F, [[-2*b - 3, 7*b - 6, -b + 3], ....: [-2*b - 3, -3*b + 2, -2*b], ....: [ b + 1, 0, -2]]) sage: C = C*C.conjugate_transpose() sage: C.is_normal() True
A matrix that is nearly normal, but for a non-real diagonal entry.
sage: A = matrix(QQbar, [[ 2, 2-I, 1+4*I], ....: [ 2+I, 3+I, 2-6*I], ....: [1-4*I, 2+6*I, 5]]) sage: A.is_normal() False sage: A[1,1] = 132 sage: A.is_normal() True
Rectangular matrices are never normal.
sage: A = matrix(QQbar, 3, 4) sage: A.is_normal() False
A square, empty matrix is trivially normal.
sage: A = matrix(QQ, 0, 0) sage: A.is_normal() True
- is_one()#
Return True if this matrix is the identity matrix.
EXAMPLES:
sage: m = matrix(QQ,2,2,range(4)) sage: m.is_one() False sage: m = matrix(QQ,2,[5,0,0,5]) sage: m.is_one() False sage: m = matrix(QQ,2,[1,0,0,1]) sage: m.is_one() True sage: m = matrix(QQ,2,[1,1,1,1]) sage: m.is_one() False
- is_permutation_of(N, check=False)#
Return
True
if there exists a permutation of rows and columns sendingself
toN
andFalse
otherwise.INPUT:
N
– a matrix.check
– boolean (default:False
). IfFalse
return Boolean indicating whether there exists a permutation of rows and columns sending
self
toN
and False otherwise. IfTrue
return a tuple of a Boolean and a permutation mappingself
toN
if such a permutation exists, and (False
,None
) if it does not.
OUTPUT:
A Boolean or a tuple of a Boolean and a permutation.
EXAMPLES:
sage: M = matrix(ZZ,[[1,2,3],[3,5,3],[2,6,4]]) sage: M [1 2 3] [3 5 3] [2 6 4] sage: N = matrix(ZZ,[[1,2,3],[2,6,4],[3,5,3]]) sage: N [1 2 3] [2 6 4] [3 5 3] sage: M.is_permutation_of(N) True
Some examples that are not permutations of each other:
sage: N = matrix(ZZ,[[1,2,3],[4,5,6],[7,8,9]]) sage: N [1 2 3] [4 5 6] [7 8 9] sage: M.is_permutation_of(N) False sage: N = matrix(ZZ,[[1,2],[3,4]]) sage: N [1 2] [3 4] sage: M.is_permutation_of(N) False
And for when
check
is True:sage: N = matrix(ZZ,[[3,5,3],[2,6,4],[1,2,3]]) sage: N [3 5 3] [2 6 4] [1 2 3] sage: r = M.is_permutation_of(N, check=True) sage: r (True, ((1,2,3), ())) sage: p = r[1] sage: M.with_permuted_rows_and_columns(*p) == N True
- is_positive_definite(certificate=False)#
Determine if a matrix is positive-definite.
A matrix \(A\) is positive definite if it
is_hermitian()
and if, for every non-zero vector \(x\),\[\left\langle Ax, x \right\rangle > 0.\]ALGORITHM:
A Hermitian matrix is positive-definite if and only if the diagonal blocks in its
block_ldlt()
factorization are all 1-by-1 and have positive entries. We first check that the matrixis_hermitian()
, and then compute this factorization.INPUT:
self
– a matrixcertificate
– (default:False
) return the lower-triangular and diagonal parts of theblock_ldlt()
factorization when the matrix is positive-definite. Deprecated.
OUTPUT:
This routine will return
True
if the matrix is Hermitian and meets the condition above for the quadratic form.The base ring for the elements of the matrix must
Have a fraction field implemented; and
Be a subring of the real numbers, complex numbers, or symbolic ring.
If
certificate
isTrue
, a triplet(b, L, d)
will be returned instead, withb
containing the result (true or false). If the matrix is positive-definite, thenL
andd
will contain the lower-triangular and diagonal parts of theblock_ldlt()
factorization, respectively. Or if the matrix is not positive-definite (that is, ifb
isFalse
), then bothL
andd
will beNone
.See also
EXAMPLES:
A real symmetric matrix that is positive-definite, as evidenced by the positive determinants of its leading principal submatrices:
sage: A = matrix(QQ, [[ 4, -2, 4, 2], ....: [-2, 10, -2, -7], ....: [ 4, -2, 8, 4], ....: [ 2, -7, 4, 7]]) sage: A.is_positive_definite() True sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)] [4, 36, 144, 144]
A real symmetric matrix that is not positive-definite and a vector
u
that makes the corresponding quadratic form negative: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.is_positive_definite() False sage: u = vector(QQ, [2, 2, 0, 1, 0]) sage: (A*u).inner_product(u) -3
Another real symmetric matrix that is not positive-definite and a vector
u
that makes the corresponding quadratic form zero:sage: A = matrix(QQ, [[21, 15, 12, -2], ....: [15, 12, 9, 6], ....: [12, 9, 7, 3], ....: [-2, 6, 3, 8]]) sage: A.is_positive_definite() False sage: u = vector(QQ, [1,1,-3,0]) sage: (A*u).inner_product(u) 0
A complex Hermitian matrix that is positive-definite, confirmed by the positive determinants of its leading principal submatrices:
sage: x = polygen(ZZ, 'x') sage: C.<I> = NumberField(x^2 + 1, embedding=CC(0,1)) # optional - sage.rings.number_field sage: A = matrix(C, [[ 23, 17*I + 3, 24*I + 25, 21*I], # optional - sage.rings.number_field ....: [ -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_positive_definite() True sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)] [23, 576, 359540, 2842600]
An Hermitian matrix that is not positive-definite and a vector
u
that makes the corresponding quadratic form negative: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: u = vector(C, [-5 + 10*I, 4 - 3*I, 0]) sage: (B*u).hermitian_inner_product(u) -50
A positive-definite matrix over an algebraically-closed field, confirmed by the positive determinants of its leading principal submatrices:
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_positive_definite() True sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)] [2, 2, 6]
- is_positive_operator_on(K1, K2=None)#
Determine if this matrix is a positive operator on a cone.
A matrix is a positive operator on a cone if the image of the cone under the matrix is itself a subset of the cone. That concept can be extended to two cones: a matrix is a positive operator on a pair of cones if the image of the first cone is contained in the second cone.
To reliably check whether or not this matrix is a positive operator, its base ring must be either exact (for example, the rationals) or the symbolic ring. An exact ring is more reliable, but in some cases a matrix whose entries contain symbolic constants like \(e\) and \(\pi\) will work. Performance is best for integer or rational matrices, for which we can check the “is a subset of the other cone” condition quickly.
INPUT:
K1
– a polyhedral closed convex cone.K2
– (default:K1
) the codomain cone; this matrix is a positive operator if the image ofK1
is a subset ofK2
.
OUTPUT:
If the base ring of this matrix is exact, then
True
will be returned if and only if this matrix is a positive operator.If the base ring of this matrix is symbolic, then the situation is more complicated:
True
will be returned if it can be proven that this matrix is a positive operator.False
will be returned if it can be proven that this matrix is not a positive operator.False
will also be returned if we can’t decide; specifically if we arrive at a symbolic inequality that cannot be resolved.
REFERENCES:
A. Berman and P. Gaiha. A generalization of irreducible monotonicity. Linear Algebra and its Applications, 5:29-38, 1972.
A. Berman and R. J. Plemmons. Nonnegative Matrices in the Mathematical Sciences. SIAM, Philadelphia, 1994.
EXAMPLES:
Nonnegative matrices are positive operators on the nonnegative orthant:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = random_matrix(QQ,3).apply_map(abs) sage: L.is_positive_operator_on(K) True
Symbolic entries also work in some easy cases:
sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: L = matrix(SR, [ [0, e, 0 ], ....: [0, 2, pi], ....: [sqrt(2), 0, 0 ] ]) sage: L.is_positive_operator_on(K) True
Your matrix can be over any exact ring, for example the ring of univariate polynomials with rational coefficients:
sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) sage: K.is_full_space() True sage: x = polygen(ZZ, 'x') sage: L = matrix(QQ[x], [[x,0],[0,1]]) sage: L.is_positive_operator_on(K) True
- is_positive_semidefinite()#
Returns whether or not this matrix is positive-semidefinite.
By SageMath convention, positive (semi)definite matrices must be either real symmetric or complex Hermitian.
ALGORITHM:
Bunch and Kaufman [BK1977] describe a fast, numerically-stable scheme for computing the “inertia” of a matrix by way Sylvester’s inertia theorem and a block-\(LDL^{T}\) factorization. We perform this factorization, and read off the signs of the eigenvalues from the resulting diagonal blocks.
REFERENCES:
See also
EXAMPLES:
A positive-definite matrix:
sage: A = matrix(QQ, [ [2,1], ....: [1,2] ] ) sage: A.eigenvalues() [3, 1] sage: A.is_positive_semidefinite() True
A positive-semidefinite (but not positive-definite) matrix:
sage: A = matrix(QQ, [ [1,1], ....: [1,1] ] ) sage: A.eigenvalues() [2, 0] sage: A.is_positive_semidefinite() True
And finally, an indefinite matrix:
sage: A = matrix(QQ, [ [0,1], ....: [1,0] ] ) sage: A.eigenvalues() [1, -1] sage: A.is_positive_semidefinite() False
A non-Hermitian matrix cannot be positive-semidefinite, regardless of its eigenvalues:
sage: A = matrix(QQ, [ [2,1], ....: [0,0] ]) sage: A.eigenvalues() [2, 0] sage: A.is_positive_semidefinite() False
Any of the preceding examples are valid over inexact rings and with complex numbers as well:
sage: A = matrix(CDF, [ [ 2, I], ....: [-I, 2] ] ) sage: A.is_positive_semidefinite() True sage: A = matrix(CDF, [ [ 1, I], ....: [-I, 1] ] ) sage: A.is_positive_semidefinite() True sage: A = matrix(CDF, [ [0,I], ....: [I,0] ] ) sage: A.is_positive_semidefinite() False sage: A = matrix(CDF, [ [2,I], ....: [0,0] ]) sage: A.is_positive_semidefinite() False
- is_scalar(a=None)#
Return True if this matrix is a scalar matrix.
INPUT:
base_ring element a, which is chosen as self[0][0] if a = None
OUTPUT:
whether self is a scalar matrix (in fact the scalar matrix aI if a is input)
EXAMPLES:
sage: m = matrix(QQ,2,2,range(4)) sage: m.is_scalar(5) False sage: m = matrix(QQ,2,[5,0,0,5]) sage: m.is_scalar(5) True sage: m = matrix(QQ,2,[1,0,0,1]) sage: m.is_scalar(1) True sage: m = matrix(QQ,2,[1,1,1,1]) sage: m.is_scalar(1) False
- is_similar(other, transformation=False)#
Return
True
ifself
andother
are similar, i.e. related by a change-of-basis matrix.INPUT:
other
– a matrix, which should be square, and of the same size asself
.transformation
– default:False
- ifTrue
, the output may include the change-of-basis matrix (also known as the similarity transformation). See below for an exact description.
OUTPUT:
Two matrices, \(A\) and \(B\) are similar if they are square matrices of the same size and there is an invertible matrix \(S\) such that \(A=S^{-1}BS\). \(S\) can be interpreted as a change-of-basis matrix if \(A\) and \(B\) are viewed as matrix representations of the same linear transformation from a vector space to itself.
When
transformation=False
this method will returnTrue
if such a matrix \(S\) exists, otherwise it will returnFalse
. Whentransformation=True
the method returns a pair. The first part of the pair isTrue
orFalse
depending on if the matrices are similar. The second part of the pair is the change-of-basis matrix when the matrices are similar andNone
when the matrices are not similar.When a similarity transformation matrix
S
is requested, it will satisfyself = S.inverse()*other*S
.rings and coefficients
Inexact rings are not supported. Only rings having a fraction field can be used as coefficients.
The base rings for the matrices are promoted to a common field for the similarity check using rational form over this field.
If the fraction fields of both matrices are the same, this field is used. Otherwise, if the fraction fields are only related by a canonical coercion, the common coercion field is used.
In all cases, the result is about similarity over a common field.
similarity transformation
For computation of the similarity transformation, the matrices are first checked to be similar over their common field.
In this case, a similarity transformation is then searched for over the common field. If this fails, the matrices are promoted to the algebraic closure of their common field (whenever it is available) and a similarity transformation is looked for over the algebraic closure.
For example, matrices over the rationals may be promoted to the field of algebraic numbers (
QQbar
) for computation of the similarity transformation.Warning
When the two matrices are similar, this routine may fail to find the similarity transformation. A technical explanation follows.
The similarity check is accomplished with rational form, which will be successful for any pair of matrices over the same field. However, the computation of rational form does not provide a transformation. So we instead compute Jordan form, which does provide a transformation. But Jordan form will require that the eigenvalues of the matrix can be represented within Sage, requiring the existence of the appropriate extension field. When this is not possible, a
RuntimeError
is raised, as demonstrated in an example below.EXAMPLES:
The two matrices in this example were constructed to be similar. The computations happen in the field of algebraic numbers, but we are able to convert the change-of-basis matrix back to the rationals (which may not always be possible).
sage: A = matrix(ZZ, [[-5, 2, -11], ....: [-6, 7, -42], ....: [0, 1, -6]]) sage: B = matrix(ZZ, [[ 1, 12, 3], ....: [-1, -6, -1], ....: [ 0, 6, 1]]) sage: A.is_similar(B) True sage: _, T = A.is_similar(B, transformation=True) sage: T [ 1.00000000000000? + 0.?e-14*I 0.?e-14 + 0.?e-14*I 0.?e-14 + 0.?e-14*I] [-0.66666666666667? + 0.?e-15*I 0.166666666666667? + 0.?e-15*I -0.83333333333334? + 0.?e-14*I] [ 0.66666666666667? + 0.?e-14*I 0.?e-14 + 0.?e-14*I -0.33333333333333? + 0.?e-14*I] sage: T.change_ring(QQ) [ 1 0 0] [-2/3 1/6 -5/6] [ 2/3 0 -1/3] sage: A == T.inverse()*B*T True
Other exact fields are supported.
sage: F.<a> = FiniteField(7^2) sage: A = matrix(F,[[2*a + 5, 6*a + 6, a + 3], ....: [ a + 3, 2*a + 2, 4*a + 2], ....: [2*a + 6, 5*a + 5, 3*a]]) sage: B = matrix(F,[[5*a + 5, 6*a + 4, a + 1], ....: [ a + 5, 4*a + 3, 3*a + 3], ....: [3*a + 5, a + 4, 5*a + 6]]) sage: A.is_similar(B) True sage: B.is_similar(A) True sage: _, T = A.is_similar(B, transformation=True) sage: T [ 1 0 0] [6*a + 1 4*a + 3 4*a + 2] [6*a + 3 3*a + 5 3*a + 6] sage: A == T.inverse()*B*T True
Two matrices with different sets of eigenvalues, so they cannot possibly be similar.
sage: A = matrix(QQ, [[ 2, 3, -3, -6], ....: [ 0, 1, -2, -8], ....: [-3, -3, 4, 3], ....: [-1, -2, 2, 6]]) sage: B = matrix(QQ, [[ 1, 1, 2, 4], ....: [-1, 2, -3, -7], ....: [-2, 3, -4, -7], ....: [ 0, -1, 0, 0]]) sage: A.eigenvalues() == B.eigenvalues() False sage: A.is_similar(B, transformation=True) (False, None)
Similarity is an equivalence relation, so this routine computes a representative of the equivalence class for each matrix, the rational form, as provided by
rational_form()
. The matrices below have identical eigenvalues (as evidenced by equal characteristic polynomials), but slightly different rational forms, and hence are not similar.sage: A = matrix(QQ, [[ 19, -7, -29], ....: [-16, 11, 30], ....: [ 15, -7, -25]]) sage: B = matrix(QQ, [[-38, -63, 42], ....: [ 14, 25, -14], ....: [-14, -21, 18]]) sage: A.charpoly() == B.charpoly() True sage: A.rational_form() [ 0 0 -48] [ 1 0 8] [ 0 1 5] sage: B.rational_form() [ 4| 0 0] [--+-----] [ 0| 0 12] [ 0| 1 1] sage: A.is_similar(B) False
Obtaining the transformation between two similar matrices requires the Jordan form, which requires computing the eigenvalues of the matrix, which may not lie in the field used for entries of the matrix. In this unfortunate case, the computation of the transformation may fail with a
RuntimeError
, EVEN when the matrices are similar. This is not the case for matrices over the integers, rationals or algebraic numbers, since the computations are done in the algebraically closed field of algebraic numbers. Here is an example where the similarity is obvious by design, but we are not able to resurrect a similarity transformation.sage: F.<a> = FiniteField(7^2) sage: C = matrix(F,[[ a + 2, 5*a + 4], ....: [6*a + 6, 6*a + 4]]) sage: S = matrix(F, [[0, 1], ....: [1, 0]]) sage: D = S.inverse()*C*S sage: C.is_similar(D) True sage: C.is_similar(D, transformation=True) Traceback (most recent call last): ... RuntimeError: unable to compute transformation for similar matrices sage: C.jordan_form() Traceback (most recent call last): ... RuntimeError: Some eigenvalue does not exist in Finite Field in a of size 7^2.
An example over a finite field of prime order, which uses the algebraic closure of this field to find the change-of-basis matrix:
sage: cox = posets.TamariLattice(3).coxeter_transformation() sage: M = cox.change_ring(GF(3)) sage: M.is_similar(M**3, True) # long time ( [1 0 0 0 0] [0 1 1 0 2] [0 0 0 0 1] [1 2 0 2 1] True, [0 0 1 0 0] )
Inexact rings and fields are not supported.
sage: A = matrix(CDF, 2, 2, range(4)) sage: B = copy(A) sage: A.is_similar(B) Traceback (most recent call last): ... TypeError: matrix entries must come from an exact field, not Complex Double Field
Base rings for the matrices need to have a fraction field. So in particular, the ring needs to be at least an integral domain.
sage: Z6 = Integers(6) sage: A = matrix(Z6, 2, 2, range(4)) sage: A.is_similar(A) Traceback (most recent call last): ... ValueError: base ring of a matrix needs a fraction field, maybe the ring is not an integral domain
If the fraction fields of the entries are unequal and do not coerce in a common field, it is an error.
sage: A = matrix(GF(3), 2, 2, range(4)) sage: B = matrix(GF(2), 2, 2, range(4)) sage: A.is_similar(B, transformation=True) Traceback (most recent call last): ... TypeError: no common canonical parent for objects with parents: 'Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 3' and 'Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 2'
A matrix over the integers and a matrix over the algebraic numbers will be compared over the algebraic numbers (by coercion of
QQ
inQQbar
).sage: A = matrix(ZZ, 2, 2, range(4)) sage: B = matrix(QQbar, 2, 2, range(4)) sage: A.is_similar(B) True
- is_triangular(side='lower')#
Return
True
if this matrix is a triangular matrix.INPUT:
side
– either"lower"
(default) or"upper"
OUTPUT:
boolean
EXAMPLES:
sage: m = matrix(QQ, 2, 2, range(4)) sage: m.is_triangular() False sage: m = matrix(QQ, 2, [5, 0, 0, 5]) sage: m.is_triangular() True sage: m = matrix(QQ, 2, [1, 2, 0, 1]) sage: m.is_triangular("upper") True sage: m.is_triangular("lower") False
- is_unitary()#
Return
True
if the columns of the matrix are an orthonormal basis.For a matrix with real entries this determines if a matrix is “orthogonal” and for a matrix with complex entries this determines if the matrix is “unitary.”
OUTPUT:
True
if the matrix is square and its conjugate-transpose is its inverse, andFalse
otherwise. In other words, a matrix is orthogonal or unitary if the product of its conjugate-transpose times the matrix is the identity matrix.For numerical matrices a specialized routine available over
RDF
andCDF
is a good choice.EXAMPLES:
sage: A = matrix(QQbar, [[(1/sqrt(5))*(1+i), (1/sqrt(55))*(3+2*I), (1/sqrt(22))*(2+2*I)], ....: [(1/sqrt(5))*(1-i), (1/sqrt(55))*(2+2*I), (1/sqrt(22))*(-3+I)], ....: [ (1/sqrt(5))*I, (1/sqrt(55))*(3-5*I), (1/sqrt(22))*(-2)]]) sage: A.is_unitary() True
A permutation matrix is always orthogonal.
sage: sigma = Permutation([1,3,4,5,2]) sage: P = sigma.to_matrix(); P [1 0 0 0 0] [0 0 0 0 1] [0 1 0 0 0] [0 0 1 0 0] [0 0 0 1 0] sage: P.is_unitary() True sage: P.change_ring(GF(3)).is_unitary() True sage: P.change_ring(GF(3)).is_unitary() True
A square matrix far from unitary.
sage: A = matrix(QQ, 4, range(16)) sage: A.is_unitary() False
Rectangular matrices are never unitary.
sage: A = matrix(QQbar, 3, 4) sage: A.is_unitary() False
- jordan_form(base_ring=None, sparse=False, subdivide=True, transformation=False, eigenvalues=None, check_input=True)#
Compute the Jordan normal form of this square matrix \(A\), if it exists.
This computation is performed in a naive way using the ranks of powers of \(A-xI\), where \(x\) is an eigenvalue of the matrix \(A\). If desired, a transformation matrix \(P\) can be returned, which is such that the Jordan canonical form is given by \(P^{-1} A P\); if the matrix is diagonalizable, this equals to eigendecomposition or spectral decomposition.
INPUT:
base_ring
- Ring in which to compute the Jordan form.sparse
- (defaultFalse
) Ifsparse=True
, return a sparse matrix.subdivide
- (defaultTrue
) Ifsubdivide=True
, the subdivisions for the Jordan blocks in the matrix are shown.transformation
- (defaultFalse
) Iftransformation=True
, computes also the transformation matrix.eigenvalues
- (defaultNone
) A complete set of roots, with multiplicity, of the characteristic polynomial of \(A\), encoded as a list of pairs, each having the form \((r, m)\) with \(r\) a root and \(m\) its multiplicity. If this isNone
, then Sage computes this list itself, but this is only possible over base rings in whose quotient fields polynomial factorization is implemented. Over all other rings, providing this list manually is the only way to compute Jordan normal forms.check_input
- (defaultTrue
) A Boolean specifying whether the listeigenvalues
(if provided) has to be checked for correctness. Set this toFalse
for a speedup if the eigenvalues are known to be correct.
Note
Currently, the Jordan normal form is not computed over inexact rings in any but the trivial cases when the matrix is either \(0 \times 0\) or \(1 \times 1\).
In the case of exact rings, this method does not compute any generalized form of the Jordan normal form, but is only able to compute the result if the characteristic polynomial of the matrix splits over the specific base ring.
Note that the base ring must be a field or a ring with an implemented fraction field.
EXAMPLES:
sage: a = matrix(ZZ,4,[1, 0, 0, 0, 0, 1, 0, 0, ....: 1, -1, 1, 0, 1, -1, 1, 2]); a [ 1 0 0 0] [ 0 1 0 0] [ 1 -1 1 0] [ 1 -1 1 2] sage: a.jordan_form() [2|0 0|0] [-+---+-] [0|1 1|0] [0|0 1|0] [-+---+-] [0|0 0|1] sage: a.jordan_form(subdivide=False) [2 0 0 0] [0 1 1 0] [0 0 1 0] [0 0 0 1] sage: b = matrix(ZZ,3,3,range(9)); b [0 1 2] [3 4 5] [6 7 8] sage: b.jordan_form() Traceback (most recent call last): ... RuntimeError: Some eigenvalue does not exist in Rational Field. sage: b.jordan_form(RealField(15)) Traceback (most recent call last): ... ValueError: Jordan normal form not implemented over inexact rings.
Here we need to specify a field, since the eigenvalues are not defined in the smallest ring containing the matrix entries (github issue #14508):
sage: c = matrix([[0,1,0],[0,0,1],[1,0,0]]) sage: c.jordan_form(CyclotomicField(3)) [ 1| 0| 0] [----------+----------+----------] [ 0| zeta3| 0] [----------+----------+----------] [ 0| 0|-zeta3 - 1]
If you need the transformation matrix as well as the Jordan form of
self
, then pass the optiontransformation=True
. For example:sage: m = matrix([[5,4,2,1],[0,1,-1,-1],[-1,-1,3,0],[1,1,-1,2]]); m [ 5 4 2 1] [ 0 1 -1 -1] [-1 -1 3 0] [ 1 1 -1 2] sage: jf, p = m.jordan_form(transformation=True) sage: jf [2|0|0 0] [-+-+---] [0|1|0 0] [-+-+---] [0|0|4 1] [0|0|0 4] sage: ~p * m * p [2 0 0 0] [0 1 0 0] [0 0 4 1] [0 0 0 4]
Note that for matrices over inexact rings, we do not attempt to compute the Jordan normal form, since it is not numerically stable:
sage: b = matrix(ZZ,3,3,range(9)) sage: jf, p = b.jordan_form(RealField(15), transformation=True) Traceback (most recent call last): ... ValueError: Jordan normal form not implemented over inexact rings.
- kernel(*args, **kwds)#
Returns the left kernel of this matrix, as a vector space or free module. This is the set of vectors
x
such thatx*self = 0
.Note
For the right kernel, use
right_kernel()
. The methodkernel()
is exactly equal toleft_kernel()
.For inexact rings use
right_kernel_matrix()
withbasis='computed'
(on the transpose of the matrix) to avoid echolonizing.INPUT:
algorithm
- default: ‘default’ - a keyword that selects the algorithm employed. Allowable values are:‘default’ - allows the algorithm to be chosen automatically
‘generic’ - naive algorithm usable for matrices over any field
‘flint’ - FLINT library code for matrices over the rationals or the integers
‘pari’ - PARI library code for matrices over number fields or the integers
‘padic’ - padic algorithm from IML library for matrices over the rationals and integers
‘pluq’ - PLUQ matrix factorization for matrices mod 2
basis
- default: ‘echelon’ - a keyword that describes the format of the basis used to construct the left kernel. Allowable values are:‘echelon’: the basis matrix is returned in echelon form
‘pivot’ : each basis vector is computed from the reduced row-echelon form of
self
by placing a single one in a non-pivot column and zeros in the remaining non-pivot columns. Only available for matrices over fields.‘LLL’: an LLL-reduced basis. Only available for matrices over the integers.
OUTPUT:
A vector space or free module whose degree equals the number of rows in
self
and which contains all the vectorsx
such thatx*self = 0
.If
self
has 0 rows, the kernel has dimension 0, while ifself
has 0 columns the kernel is the entire ambient vector space.The result is cached. Requesting the left kernel a second time, but with a different basis format, will return the cached result with the format from the first computation.
Note
For much more detailed documentation of the various options see
right_kernel()
, since this method just computes the right kernel of the transpose ofself
.EXAMPLES:
Over the rationals with a basis matrix in echelon form.
sage: A = matrix(QQ, [[1, 2, 4, -7, 4], ....: [1, 1, 0, 2, -1], ....: [1, 0, 3, -3, 1], ....: [0, -1, -1, 3, -2], ....: [0, 0, -1, 2, -1]]) sage: A.left_kernel() Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: [ 1 0 -1 2 -1] [ 0 1 -1 1 -4]
Over a finite field, with a basis matrix in “pivot” format.
sage: A = matrix(FiniteField(7), [[5, 0, 5, 2, 4], ....: [1, 3, 2, 3, 6], ....: [1, 1, 6, 5, 3], ....: [2, 5, 6, 0, 0]]) sage: A.kernel(basis='pivot') Vector space of degree 4 and dimension 2 over Finite Field of size 7 User basis matrix: [5 2 1 0] [6 3 0 1]
The left kernel of a zero matrix is the entire ambient vector space whose degree equals the number of rows of
self
(i.e. everything).sage: A = MatrixSpace(QQ, 3, 4)(0) sage: A.kernel() Vector space of degree 3 and dimension 3 over Rational Field Basis matrix: [1 0 0] [0 1 0] [0 0 1]
We test matrices with no rows or columns.
sage: A = matrix(QQ, 2, 0) sage: A.left_kernel() Vector space of degree 2 and dimension 2 over Rational Field Basis matrix: [1 0] [0 1] sage: A = matrix(QQ, 0, 2) sage: A.left_kernel() Vector space of degree 0 and dimension 0 over Rational Field Basis matrix: []
The results are cached. Note that requesting a new format for the basis is ignored and the cached copy is returned. Work with a copy if you need a new left kernel, or perhaps investigate the
right_kernel_matrix()
method on the transpose, which does not cache its results and is more flexible.sage: A = matrix(QQ, [[1,1],[2,2]]) sage: K1 = A.left_kernel() sage: K1 Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 -1/2] sage: K2 = A.left_kernel() sage: K1 is K2 True sage: K3 = A.left_kernel(basis='pivot') sage: K3 Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 -1/2] sage: B = copy(A) sage: K3 = B.left_kernel(basis='pivot') sage: K3 Vector space of degree 2 and dimension 1 over Rational Field User basis matrix: [-2 1] sage: K3 is K1 False sage: K3 == K1 True
- kernel_on(V, poly=None, check=True)#
Return the kernel of self restricted to the invariant subspace V. The result is a vector subspace of V, which is also a subspace of the ambient space.
INPUT:
V
- vector subspacecheck
- (optional) default: True; whether to check that V is invariant under the action of self.poly
- (optional) default: None; if not None, compute instead the kernel of poly(self) on V.
OUTPUT:
a subspace
Warning
This function does not check that V is in fact invariant under self if check is False. With check False this function is much faster.
EXAMPLES:
sage: t = matrix(QQ, 4, [39, -10, 0, -12, 0, 2, 0, -1, 0, 1, -2, 0, 0, 2, 0, -2]); t [ 39 -10 0 -12] [ 0 2 0 -1] [ 0 1 -2 0] [ 0 2 0 -2] sage: t.fcp() (x - 39) * (x + 2) * (x^2 - 2) sage: s = (t-39)*(t^2-2) sage: V = s.kernel(); V Vector space of degree 4 and dimension 3 over Rational Field Basis matrix: [1 0 0 0] [0 1 0 0] [0 0 0 1] sage: s.restrict(V) [0 0 0] [0 0 0] [0 0 0] sage: s.kernel_on(V) Vector space of degree 4 and dimension 3 over Rational Field Basis matrix: [1 0 0 0] [0 1 0 0] [0 0 0 1] sage: k = t-39 sage: k.restrict(V) [ 0 -10 -12] [ 0 -37 -1] [ 0 2 -41] sage: ker = k.kernel_on(V); ker Vector space of degree 4 and dimension 1 over Rational Field Basis matrix: [ 1 -2/7 0 -2/7] sage: ker.0 * k (0, 0, 0, 0)
Test that github issue #9425 is fixed.
sage: V = span([[1/7,0,0] ,[0,1,0]], ZZ); V Free module of degree 3 and rank 2 over Integer Ring Echelon basis matrix: [1/7 0 0] [ 0 1 0] sage: T = matrix(ZZ,3,[1,0,0,0,0,0,0,0,0]); T [1 0 0] [0 0 0] [0 0 0] sage: W = T.kernel_on(V); W.basis() [ (0, 1, 0) ] sage: W.is_submodule(V) True
- left_eigenmatrix(other=None)#
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
self
is a square matrix \(A\), then the output is a diagonal matrix \(D\) and a matrix \(P\) such that\[P A = D P,\]where the rows of \(P\) are eigenvectors of \(A\) and the diagonal entries of \(D\) are the corresponding eigenvalues.
If a matrix \(B\) is passed as optional argument, the output is a solution to the generalized eigenvalue problem such that
\[P A = D P B.\]The ordinary eigenvalue problem is equivalent to the generalized one if \(B\) is the identity matrix.
The generalized eigenvector decomposition is currently only implemented for matrices over
RDF
andCDF
.EXAMPLES:
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
Because \(P\) is invertible, \(A\) is diagonalizable.
sage: A == (~P)*D*P 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: 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
A generalized eigenvector decomposition:
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
The matrix \(B\) in a generalized eigenvalue problem may be singular:
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]
In this case, we can still verify the eigenvector equation for the first eigenvalue and first eigenvector:
sage: l = D[0, 0] sage: v = P[0, :] sage: (v * A - l * v * B).norm() < 1e-14 True
The second eigenvector is contained in the left kernel of \(B\):
sage: (P[1, :] * B).norm() < 1e-14 True
- left_eigenspaces(format='all', var='a', algebraic_multiplicity=False)#
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 classes.format
- default:None
'all'
- attempts to create every eigenspace. This will always be possible for matrices with rational entries.'galois'
- for each irreducible factor of the characteristic polynomial, a single eigenspace will be output for a single root/eigenvalue for the irreducible factor.None
- Uses the ‘all’ format if the base ring is contained in an algebraically closed field which is implemented. Otherwise, uses the ‘galois’ format.
var
- default: ‘a’ - variable name used to represent elements of the root field of each irreducible factor of the characteristic polynomial. If var=’a’, then the root fields will be in terms of a0, a1, a2, …., where the numbering runs across all the irreducible factors of the characteristic polynomial, even for linear factors.algebraic_multiplicity
- default: False - whether or not to include the algebraic multiplicity of each eigenvalue in the output. See the discussion below.
OUTPUT:
If algebraic_multiplicity=False, return a list of pairs (e, V) where e is an eigenvalue of the matrix, and V is the corresponding left eigenspace. For Galois conjugates of eigenvalues, there may be just one representative eigenspace, depending on the
format
keyword.If algebraic_multiplicity=True, return a list of triples (e, V, n) where e and V are as above and n is the algebraic multiplicity of the eigenvalue.
Warning
Uses a somewhat naive algorithm (simply factors the characteristic polynomial and computes kernels directly over the extension field).
EXAMPLES:
We compute the left eigenspaces of a \(3\times 3\) rational matrix. First, we request \(all\) of the eigenvalues, so the results are in the field of algebraic numbers, \(QQbar\). Then we request just one eigenspace per irreducible factor of the characteristic polynomial with the \(galois\) keyword.
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 [ (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: 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
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') [ (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: 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) ]
Next we compute the left eigenspaces over the finite field of order 11.
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]) ]
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: 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: 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]
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: 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') 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') [ (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]) ]
- left_eigenvectors(other=None, extend=True)#
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 ofself
isRDF
orCDF
)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: 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
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() [(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) [(2, [ (0, 0, 1) ], 1)]
- left_kernel(*args, **kwds)#
Returns the left kernel of this matrix, as a vector space or free module. This is the set of vectors
x
such thatx*self = 0
.Note
For the right kernel, use
right_kernel()
. The methodkernel()
is exactly equal toleft_kernel()
.For inexact rings use
right_kernel_matrix()
withbasis='computed'
(on the transpose of the matrix) to avoid echolonizing.INPUT:
algorithm
- default: ‘default’ - a keyword that selects the algorithm employed. Allowable values are:‘default’ - allows the algorithm to be chosen automatically
‘generic’ - naive algorithm usable for matrices over any field
‘flint’ - FLINT library code for matrices over the rationals or the integers
‘pari’ - PARI library code for matrices over number fields or the integers
‘padic’ - padic algorithm from IML library for matrices over the rationals and integers
‘pluq’ - PLUQ matrix factorization for matrices mod 2
basis
- default: ‘echelon’ - a keyword that describes the format of the basis used to construct the left kernel. Allowable values are:‘echelon’: the basis matrix is returned in echelon form
‘pivot’ : each basis vector is computed from the reduced row-echelon form of
self
by placing a single one in a non-pivot column and zeros in the remaining non-pivot columns. Only available for matrices over fields.‘LLL’: an LLL-reduced basis. Only available for matrices over the integers.
OUTPUT:
A vector space or free module whose degree equals the number of rows in
self
and which contains all the vectorsx
such thatx*self = 0
.If
self
has 0 rows, the kernel has dimension 0, while ifself
has 0 columns the kernel is the entire ambient vector space.The result is cached. Requesting the left kernel a second time, but with a different basis format, will return the cached result with the format from the first computation.
Note
For much more detailed documentation of the various options see
right_kernel()
, since this method just computes the right kernel of the transpose ofself
.EXAMPLES:
Over the rationals with a basis matrix in echelon form.
sage: A = matrix(QQ, [[1, 2, 4, -7, 4], ....: [1, 1, 0, 2, -1], ....: [1, 0, 3, -3, 1], ....: [0, -1, -1, 3, -2], ....: [0, 0, -1, 2, -1]]) sage: A.left_kernel() Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: [ 1 0 -1 2 -1] [ 0 1 -1 1 -4]
Over a finite field, with a basis matrix in “pivot” format.
sage: A = matrix(FiniteField(7), [[5, 0, 5, 2, 4], ....: [1, 3, 2, 3, 6], ....: [1, 1, 6, 5, 3], ....: [2, 5, 6, 0, 0]]) sage: A.kernel(basis='pivot') Vector space of degree 4 and dimension 2 over Finite Field of size 7 User basis matrix: [5 2 1 0] [6 3 0 1]
The left kernel of a zero matrix is the entire ambient vector space whose degree equals the number of rows of
self
(i.e. everything).sage: A = MatrixSpace(QQ, 3, 4)(0) sage: A.kernel() Vector space of degree 3 and dimension 3 over Rational Field Basis matrix: [1 0 0] [0 1 0] [0 0 1]
We test matrices with no rows or columns.
sage: A = matrix(QQ, 2, 0) sage: A.left_kernel() Vector space of degree 2 and dimension 2 over Rational Field Basis matrix: [1 0] [0 1] sage: A = matrix(QQ, 0, 2) sage: A.left_kernel() Vector space of degree 0 and dimension 0 over Rational Field Basis matrix: []
The results are cached. Note that requesting a new format for the basis is ignored and the cached copy is returned. Work with a copy if you need a new left kernel, or perhaps investigate the
right_kernel_matrix()
method on the transpose, which does not cache its results and is more flexible.sage: A = matrix(QQ, [[1,1],[2,2]]) sage: K1 = A.left_kernel() sage: K1 Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 -1/2] sage: K2 = A.left_kernel() sage: K1 is K2 True sage: K3 = A.left_kernel(basis='pivot') sage: K3 Vector space of degree 2 and dimension 1 over Rational Field Basis matrix: [ 1 -1/2] sage: B = copy(A) sage: K3 = B.left_kernel(basis='pivot') sage: K3 Vector space of degree 2 and dimension 1 over Rational Field User basis matrix: [-2 1] sage: K3 is K1 False sage: K3 == K1 True
- left_kernel_matrix(*args, **kwds)#
Returns a matrix whose rows form a basis for the left kernel of
self
.This method is a thin wrapper around
right_kernel_matrix()
. For supported parameters and input/output formats, see there.EXAMPLES:
sage: M = matrix([[1,2],[3,4],[5,6]]) sage: K = M.left_kernel_matrix(); K [ 1 -2 1] sage: K * M [0 0]
- left_nullity()#
Return the (left) nullity of this matrix, which is the dimension of the (left) kernel of this matrix acting from the right on row vectors.
EXAMPLES:
sage: M = Matrix(QQ,[[1,0,0,1],[0,1,1,0],[1,1,1,0]]) sage: M.nullity() 0 sage: M.left_nullity() 0
sage: A = M.transpose() sage: A.nullity() 1 sage: A.left_nullity() 1
sage: M = M.change_ring(ZZ) sage: M.nullity() 0 sage: A = M.transpose() sage: A.nullity() 1
- matrix_window(row=0, col=0, nrows=-1, ncols=-1, check=1)#
Return the requested matrix window.
EXAMPLES:
sage: A = matrix(QQ, 3, 3, range(9)) sage: A.matrix_window(1,1, 2, 1) Matrix window of size 2 x 1 at (1,1): [0 1 2] [3 4 5] [6 7 8]
We test the optional check flag.
sage: matrix([1]).matrix_window(0,1,1,1, check=False) Matrix window of size 1 x 1 at (0,1): [1] sage: matrix([1]).matrix_window(0,1,1,1) Traceback (most recent call last): ... IndexError: matrix window index out of range
Another test of bounds checking:
sage: matrix([1]).matrix_window(1,1,1,1) Traceback (most recent call last): ... IndexError: matrix window index out of range
- maxspin(v)#
Computes the largest integer n such that the list of vectors \(S=[v, v*A, ..., v * A^n]\) are linearly independent, and returns that list.
INPUT:
self
- Matrixv
- Vector
OUTPUT:
list
- list of Vectors
ALGORITHM: The current implementation just adds vectors to a vector space until the dimension doesn’t grow. This could be optimized by directly using matrices and doing an efficient Echelon form. Also, when the base is Q, maybe we could simultaneously keep track of what is going on in the reduction modulo p, which might make things much faster.
EXAMPLES:
sage: t = matrix(QQ, 3, 3, range(9)); t [0 1 2] [3 4 5] [6 7 8] sage: v = (QQ^3).0 sage: t.maxspin(v) [(1, 0, 0), (0, 1, 2), (15, 18, 21)] sage: k = t.kernel(); k Vector space of degree 3 and dimension 1 over Rational Field Basis matrix: [ 1 -2 1] sage: t.maxspin(k.0) [(1, -2, 1)]
- minimal_polynomial(var='x', **kwds)#
This is a synonym for
self.minpoly
EXAMPLES:
sage: a = matrix(QQ, 4, 4, range(16)) sage: a.minimal_polynomial('z') z^3 - 30*z^2 - 80*z sage: a.minpoly() x^3 - 30*x^2 - 80*x
- minors(k)#
Return the list of all \(k \times k\) minors of self.
Let \(A\) be an \(m \times n\) matrix and \(k\) an integer with \(0 \leq k\), \(k \leq m\) and \(k \leq n\). A \(k \times k\) minor of \(A\) is the determinant of a \(k \times k\) matrix obtained from \(A\) by deleting \(m - k\) rows and \(n - k\) columns. There are no \(k \times k\) minors of \(A\) if \(k\) is larger than either \(m\) or \(n\).
The returned list is sorted in lexicographical row major ordering, e.g., if A is a \(3 \times 3\) matrix then the minors returned are with these rows/columns: [ [0, 1]x[0, 1], [0, 1]x[0, 2], [0, 1]x[1, 2], [0, 2]x[0, 1], [0, 2]x[0, 2], [0, 2]x[1, 2], [1, 2]x[0, 1], [1, 2]x[0, 2], [1, 2]x[1, 2] ].
INPUT:
k
– integer
EXAMPLES:
sage: A = Matrix(ZZ,2,3,[1,2,3,4,5,6]); A [1 2 3] [4 5 6] sage: A.minors(2) [-3, -6, -3] sage: A.minors(1) [1, 2, 3, 4, 5, 6] sage: A.minors(0) [1] sage: A.minors(5) []
sage: k = GF(37) sage: P.<x0,x1,x2> = PolynomialRing(k) sage: A = Matrix(P,2,3,[x0*x1, x0, x1, x2, x2 + 16, x2 + 5*x1 ]) sage: A.minors(2) [x0*x1*x2 + 16*x0*x1 - x0*x2, 5*x0*x1^2 + x0*x1*x2 - x1*x2, 5*x0*x1 + x0*x2 - x1*x2 - 16*x1]
This test addresses an issue raised at github issue #20512:
sage: A.minors(0)[0].parent() == P True
- minpoly(var='x', **kwds)#
Return the minimal polynomial of self.
This uses a simplistic - and potentially very very slow - algorithm that involves computing kernels to determine the powers of the factors of the charpoly that divide the minpoly.
EXAMPLES:
sage: A = matrix(GF(9,'c'), 4, [1, 1, 0,0, 0,1,0,0, 0,0,5,0, 0,0,0,5]) sage: factor(A.minpoly()) (x + 1) * (x + 2)^2 sage: A.minpoly()(A) == 0 True sage: factor(A.charpoly()) (x + 1)^2 * (x + 2)^2
The default variable name is \(x\), but you can specify another name:
sage: factor(A.minpoly('y')) (y + 1) * (y + 2)^2
- norm(p=2)#
Return the p-norm of this matrix, where \(p\) can be 1, 2, \(\inf\), or the Frobenius norm.
INPUT:
self
- a matrix whose entries are coercible into CDFp
- one of the following options:1
- the largest column-sum norm2 (default)
- the Euclidean normInfinity
- the largest row-sum norm'frob'
- the Frobenius (sum of squares) norm
OUTPUT: RDF number
See also
EXAMPLES:
sage: A = matrix(ZZ, [[1,2,4,3],[-1,0,3,-10]]) sage: A.norm(1) 13.0 sage: A.norm(Infinity) 14.0 sage: B = random_matrix(QQ, 20, 21) sage: B.norm(Infinity) == (B.transpose()).norm(1) True
sage: Id = identity_matrix(12) sage: Id.norm(2) 1.0 sage: A = matrix(RR, 2, 2, [13,-4,-4,7]) sage: A.norm() # rel tol 2e-16 14.999999999999998
Norms of numerical matrices over high-precision reals are computed by this routine. Faster routines for double precision entries from \(RDF\) or \(CDF\) are provided by the
Matrix_double_dense
class.sage: A = matrix(CC, 2, 3, [3*I,4,1-I,1,2,0]) sage: A.norm('frob') 5.656854249492381 sage: A.norm(2) 5.470684443210... sage: A.norm(1) 6.0 sage: A.norm(Infinity) 8.414213562373096 sage: a = matrix([[],[],[],[]]) sage: a.norm() 0.0 sage: a.norm(Infinity) == a.norm(1) True
- nullity()#
Return the (left) nullity of this matrix, which is the dimension of the (left) kernel of this matrix acting from the right on row vectors.
EXAMPLES:
sage: M = Matrix(QQ,[[1,0,0,1],[0,1,1,0],[1,1,1,0]]) sage: M.nullity() 0 sage: M.left_nullity() 0
sage: A = M.transpose() sage: A.nullity() 1 sage: A.left_nullity() 1
sage: M = M.change_ring(ZZ) sage: M.nullity() 0 sage: A = M.transpose() sage: A.nullity() 1
- numerical_approx(prec=None, digits=None, algorithm=None)#
Return a numerical approximation of
self
withprec
bits (or decimaldigits
) of precision.INPUT:
prec
– precision in bitsdigits
– precision in decimal digits (only used ifprec
is not given)algorithm
– ignored for matrices
OUTPUT: A matrix converted to a real or complex field
EXAMPLES:
sage: d = matrix([[3, 0],[0,sqrt(2)]]) sage: b = matrix([[1, -1], [2, 2]]) ; e = b * d * b.inverse();e [ 1/2*sqrt(2) + 3/2 -1/4*sqrt(2) + 3/4] [ -sqrt(2) + 3 1/2*sqrt(2) + 3/2]
sage: e.numerical_approx(53) [ 2.20710678118655 0.396446609406726] [ 1.58578643762690 2.20710678118655]
sage: e.numerical_approx(20) [ 2.2071 0.39645] [ 1.5858 2.2071]
sage: (e-I).numerical_approx(20) [2.2071 - 1.0000*I 0.39645] [ 1.5858 2.2071 - 1.0000*I]
sage: M=matrix(QQ,4,[i/(i+1) for i in range(12)]);M [ 0 1/2 2/3] [ 3/4 4/5 5/6] [ 6/7 7/8 8/9] [ 9/10 10/11 11/12]
sage: M.numerical_approx() [0.000000000000000 0.500000000000000 0.666666666666667] [0.750000000000000 0.800000000000000 0.833333333333333] [0.857142857142857 0.875000000000000 0.888888888888889] [0.900000000000000 0.909090909090909 0.916666666666667]
sage: matrix(SR, 2, 2, range(4)).n() # optional - sage.symbolic [0.000000000000000 1.00000000000000] [ 2.00000000000000 3.00000000000000]
sage: numerical_approx(M) [0.000000000000000 0.500000000000000 0.666666666666667] [0.750000000000000 0.800000000000000 0.833333333333333] [0.857142857142857 0.875000000000000 0.888888888888889] [0.900000000000000 0.909090909090909 0.916666666666667]
We check that github issue #29700 is fixed:
sage: M = matrix(3,[1,1,1,1,0,0,0,1,0]) sage: A,B = M.diagonalization(QQbar) sage: _ = A.n()
- permanent(algorithm='Ryser')#
Return the permanent of this matrix.
Let \(A = (a_{i,j})\) be an \(m \times n\) matrix over any commutative ring with \(m \le n\). The permanent of \(A\) is
\[\mathrm{per}(A) = \sum_\pi a_{1,\pi(1)} a_{2,\pi(2)} \cdots a_{m,\pi(m)}\]where the summation extends over all one-to-one functions \(\pi\) from \(\{1, \ldots, m\}\) to \(\{1, \ldots, n\}\).
The product \(a_{1,\pi(1)} a_{2,\pi(2)} \cdots a_{m,\pi(m)}\) is called diagonal product. So the permanent of an \(m \times n\) matrix \(A\) is the sum of all the diagonal products of \(A\).
By default, this method uses Ryser’s algorithm, but setting
algorithm
to “ButeraPernici” you can use the algorithm of Butera and Pernici (which is well suited for band matrices, i.e. matrices whose entries are concentrated near the diagonal).INPUT:
A
– matrix of size \(m \times n\) with \(m \leq n\)algorithm
– either “Ryser” (default) or “ButeraPernici”. The Butera-Pernici algorithm takes advantage of presence of zeros and is very well suited for sparse matrices.
ALGORITHM:
The Ryser algorithm is implemented in the method
_permanent_ryser()
. It is a modification of theorem 7.1.1. from Brualdi and Ryser: Combinatorial Matrix Theory. Instead of deleting columns from \(A\), we choose columns from \(A\) and calculate the product of the row sums of the selected submatrix.The Butera-Pernici algorithm is implemented in the function
permanental_minor_polynomial()
. It takes advantage of cancellations that may occur in the computations.EXAMPLES:
sage: A = ones_matrix(4,4) sage: A.permanent() 24 sage: A = matrix(3,6,[1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1]) sage: A.permanent() 36 sage: B = A.change_ring(RR) sage: B.permanent() 36.0000000000000
The permanent above is directed to the Sloane’s sequence OEIS sequence A079908 (“The Dancing School Problems”) for which the third term is 36:
sage: oeis(79908) # optional -- internet A079908: Solution to the Dancing School Problem with 3 girls and n+3 boys: f(3,n). sage: _(3) # optional -- internet 36
sage: A = matrix(4,5,[1,1,0,1,1,0,1,1,1,1,1,0,1,0,1,1,1,0,1,0]) sage: A.permanent() 32
A huge permanent that cannot be reasonably computed with the Ryser algorithm (a \(50 \times 50\) band matrix with width \(5\)):
sage: n, w = 50, 5 sage: A = matrix(ZZ, n, n, lambda i,j: (i+j)%5 + 1 if abs(i-j) <= w else 0) sage: A.permanent(algorithm="ButeraPernici") 57766972735511097036962481710892268404670105604676932908
See Minc: Permanents, Example 2.1, p. 5.
sage: A = matrix(QQ,2,2,[1/5,2/7,3/2,4/5]) sage: A.permanent() 103/175
sage: R.<a> = PolynomialRing(ZZ) sage: A = matrix(R,2,2,[a,1,a,a+1]) sage: A.permanent() a^2 + 2*a
sage: R.<x,y> = PolynomialRing(ZZ,2) sage: A = matrix(R,2,2,[x, y, x^2, y^2]) sage: A.permanent() x^2*y + x*y^2
- permanental_minor(k, algorithm='Ryser')#
Return the permanental \(k\)-minor of this matrix.
The permanental \(k\)-minor of a matrix \(A\) is the sum of the permanents of all possible \(k\) by \(k\) submatrices of \(A\). Note that the maximal permanental minor is just the permanent.
For a (0,1)-matrix \(A\) the permanental \(k\)-minor counts the number of different selections of \(k\) 1’s of \(A\) with no two of the 1’s on the same row and no two of the 1’s on the same column.
See Brualdi and Ryser: Combinatorial Matrix Theory, p. 203. Note the typo \(p_0(A) = 0\) in that reference! For applications see Theorem 7.2.1 and Theorem 7.2.4.
See also
The method
rook_vector()
returns the list of all permanental minors.INPUT:
k
– the size of the minoralgorithm
– either “Ryser” (default) or “ButeraPernici”. The Butera-Pernici algorithm is well suited for band matrices.
EXAMPLES:
sage: A = matrix(4,[1,0,1,0,1,0,1,0,1,0,10,10,1,0,1,1]) sage: A.permanental_minor(2) 114
sage: A = matrix(3,6,[1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1]) sage: A.permanental_minor(0) 1 sage: A.permanental_minor(1) 12 sage: A.permanental_minor(2) 40 sage: A.permanental_minor(3) 36
Note that if \(k = m = n\), the permanental \(k\)-minor equals \(\mathrm{per}(A)\):
sage: A.permanent() 36
The permanental minors of the “complement” matrix of \(A\) is related to the permanent of \(A\):
sage: m, n = 3, 6 sage: C = matrix(m, n, lambda i,j: 1 - A[i,j]) sage: sum((-1)^k * C.permanental_minor(k)*factorial(n-k)/factorial(n-m) for k in range(m+1)) 36
See Theorem 7.2.1 of Brualdi and Ryser: Combinatorial Matrix Theory: per(A)
- permutation_normal_form(check=False)#
Take the set of matrices that are
self
permuted by any row and column permutation, and return the maximal one of the set where matrices are ordered lexicographically going along each row.INPUT:
check
– (default:False
) IfTrue
return a tuple ofthe maximal matrix and the permutations taking
self
to the maximal matrix. IfFalse
, return only the maximal matrix.
OUTPUT:
The maximal matrix.
EXAMPLES:
sage: M = matrix(ZZ, [[0, 0, 1], [1, 0, 2], [0, 0, 0]]) sage: M [0 0 1] [1 0 2] [0 0 0] sage: M.permutation_normal_form() [2 1 0] [1 0 0] [0 0 0] sage: M = matrix(ZZ, [[-1, 3], [-1, 5], [2, 4]]) sage: M [-1 3] [-1 5] [ 2 4] sage: M.permutation_normal_form(check=True) ( [ 5 -1] [ 4 2] [ 3 -1], ((1,2,3), (1,2)) )
- pfaffian(algorithm=None, check=True)#
Return the Pfaffian of
self
, assuming thatself
is an alternating matrix.The result is cached.
INPUT:
algorithm
(default:None
) – string, the algorithm to use; currently the following algorithms have been implemented:'bfl'
- using the Bär-Faddeev-LeVerrier algorithm'definition'
- using the definition given by perfect matchings
check
(default:True
) – boolean determining whether to checkself
for alternatingness and squareness. This has to be set toFalse
ifself
is defined over a non-discrete ring.
The Pfaffian of an alternating matrix is defined as follows:
Let \(A\) be an alternating \(k \times k\) matrix over a commutative ring. (Here, “alternating” means that \(A^T = -A\) and that the diagonal entries of \(A\) are zero.) If \(k\) is odd, then the Pfaffian of the matrix \(A\) is defined to be \(0\). Let us now define it when \(k\) is even. In this case, set \(n = k/2\) (this is an integer). For every \(i\) and \(j\), we denote the \((i, j)\)-th entry of \(A\) by \(a_{i, j}\). Let \(M\) denote the set of all perfect matchings of the set \(\{ 1, 2, \ldots, 2n \}\) (see
sage.combinat.perfect_matching.PerfectMatchings
). For every matching \(m \in M\), define the sign \(\mathrm{sign}(m)\) of \(m\) by writing \(m\) as \(\{ \{ i_1, j_1 \}, \{ i_2, j_2 \}, \ldots, \{ i_n, j_n \} \}\) with \(i_k < j_k\) for all \(k\), and setting \(\mathrm{sign}(m)\) to be the sign of the permutation \(( i_1, j_1, i_2, j_2, \ldots, i_n, j_n )\) (written here in one-line notation). For every matching \(m \in M\), define the weight \(w(m)\) of \(m\) by writing \(m\) as \(\{ \{ i_1, j_1 \}, \{ i_2, j_2 \}, \ldots, \{ i_n, j_n \} \}\) with \(i_k < j_k\) for all \(k\), and setting \(w(m) = a_{i_1, j_1} a_{i_2, j_2} \cdots a_{i_n, j_n}\). Now, the Pfaffian of the matrix \(A\) is defined to be the sum\[\sum_{m \in M} \mathrm{sign}(m) w(m).\]The Pfaffian of \(A\) is commonly denoted by \(\mathrm{Pf}(A)\). It is well-known that \((\mathrm{Pf}(A))^2 = \det A\) for every alternating matrix \(A\), and that \(\mathrm{Pf} (U^T A U) = \det U \cdot \mathrm{Pf}(A)\) for any \(n \times n\) matrix \(U\) and any alternating \(n \times n\) matrix \(A\).
See [Knu1995], [DW1995] and [Rot2001], [Baer2020], just to name a few sources, for further properties of Pfaffians.
ALGORITHM:
If the matrix is small, namely up to size \(4 \times 4\), the naive formulas are always used.
The Bär-Faddeev-LeVerrier algorithm can be accessed using
'bfl'
. It works over any \(\QQ\)-algebra or ring whose fraction field is an \(\QQ\)-algebra (see [Baer2020] for details). If that check fails, the implementation raises an error because correct results cannot be guaranteed.To access the algorithm using the above definition, use
'definition'
. However, notice that this algorithm is usually very slow.By default, i.e. if no options are set, the implementation tries to apply the BFL algorithm first. If BFL is not applicable, it uses the definition by perfect matchings.
The alternatingness of the matrix
self
is checked only ifcheck
isTrue
(this is important because even ifself
is alternating, a non-discrete base ring might prevent Sage from being able to check this).EXAMPLES:
A \(3 \times 3\) alternating matrix has Pfaffian 0 independently of its entries:
sage: MSp = MatrixSpace(Integers(27), 3) sage: A = MSp([0, 2, -3, -2, 0, 8, 3, -8, 0]) sage: A.pfaffian() 0 sage: parent(A.pfaffian()) Ring of integers modulo 27
The Pfaffian of a \(2 \times 2\) alternating matrix is just its northeast entry:
sage: MSp = MatrixSpace(QQ, 2) sage: A = MSp([0, 4, -4, 0]) sage: A.pfaffian() 4 sage: parent(A.pfaffian()) Rational Field
The Pfaffian of a \(0 \times 0\) alternating matrix is \(1\):
sage: MSp = MatrixSpace(ZZ, 0) sage: A = MSp([]) sage: A.pfaffian() 1 sage: parent(A.pfaffian()) Integer Ring
Let us compute the Pfaffian of a generic \(4 \times 4\) alternating matrix:
sage: R = PolynomialRing(QQ, 'x12,x13,x14,x23,x24,x34') sage: x12, x13, x14, x23, x24, x34 = R.gens() sage: A = matrix(R, [[ 0, x12, x13, x14], ....: [-x12, 0, x23, x24], ....: [-x13, -x23, 0, x34], ....: [-x14, -x24, -x34, 0]]) sage: A.pfaffian() x14*x23 - x13*x24 + x12*x34 sage: parent(A.pfaffian()) Multivariate Polynomial Ring in x12, x13, x14, x23, x24, x34 over Rational Field
The Pfaffian of an alternating matrix squares to its determinant:
sage: A = [[0] * 6 for i in range(6)] sage: for i in range(6): ....: for j in range(i): ....: u = floor(random() * 10) ....: A[i][j] = u ....: A[j][i] = -u ....: A[i][i] = 0 sage: AA = Matrix(ZZ, A) sage: AA.pfaffian() ** 2 == AA.det() True
In order to use the Bär-Faddeev-LeVerrier algorithm, the base ring must have characteristic zero:
sage: A = matrix(GF(5), [(0, 3, 4, 1, 3, 4), ....: (2, 0, 2, 0, 1, 0), ....: (1, 3, 0, 4, 1, 0), ....: (4, 0, 1, 0, 2, 0), ....: (2, 4, 4, 3, 0, 0), ....: (1, 0, 0, 0, 0, 0)]) sage: A.pfaffian(algorithm='bfl') Traceback (most recent call last): ... TypeError: Bär-Faddeev-LeVerrier algorithm not applicable, use another algorithm instead
In that case, the definition by perfect matchings is used instead:
sage: A.pfaffian() 2
- pivot_rows()#
Return the pivot row positions for this matrix, which are a topmost subset of the rows that span the row space and are linearly independent.
OUTPUT: a tuple of integers
EXAMPLES:
sage: A = matrix(QQ,3,3, [0,0,0,1,2,3,2,4,6]); A [0 0 0] [1 2 3] [2 4 6] sage: A.pivot_rows() (1,) sage: A.pivot_rows() # testing cached value (1,)
- plot(*args, **kwds)#
A plot of this matrix.
Each (ith, jth) matrix element is given a different color value depending on its relative size compared to the other elements in the matrix.
The tick marks drawn on the frame axes denote the (ith, jth) element of the matrix.
This method just calls
matrix_plot
.*args
and**kwds
are passed tomatrix_plot
.EXAMPLES:
A matrix over ZZ colored with different grey levels:
sage: A = matrix([[1,3,5,1],[2,4,5,6],[1,3,5,7]]) sage: A.plot() Graphics object consisting of 1 graphics primitive
Here we make a random matrix over RR and use cmap=’hsv’ to color the matrix elements different RGB colors (see documentation for
matrix_plot
for more information on cmaps):sage: A = random_matrix(RDF, 50) sage: plot(A, cmap='hsv') Graphics object consisting of 1 graphics primitive
Another random plot, but over GF(389):
sage: A = random_matrix(GF(389), 10) sage: A.plot(cmap='Oranges') Graphics object consisting of 1 graphics primitive
- principal_square_root(check_positivity=True)#
Return the principal square root of a positive definite matrix.
A positive definite matrix \(A\) has a unique positive definite matrix \(M\) such that \(M^2 = A\).
See Wikipedia article Square_root_of_a_matrix.
EXAMPLES:
sage: A = Matrix([[1,-1/2,0],[-1/2,1,-1/2],[0,-1/2,1]]) sage: B = A.principal_square_root() sage: A == B^2 True
- prod_of_row_sums(cols)#
Calculate the product of all row sums of a submatrix of \(A\) for a list of selected columns
cols
.EXAMPLES:
sage: a = matrix(QQ, 2,2, [1,2,3,2]); a [1 2] [3 2] sage: a.prod_of_row_sums([0,1]) 15
Another example:
sage: a = matrix(QQ, 2,3, [1,2,3,2,5,6]); a [1 2 3] [2 5 6] sage: a.prod_of_row_sums([1,2]) 55
- pseudoinverse(algorithm=None)#
Return the Moore-Penrose pseudoinverse of this matrix.
INPUT:
algorithm
(default: guess) – one of the following:"numpy"
– Use numpy’slinalg.pinv()
which is suitable over real or complex fields."exact"
– Use a simple algorithm which is not numerically stable but useful over exact fields. Assume that no conjugation is needed, that the conjugate transpose is just the transpose."exactconj"
– Likeexact
but use the conjugate transpose.
OUTPUT: a matrix
EXAMPLES:
sage: M = diagonal_matrix(CDF, [0, I, 1+I]) sage: M [ 0.0 0.0 0.0] [ 0.0 1.0*I 0.0] [ 0.0 0.0 1.0 + 1.0*I] sage: M.pseudoinverse() # tol 1e-15 [ 0.0 0.0 0.0] [ 0.0 -1.0*I 0.0] [ 0.0 0.0 0.5 - 0.5*I]
We check the properties of the pseudoinverse over an exact field:
sage: M = random_matrix(QQ, 6, 3) * random_matrix(QQ, 3, 5) sage: Mx = M.pseudoinverse() sage: M * Mx * M == M True sage: Mx * M * Mx == Mx True sage: (M * Mx).is_symmetric() True sage: (Mx * M).is_symmetric() True
Beware that the
exact
algorithm is not numerically stable, but the defaultnumpy
algorithm is:sage: M = matrix.hilbert(12,ring=RR) sage: (~M*M).norm() # a considerable error 1.3... sage: Mx = M.pseudoinverse(algorithm="exact") sage: (Mx*M).norm() # huge error 11.5... sage: Mx = M.pseudoinverse(algorithm="numpy") sage: (Mx*M).norm() # still OK 1.00...
When multiplying the given matrix with the pseudoinverse, the result is symmetric for the
exact
algorithm or hermitian for theexactconj
algorithm:sage: M = matrix(QQbar, 2, 2, [1, sqrt(-3), -sqrt(-3), 3]) sage: M * M.pseudoinverse() [ 0.2500000000000000? 0.4330127018922193?*I] [-0.4330127018922193?*I 0.750000000000000?] sage: M * M.pseudoinverse(algorithm="exactconj") [ 1/4 0.4330127018922193?*I] [-0.4330127018922193?*I 3/4] sage: M * M.pseudoinverse(algorithm="exact") [ -1/2 0.866025403784439?*I] [0.866025403784439?*I 3/2]
For an invertible matrix, the pseudoinverse is just the inverse:
sage: M = matrix([[1,2], [3,4]]) sage: ~M [ -2 1] [ 3/2 -1/2] sage: M.pseudoinverse() [ -2 1] [ 3/2 -1/2]
Numpy gives a strange answer due to rounding errors:
sage: M.pseudoinverse(algorithm="numpy") # random [-1286742750677287/643371375338643 1000799917193445/1000799917193444] [ 519646110850445/346430740566963 -300239975158034/600479950316067]
Although it is not too far off:
sage: (~M-M.pseudoinverse(algorithm="numpy")).norm() < 1e-14 True
- qdet(q=None)#
Return the quantum deteminant of
self
.The quantum determinant of a matrix \(M = (m_{ij})_{i,j=1}^n\) is defined by
\[\det_q(M) = \sum_{\sigma \in S_n} (-q)^{\ell(\sigma)} M_{\sigma(i),j},\]where \(S_n\) is the symmetric group on \(\{1, \ldots, n\}\) and \(\ell(\sigma)\) denotes the length of \(\sigma\) written as simple transpositions (equivalently the number of inversions when written in one-line notation).
INPUT:
q
– the parameter \(q\); the default is \(q \in F[q]\), where \(F\) is the base ring ofself
EXAMPLES:
sage: A = matrix(SR, 2, lambda i, j: f'a{i}{j}'); A # optional - sage.symbolic [a00 a01] [a10 a11] sage: A.quantum_determinant() # optional - sage.symbolic -a01*a10*q + a00*a11 sage: A = matrix(SR, 3, lambda i, j: f'a{i}{j}') # optional - sage.symbolic sage: A.quantum_determinant() # optional - sage.symbolic -a02*a11*a20*q^3 + (a01*a12*a20 + a02*a10*a21)*q^2 + (-a00*a12*a21 - a01*a10*a22)*q + a00*a11*a22 sage: R.<q> = LaurentPolynomialRing(ZZ) sage: MS = MatrixSpace(Integers(8), 3) sage: A = MS([1,7,3, 1,1,1, 3,4,5]) sage: A.det() 6 sage: A.quantum_determinant(q^-2) 7*q^-6 + q^-4 + q^-2 + 5 sage: S.<x,y> = PolynomialRing(GF(7)) sage: R.<q> = LaurentPolynomialRing(S) sage: MS = MatrixSpace(S, 3, sparse=True) sage: A = MS([[x, y, 3], [4, 2+y, x^2], [0, 1-x, x+y]]) sage: A.det() x^4 - x^3 + x^2*y + x*y^2 + 2*x^2 - 2*x*y + 3*y^2 + 2*x - 2 sage: A.quantum_determinant() (2*x - 2)*q^2 + (x^4 - x^3 + 3*x*y + 3*y^2)*q + x^2*y + x*y^2 + 2*x^2 + 2*x*y sage: A.quantum_determinant(int(2)) 2*x^4 - 2*x^3 + x^2*y + x*y^2 + 2*x^2 + x*y - y^2 + x - 1 sage: A.quantum_determinant(q*x + q^-1*y) (2*x*y^2 - 2*y^2)*q^-2 + (x^4*y - x^3*y + 3*x*y^2 + 3*y^3)*q^-1 + (-2*x^2*y + x*y^2 + 2*x^2 - 2*x*y) + (x^5 - x^4 + 3*x^2*y + 3*x*y^2)*q + (2*x^3 - 2*x^2)*q^2
- quantum_determinant(q=None)#
Return the quantum deteminant of
self
.The quantum determinant of a matrix \(M = (m_{ij})_{i,j=1}^n\) is defined by
\[\det_q(M) = \sum_{\sigma \in S_n} (-q)^{\ell(\sigma)} M_{\sigma(i),j},\]where \(S_n\) is the symmetric group on \(\{1, \ldots, n\}\) and \(\ell(\sigma)\) denotes the length of \(\sigma\) written as simple transpositions (equivalently the number of inversions when written in one-line notation).
INPUT:
q
– the parameter \(q\); the default is \(q \in F[q]\), where \(F\) is the base ring ofself
EXAMPLES:
sage: A = matrix(SR, 2, lambda i, j: f'a{i}{j}'); A # optional - sage.symbolic [a00 a01] [a10 a11] sage: A.quantum_determinant() # optional - sage.symbolic -a01*a10*q + a00*a11 sage: A = matrix(SR, 3, lambda i, j: f'a{i}{j}') # optional - sage.symbolic sage: A.quantum_determinant() # optional - sage.symbolic -a02*a11*a20*q^3 + (a01*a12*a20 + a02*a10*a21)*q^2 + (-a00*a12*a21 - a01*a10*a22)*q + a00*a11*a22 sage: R.<q> = LaurentPolynomialRing(ZZ) sage: MS = MatrixSpace(Integers(8), 3) sage: A = MS([1,7,3, 1,1,1, 3,4,5]) sage: A.det() 6 sage: A.quantum_determinant(q^-2) 7*q^-6 + q^-4 + q^-2 + 5 sage: S.<x,y> = PolynomialRing(GF(7)) sage: R.<q> = LaurentPolynomialRing(S) sage: MS = MatrixSpace(S, 3, sparse=True) sage: A = MS([[x, y, 3], [4, 2+y, x^2], [0, 1-x, x+y]]) sage: A.det() x^4 - x^3 + x^2*y + x*y^2 + 2*x^2 - 2*x*y + 3*y^2 + 2*x - 2 sage: A.quantum_determinant() (2*x - 2)*q^2 + (x^4 - x^3 + 3*x*y + 3*y^2)*q + x^2*y + x*y^2 + 2*x^2 + 2*x*y sage: A.quantum_determinant(int(2)) 2*x^4 - 2*x^3 + x^2*y + x*y^2 + 2*x^2 + x*y - y^2 + x - 1 sage: A.quantum_determinant(q*x + q^-1*y) (2*x*y^2 - 2*y^2)*q^-2 + (x^4*y - x^3*y + 3*x*y^2 + 3*y^3)*q^-1 + (-2*x^2*y + x*y^2 + 2*x^2 - 2*x*y) + (x^5 - x^4 + 3*x^2*y + 3*x*y^2)*q + (2*x^3 - 2*x^2)*q^2
- randomize(density=1, nonzero=False, *args, **kwds)#
Replace a proportion of the entries of a matrix by random elements, leaving the remaining entries unchanged.
Note
The locations of the entries of the matrix to change are determined randomly, with the total number of locations determined by the
density
keyword. These locations are not guaranteed to be distinct. So it is possible that the same position can be chosen multiple times, especially for a very small matrix. The exception is whendensity = 1
, in which case every entry of the matrix will be changed.INPUT:
density
-float
(default:1
); upper bound for the proportion of entries that are changednonzero
- Bool (default:False
); ifTrue
, then new entries will be nonzero*args, **kwds
- Remaining parameters may be passed to therandom_element
function of the base ring
EXAMPLES:
We construct the zero matrix over a polynomial ring.
sage: a = matrix(QQ['x'], 3); a [0 0 0] [0 0 0] [0 0 0]
We then randomize roughly half the entries:
sage: a.randomize(0.5) sage: a.density() <= 0.5 True
Now we randomize all the entries of the resulting matrix:
sage: while a.density() < 0.9: ....: a = matrix(QQ['x'], 3) ....: a.randomize()
We create the zero matrix over the integers:
sage: a = matrix(ZZ, 2); a [0 0] [0 0]
Then we randomize it; the
x
andy
keywords, which determine the size of the random elements, are passed on to therandom_element
method forZZ
.sage: a.randomize(x=-2^64, y=2^64) sage: while all(abs(b) < 2^63 for b in a.list()): ....: a.randomize(x=-2^64, y=2^64) sage: all(abs(b) < 2^64 for b in a.list()) True
- rational_form(format='right', subdivide=True)#
Returns the rational canonical form, also known as Frobenius form.
INPUT:
self
- a square matrix with entries from an exact field.format
- default: ‘right’ - one of ‘right’, ‘bottom’, ‘left’, ‘top’ or ‘invariants’. The first four will cause a matrix to be returned with companion matrices dictated by the keyword. The value ‘invariants’ will cause a list of lists to be returned, where each list contains coefficients of a polynomial associated with a companion matrix.subdivide
- default: ‘True’ - if ‘True’ and a matrix is returned, then it contains subdivisions delineating the companion matrices along the diagonal.
OUTPUT:
The rational form of a matrix is a similar matrix composed of submatrices (“blocks”) placed on the main diagonal. Each block is a companion matrix. Associated with each companion matrix is a polynomial. In rational form, the polynomial of one block will divide the polynomial of the next block (and thus, the polynomials of all subsequent blocks).
Rational form, also known as Frobenius form, is a canonical form. In other words, two matrices are similar if and only if their rational canonical forms are equal. The algorithm used does not provide the similarity transformation matrix (also known as the change-of-basis matrix).
Companion matrices may be written in one of four styles, and any such style may be selected with the
format
keyword. See the companion matrix constructor,sage.matrix.constructor.companion_matrix()
, for more information about companion matrices.If the ‘invariants’ value is used for the
format
keyword, then the return value is a list of lists, where each list is the coefficients of the polynomial associated with one of the companion matrices on the diagonal. These coefficients include the leading one of the monic polynomial and are ready to be coerced into any polynomial ring over the same field (see examples of this below). This return value is intended to be the most compact representation and the easiest to use for testing equality of rational forms.Because the minimal and characteristic polynomials of a companion matrix are the associated polynomial, it is easy to see that the product of the polynomials of the blocks will be the characteristic polynomial and the final polynomial will be the minimal polynomial of the entire matrix.
ALGORITHM:
We begin with ZigZag form, which is due to Arne Storjohann and is documented at
zigzag_form()
. Then we eliminate ‘’corner’’ entries enroute to rational form via an additional algorithm of Storjohann’s [Sto2011].EXAMPLES:
The lists of coefficients returned with the
invariants
keyword are designed to easily convert to the polynomials associated with the companion matrices. This is illustrated by the construction below of thepolys
list. Then we can test the divisibility condition on the list of polynomials. Also the minimal and characteristic polynomials are easy to determine from this list.sage: A = matrix(QQ, [[ 11, 14, -15, -4, -38, -29, 1, 23, 14, -63, 17, 24, 36, 32], ....: [ 18, 6, -17, -11, -31, -43, 12, 26, 0, -69, 11, 13, 17, 24], ....: [ 11, 16, -22, -8, -48, -34, 0, 31, 16, -82, 26, 31, 39, 37], ....: [ -8, -18, 22, 10, 46, 33, 3, -27, -12, 70, -19, -20, -42, -31], ....: [-13, -21, 16, 10, 52, 43, 4, -28, -25, 89, -37, -20, -53, -62], ....: [ -2, -6, 0, 0, 6, 10, 1, 1, -7, 14, -11, -3, -10, -18], ....: [ -9, -19, -3, 4, 23, 30, 8, -3, -27, 55, -40, -5, -40, -69], ....: [ 4, -8, -1, -1, 5, -4, 9, 5, -11, 4, -14, -2, -13, -17], ....: [ 1, -2, 16, -1, 19, -2, -1, -17, 2, 19, 5, -25, -7, 14], ....: [ 7, 7, -13, -4, -26, -21, 3, 18, 5, -40, 7, 15, 20, 14], ....: [ -6, -7, -12, 4, -1, 18, 3, 8, -11, 15, -18, 17, -15, -41], ....: [ 5, 11, -11, -3, -26, -19, -1, 14, 10, -42, 14, 17, 25, 23], ....: [-16, -15, 3, 10, 29, 45, -1, -13, -19, 71, -35, -2, -35, -65], ....: [ 4, 2, 3, -2, -2, -10, 1, 0, 3, -11, 6, -4, 6, 17]]) sage: A.rational_form() [ 0 -4| 0 0 0 0| 0 0 0 0 0 0 0 0] [ 1 4| 0 0 0 0| 0 0 0 0 0 0 0 0] [---------+-------------------+---------------------------------------] [ 0 0| 0 0 0 12| 0 0 0 0 0 0 0 0] [ 0 0| 1 0 0 -4| 0 0 0 0 0 0 0 0] [ 0 0| 0 1 0 -9| 0 0 0 0 0 0 0 0] [ 0 0| 0 0 1 6| 0 0 0 0 0 0 0 0] [---------+-------------------+---------------------------------------] [ 0 0| 0 0 0 0| 0 0 0 0 0 0 0 -216] [ 0 0| 0 0 0 0| 1 0 0 0 0 0 0 108] [ 0 0| 0 0 0 0| 0 1 0 0 0 0 0 306] [ 0 0| 0 0 0 0| 0 0 1 0 0 0 0 -271] [ 0 0| 0 0 0 0| 0 0 0 1 0 0 0 -41] [ 0 0| 0 0 0 0| 0 0 0 0 1 0 0 134] [ 0 0| 0 0 0 0| 0 0 0 0 0 1 0 -64] [ 0 0| 0 0 0 0| 0 0 0 0 0 0 1 13] sage: R = PolynomialRing(QQ, 'x') sage: invariants = A.rational_form(format='invariants') sage: invariants [[4, -4, 1], [-12, 4, 9, -6, 1], [216, -108, -306, 271, 41, -134, 64, -13, 1]] sage: polys = [R(p) for p in invariants] sage: [p.factor() for p in polys] [(x - 2)^2, (x - 3) * (x + 1) * (x - 2)^2, (x + 1)^2 * (x - 3)^3 * (x - 2)^3] sage: all(polys[i].divides(polys[i+1]) for i in range(len(polys)-1)) True sage: polys[-1] == A.minimal_polynomial(var='x') True sage: prod(polys) == A.characteristic_polynomial(var='x') True
Rational form is a canonical form. Any two matrices are similar if and only if their rational forms are equal. By starting with Jordan canonical forms, the matrices
C
andD
below were built as similar matrices, whileE
was built to be just slightly different. All three matrices have equal characteristic polynomials thoughE
’s minimal polynomial differs.sage: C = matrix(QQ, [[2, 31, -10, -9, -125, 13, 62, -12], ....: [0, 48, -16, -16, -188, 20, 92, -16], ....: [0, 9, -1, 2, -33, 5, 18, 0], ....: [0, 15, -5, 0, -59, 7, 30, -4], ....: [0, -21, 7, 2, 84, -10, -42, 5], ....: [0, -42, 14, 8, 167, -17, -84, 13], ....: [0, -50, 17, 10, 199, -23, -98, 14], ....: [0, 15, -5, -2, -59, 7, 30, -2]]) sage: C.minimal_polynomial().factor() (x - 2)^2 sage: C.characteristic_polynomial().factor() (x - 2)^8 sage: C.rational_form() [ 0 -4| 0 0| 0 0| 0 0] [ 1 4| 0 0| 0 0| 0 0] [-----+-----+-----+-----] [ 0 0| 0 -4| 0 0| 0 0] [ 0 0| 1 4| 0 0| 0 0] [-----+-----+-----+-----] [ 0 0| 0 0| 0 -4| 0 0] [ 0 0| 0 0| 1 4| 0 0] [-----+-----+-----+-----] [ 0 0| 0 0| 0 0| 0 -4] [ 0 0| 0 0| 0 0| 1 4] sage: D = matrix(QQ, [[ -4, 3, 7, 2, -4, 5, 7, -3], ....: [ -6, 5, 7, 2, -4, 5, 7, -3], ....: [ 21, -12, 89, 25, 8, 27, 98, -95], ....: [ -9, 5, -44, -11, -3, -13, -48, 47], ....: [ 23, -13, 74, 21, 12, 22, 85, -84], ....: [ 31, -18, 135, 38, 12, 47, 155, -147], ....: [-33, 19, -138, -39, -13, -45, -156, 151], ....: [ -7, 4, -29, -8, -3, -10, -34, 34]]) sage: D.minimal_polynomial().factor() (x - 2)^2 sage: D.characteristic_polynomial().factor() (x - 2)^8 sage: D.rational_form() [ 0 -4| 0 0| 0 0| 0 0] [ 1 4| 0 0| 0 0| 0 0] [-----+-----+-----+-----] [ 0 0| 0 -4| 0 0| 0 0] [ 0 0| 1 4| 0 0| 0 0] [-----+-----+-----+-----] [ 0 0| 0 0| 0 -4| 0 0] [ 0 0| 0 0| 1 4| 0 0] [-----+-----+-----+-----] [ 0 0| 0 0| 0 0| 0 -4] [ 0 0| 0 0| 0 0| 1 4] sage: E = matrix(QQ, [[ 0, -8, 4, -6, -2, 5, -3, 11], ....: [-2, -4, 2, -4, -2, 4, -2, 6], ....: [ 5, 14, -7, 12, 3, -8, 6, -27], ....: [-3, -8, 7, -5, 0, 2, -6, 17], ....: [ 0, 5, 0, 2, 4, -4, 1, 2], ....: [-3, -7, 5, -6, -1, 5, -4, 14], ....: [ 6, 18, -10, 14, 4, -10, 10, -28], ....: [-2, -6, 4, -5, -1, 3, -3, 13]]) sage: E.minimal_polynomial().factor() (x - 2)^3 sage: E.characteristic_polynomial().factor() (x - 2)^8 sage: E.rational_form() [ 2| 0 0| 0 0| 0 0 0] [---+-------+-------+-----------] [ 0| 0 -4| 0 0| 0 0 0] [ 0| 1 4| 0 0| 0 0 0] [---+-------+-------+-----------] [ 0| 0 0| 0 -4| 0 0 0] [ 0| 0 0| 1 4| 0 0 0] [---+-------+-------+-----------] [ 0| 0 0| 0 0| 0 0 8] [ 0| 0 0| 0 0| 1 0 -12] [ 0| 0 0| 0 0| 0 1 6]
The principal feature of rational canonical form is that it can be computed over any field using only field operations. Other forms, such as Jordan canonical form, are complicated by the need to determine the eigenvalues of the matrix, which can lie outside the field. The following matrix has all of its eigenvalues outside the rationals - some are irrational (\(\pm\sqrt{2}\)) and the rest are complex (\(-1\pm 2i\)).
sage: A = matrix(QQ, ....: [[-154, -3, -54, 44, 48, -244, -19, 67, -326, 85, 355, 581], ....: [ 504, 25, 156, -145, -171, 793, 99, -213, 1036, -247, -1152, -1865], ....: [ 294, -1, 112, -89, -90, 469, 36, -128, 634, -160, -695, -1126], ....: [ -49, -32, 25, 7, 37, -64, -58, 12, -42, -14, 72, 106], ....: [-261, -123, 65, 47, 169, -358, -254, 70, -309, -29, 454, 673], ....: [-448, -123, -10, 109, 227, -668, -262, 163, -721, 95, 896, 1410], ....: [ 38, 7, 8, -14, -17, 66, 6, -23, 73, -29, -78, -143], ....: [ -96, 10, -55, 37, 24, -168, 17, 56, -231, 88, 237, 412], ....: [ 310, 67, 31, -81, -143, 473, 143, -122, 538, -98, -641, -1029], ....: [ 139, -35, 99, -49, -18, 236, -41, -70, 370, -118, -377, -619], ....: [ 243, 9, 81, -72, -81, 386, 43, -105, 508, -124, -564, -911], ....: [-155, -3, -55, 45, 50, -245, -27, 65, -328, 77, 365, 583]]) sage: A.characteristic_polynomial().factor() (x^2 - 2)^2 * (x^2 + 2*x + 5)^4 sage: A.eigenvalues(extend=False) [] sage: A.rational_form() [ 0 -5| 0 0 0 0| 0 0 0 0 0 0] [ 1 -2| 0 0 0 0| 0 0 0 0 0 0] [-------+---------------+-----------------------] [ 0 0| 0 0 0 10| 0 0 0 0 0 0] [ 0 0| 1 0 0 4| 0 0 0 0 0 0] [ 0 0| 0 1 0 -3| 0 0 0 0 0 0] [ 0 0| 0 0 1 -2| 0 0 0 0 0 0] [-------+---------------+-----------------------] [ 0 0| 0 0 0 0| 0 0 0 0 0 50] [ 0 0| 0 0 0 0| 1 0 0 0 0 40] [ 0 0| 0 0 0 0| 0 1 0 0 0 3] [ 0 0| 0 0 0 0| 0 0 1 0 0 -12] [ 0 0| 0 0 0 0| 0 0 0 1 0 -12] [ 0 0| 0 0 0 0| 0 0 0 0 1 -4] sage: F.<x> = QQ[] sage: polys = A.rational_form(format='invariants') sage: [F(p).factor() for p in polys] [x^2 + 2*x + 5, (x^2 - 2) * (x^2 + 2*x + 5), (x^2 - 2) * (x^2 + 2*x + 5)^2]
Rational form may be computed over any field. The matrix below is an example where the eigenvalues lie outside the field.
sage: F.<a> = FiniteField(7^2) sage: A = matrix(F, ....: [[5*a + 3, 4*a + 1, 6*a + 2, 2*a + 5, 6, 4*a + 5, 4*a + 5, 5, a + 6, 5, 4*a + 4], ....: [6*a + 3, 2*a + 4, 0, 6, 5*a + 5, 2*a, 5*a + 1, 1, 5*a + 2, 4*a, 5*a + 6], ....: [3*a + 1, 6