Base class for matrices, part 2
For design documentation see matrix/docs.py.
AUTHORS:
TESTS:
sage: m = matrix(ZZ['x'], 2, 3, [1..6])
sage: TestSuite(m).run()
Bases: sage.matrix.matrix1.Matrix
The initialization routine of the Matrix base class ensures that it sets the attributes self._parent, self._base_ring, self._nrows, self._ncols. It sets the latter ones by accessing the relevant information on parent, which is often slower than what a more specific subclass can do.
Subclasses of Matrix can safely skip calling Matrix.__init__ provided they take care of initializing these attributes themselves.
The private attributes self._is_immutable and self._cache are implicitly initialized to valid values upon memory allocation.
EXAMPLES:
sage: import sage.matrix.matrix0
sage: A = sage.matrix.matrix0.Matrix(MatrixSpace(QQ,2))
sage: type(A)
<type 'sage.matrix.matrix0.Matrix'>
Returns the conjugate matrix.
EXAMPLE:
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]
Returns the conjugatetranspose (Hermitian) matrix.
EXAMPLE:
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]
Returns the inverse of the matrix, if it exists.
EXAMPLES:
sage: A = matrix(QQ, [[5, 3, 1, 7],
... [ 1, 1, 1, 0],
... [1, 2, 2, 0],
... [2, 1, 0, 4]])
sage: A.I
[ 0 2 1 0]
[4 8 2 7]
[ 4 7 1 7]
[ 1 1 0 2]
sage: B = matrix(QQ, [[11, 5, 18, 6],
... [ 1, 2, 6, 8],
... [ 4, 2, 7, 3],
... [ 1, 2, 5, 11]])
sage: B.I
Traceback (most recent call last):
...
ZeroDivisionError: input matrix must be nonsingular
Finds a decomposition into a lowertriangular matrix and an uppertriangular matrix.
INPUT:
OUTPUT:
Suppose that \(A\) is an \(m\times n\) matrix, then an LU decomposition is a lowertriangular \(m\times m\) matrix \(L\) with every diagonal element equal to 1, and an uppertriangular \(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
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 nondiagonal 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 rowechelon 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 [TREFETHENBAU].
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 and CDF). Matrices returned are over the original field, or the fraction field of the ring. If the field is not ordered (i.e. the absolute value function is not implemented), then the pivot strategy needs to be ‘nonzero’.
sage: A = matrix(RealField(100), 3, 3, range(9))
sage: P, L, U = A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix must be exact, not Real Field with 100 bits of precision
sage: A = matrix(Integers(6), 3, 2, range(6))
sage: A.LU()
Traceback (most recent call last):
...
TypeError: base ring of the matrix needs a field of fractions, not Ring of integers modulo 6
sage: R.<y> = PolynomialRing(QQ, 'y')
sage: B = matrix(R, [[y+1, y^2+y], [y^2, y^3]])
sage: P, L, U = B.LU(pivot='partial')
Traceback (most recent call last):
...
TypeError: cannot take absolute value of matrix entries, try 'pivot=nonzero'
sage: P, L, U = B.LU(pivot='nonzero')
sage: P
[1 0]
[0 1]
sage: L
[ 1 0]
[y^2/(y + 1) 1]
sage: U
[ y + 1 y^2 + y]
[ 0 0]
sage: L.base_ring()
Fraction Field of Univariate Polynomial Ring in y over Rational Field
sage: B == P*L*U
True
sage: 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 uppertriangular. 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)
TESTS:
Unknown keywords are caught.
sage: A = matrix(ZZ, 2, range(4))
sage: A.LU(pivot='junk')
Traceback (most recent call last):
...
ValueError: pivot strategy must be None, 'partial' or 'nonzero', not junk
sage: A.LU(format='garbage')
Traceback (most recent call last):
...
ValueError: format must be 'plu' or 'compact', not garbage
Components of the ‘compact’ format are immutable, while components of the ‘plu’ format are not.
sage: A = matrix(ZZ, 2, range(4))
sage: perm, M = A.LU(format='compact')
sage: perm[0] = 25
Traceback (most recent call last):
...
TypeError: 'tuple' object does not support item assignment
sage: M.is_immutable()
True
sage: P, L, U = A.LU(format='plu')
sage: all([A.is_mutable() for A in [P, L, U]])
True
Partial pivoting is based on the absolute values of entries of a column. Trac #12208 shows that the return value of the absolute value must be handled carefully. This tests that situation in the case of cylotomic fields.
sage: C = SymmetricGroup(5).character_table()
sage: C.base_ring()
Cyclotomic Field of order 1 and degree 1
sage: P, L, U = C.LU(pivot='partial')
sage: C == P*L*U
True
AUTHOR:
Return a numerical approximation of self as either a real or complex number with at least the requested number of bits or digits of precision.
INPUT:
OUTPUT: A matrix coerced to a real or complex field with prec bits of precision.
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: (eI).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()
[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]
Returns a factorization of self as a unitary matrix and an uppertriangular matrix.
INPUT:
OUTPUT:
If self is an \(m\times n\) matrix and full=True then this method returns a pair of matrices: \(Q\) is an \(m\times m\) unitary matrix (meaning its inverse is its conjugatetranspose) and \(R\) is an \(m\times n\) uppertriangular matrix with nonnegative entries on the diagonal. For a matrix of full rank this factorization is unique (due to the restriction to positive entries on the diagonal).
If full=False then \(Q\) has \(m\) rows and the columns form an orthonormal basis for the column space of self. So, in particular, the conjugatetranspose of \(Q\) times \(Q\) will be an identity matrix. The matrix \(R\) will still be uppertriangular but will also have full rank, in particular it will lack the zero rows present in a full factorization of a rankdeficient 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 GramSchmidt,” Algorithm 8.1 of [TREFETHENBAU].
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.687440062597?]
[ 0.4588314677411235? 0.4726901187474409? 0.05198346588033394? 0.717294125164660? 0.220962877263?]
[ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? 0.180872093737548? 0.1964114464561?]
[ 0.6882472016116853? 0.1890760474989764? 0.2044682991293135? 0.096630296654307? 0.662888631790?]
[ 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.682716185228386?]
[ 0 0 0 1.027626039419836? 3.61930014968662?]
[ 0 0 0 0 0.02455143080702?]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? 0.?e18 0.?e17 0.?e15 0.?e12]
[ 0.?e18 1.000000000000000? 0.?e16 0.?e15 0.?e12]
[ 0.?e17 0.?e16 1.000000000000000? 0.?e15 0.?e12]
[ 0.?e15 0.?e15 0.?e15 1.000000000000000? 0.?e12]
[ 0.?e12 0.?e12 0.?e12 0.?e12 1.000000000000?]
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.3786559533863032?  0.1952221495524667?*I 0.701244450214469?  0.364371165098660?*I]
[ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I 0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506763?  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.?e18  1.917028951268082?*I 5.385938482134133?  2.190890230020665?*I 0.2738612787525831?  2.190890230020665?*I]
[ 0 4.829596256417300? + 0.?e17*I 0.869637911123373?  5.864879483945125?*I 0.993871898426712?  0.3054085521207082?*I]
[ 0 0 12.00160760935814? + 0.?e16*I 0.2709533402297273? + 0.4420629644486323?*I]
[ 0 0 0 1.942963944258992? + 0.?e16*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e19*I 0.?e18 + 0.?e17*I 0.?e17 + 0.?e17*I 0.?e16 + 0.?e16*I]
[ 0.?e18 + 0.?e17*I 1.000000000000000? + 0.?e17*I 0.?e17 + 0.?e17*I 0.?e16 + 0.?e16*I]
[ 0.?e17 + 0.?e17*I 0.?e17 + 0.?e17*I 1.000000000000000? + 0.?e16*I 0.?e16 + 0.?e16*I]
[ 0.?e16 + 0.?e16*I 0.?e16 + 0.?e16*I 0.?e16 + 0.?e16*I 1.000000000000000? + 0.?e15*I]
sage: Q*R  A
[ 0.?e17 0.?e17 + 0.?e17*I 0.?e16 + 0.?e16*I 0.?e16 + 0.?e16*I]
[ 0.?e18 0.?e17 + 0.?e17*I 0.?e16 + 0.?e16*I 0.?e15 + 0.?e15*I]
[0.?e17 + 0.?e18*I 0.?e17 + 0.?e17*I 0.?e16 + 0.?e16*I 0.?e16 + 0.?e16*I]
[0.?e18 + 0.?e18*I 0.?e18 + 0.?e17*I 0.?e16 + 0.?e16*I 0.?e15 + 0.?e16*I]
A rankdeficient 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.?e18 0.?e18 0.?e18]
[ 0.?e18 1 0.?e18 0.?e18]
[ 0.?e18 0.?e18 1.000000000000000? 0.?e18]
[ 0.?e18 0.?e18 0.?e18 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.?e18]
[0.?e18 1]
Another rankdeficient 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.000000000000000? + 0.?e18*I 0.?e18 + 0.?e18*I]
[ 0.?e17 + 0.?e17*I 1.000000000000000? + 0.?e17*I]
sage: Q*RA
[0.?e18 + 0.?e18*I 0.?e18 + 0.?e18*I 0.?e17 + 0.?e17*I 0.?e18 + 0.?e18*I]
[0.?e18 + 0.?e18*I 0.?e18 + 0.?e18*I 0.?e18 + 0.?e17*I 0.?e18 + 0.?e18*I]
[0.?e18 + 0.?e18*I 0.?e17 + 0.?e18*I 0.?e17 + 0.?e17*I 0.?e18 + 0.?e18*I]
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]
TESTS:
Inexact rings are caught and CDF suggested.
sage: A = matrix(RealField(100), 2, range(4))
sage: A.QR()
Traceback (most recent call last):
...
NotImplementedError: QR decomposition is implemented over exact rings, try CDF for numerical results, not Real Field with 100 bits of precision
Without a fraction field, we cannot hope to run the algorithm.
sage: A = matrix(Integers(6), 2, range(4))
sage: A.QR()
Traceback (most recent call last):
...
ValueError: QR decomposition needs a fraction field of Ring of integers modulo 6
The biggest obstacle is making unit vectors, thus requiring square roots, though some small cases pass through.
sage: A = matrix(ZZ, 3, range(9))
sage: A.QR()
Traceback (most recent call last):
...
TypeError: QR decomposition unable to compute square roots in Rational Field
sage: A = matrix(ZZ, 2, range(4))
sage: Q, R = A.QR()
sage: Q
[0 1]
[1 0]
sage: R
[2 3]
[0 1]
REFERENCES:
[TREFETHENBAU]  (1, 2) Trefethen, Lloyd N., Bau, David, III “Numerical Linear Algebra” SIAM, Philadelphia, 1997. 
AUTHOR:
Returns the transpose of a matrix.
EXAMPLE:
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]
Returns the adjoint matrix of self (matrix of cofactors).
OUTPUT:
ALGORITHM:
Use PARI whenever the method self._adjoint is included to do so in an inheriting class. Otherwise, use a generic divisionfree algorithm to compute the characteristic polynomial and hence the adjoint.
The result is cached.
EXAMPLES:
sage: M = Matrix(ZZ,2,2,[5,2,3,4]) ; M
[5 2]
[3 4]
sage: N = M.adjoint() ; 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.adjoint() ; N
[ 41/10 1/28]
[33/13 5/3]
sage: M * N
[7363/1092 0]
[ 0 7363/1092]
AUTHORS:
Returns the current matrix as a sum of permutation matrices
According to the Birkhoffvon 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 nonbistochastic 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.
EXAMPLE:
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
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
Returns the characteristic polynomial of self, as a polynomial over the base ring.
ALGORITHM:
In the generic case of matrices over a ring (commutative and with unity), there is a divisionfree 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 divisionfree algorithm. If an algorithm is specified explicitly, if algorithm == "hessenberg", use the Hessenberg form, or otherwise use the generic divisionfree algorithm.
The result is cached.
INPUT:
var  a variable name (default: ‘x’)
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 multivariate 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, [[a1, 2], [a, a+1]])
sage: m.charpoly('Z')
Z^2  2*a*Z  2*a + 1
sage: m.charpoly('a')(m) == 0
True
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
TESTS:
sage: P.<a,b,c> = PolynomialRing(Rationals())
sage: u = MatrixSpace(P,3)([[0,0,a],[1,0,b],[0,1,c]])
sage: Q.<x> = PolynomialRing(P)
sage: u.charpoly('x')
x^3  c*x^2  b*x  a
A test case from trac ticket #6442. Prior to trac ticket #12292, the call to A.det() would attempt to use the cached charpoly, and crash if an empty dictionary was cached. We don’t cache dictionaries anymore, but this test should still pass:
sage: z = Zp(p=5)
sage: A = matrix(z, [ [3 + O(5^1), 4 + O(5^1), 4 + O(5^1)],
... [2*5^2 + O(5^3), 2 + O(5^1), 1 + O(5^1)],
... [5 + O(5^2), 1 + O(5^1), 1 + O(5^1)] ])
sage: A.charpoly(algorithm='hessenberg')
Traceback (most recent call last):
...
ValueError: negative valuation
sage: A.det()
3 + O(5)
The cached polynomial should be independent of the var argument (trac ticket #12292). We check (indirectly) that the second call uses the cached value by noting that its result is not cached:
sage: M = MatrixSpace(RR, 2)
sage: A = M(range(0, 2^2))
sage: type(A)
<type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
sage: A.charpoly('x')
x^2  3.00000000000000*x  2.00000000000000
sage: A.charpoly('y')
y^2  3.00000000000000*y  2.00000000000000
sage: A._cache['charpoly']
x^2  3.00000000000000*x  2.00000000000000
AUTHORS:
Returns the Cholesky decomposition of a symmetric or Hermitian matrix.
INPUT:
A square matrix that is real, symmetric and positive definite. Or a square matrix that is complex, Hermitian and positive definite. 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.
OUTPUT:
For a matrix \(A\) the routine returns a lower triangular matrix \(L\) such that,
where \(L^\ast\) is the conjugatetranspose in the complex case, and just the transpose in the real case. If the matrix fails to be positive definite (perhaps because it is not symmetric or Hermitian), then a ValueError results.
ALGORITHM:
Whether or not the matrix is positive definite is checked first in every case. This is accomplished with an indefinite factorization (see indefinite_factorization()) which caches its result. This algorithm is of an order \(n^3/3\). If the matrix is positive definite, this computation always succeeds, using just field operations. The transistion to a Cholesky decomposition “only” requires computing square roots of the positive (real) entries of the diagonal matrix produced in the indefinite factorization. Hence, there is no real penalty in the positive definite check (here, or prior to calling this routine), but a field extension with square roots may not be implemented in all reasonable cases.
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 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.change_ring(QQbar)).norm() < 10^10
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]
There are a variety of situations which will prevent the computation of a Cholesky decomposition.
The base ring must be exact. For numerical work, create a matrix with a base ring of RDF or CDF and use the cholesky() method for matrices of that type.
sage: F = RealField(100)
sage: A = matrix(F, [[1.0, 3.0], [3.0, 6.0]])
sage: A.cholesky()
Traceback (most recent call last):
...
TypeError: base ring of the matrix must be exact, not Real Field with 100 bits of precision
The base ring may not have a fraction field.
sage: A = matrix(Integers(6), [[2, 0], [0, 4]])
sage: A.cholesky()
Traceback (most recent call last):
...
ValueError: unable to check positive definiteness because
Unable to create the fraction field of Ring of integers modulo 6
The base field may not have elements that are comparable to zero.
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: unable to check positive definiteness because
cannot convert computations from Finite Field in a of size 5^4 into real numbers
The algebraic closure of the fraction field of the base ring may not be implemented.
sage: F = Integers(7)
sage: A = matrix(F, [[4, 0], [0, 3]])
sage: A.cholesky()
Traceback (most recent call last):
...
TypeError: base field needs an algebraic closure with square roots,
not Ring of integers modulo 7
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,
so cannot compute Cholesky decomposition
The matrix could be positive semidefinite, and thus lack a Cholesky decomposition.
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[:i,:i].determinant() for i in range(1,A.nrows()+1)]
[21, 27, 0, 0]
sage: A.cholesky()
Traceback (most recent call last):
...
ValueError: matrix is not positive definite,
so cannot compute Cholesky decomposition
Even in light of the above, you can sometimes get lucky and arrive at a situation where a particular matrix has a Cholesky decomposition when the general characteristics of the matrix suggest this routine would fail. In this example, the indefinite factorization produces a diagonal matrix whose elements from the finite field convert naturally to positive integers and are also perfect squares.
sage: F.<a> = FiniteField(5^3)
sage: A = matrix(F, [[ 4, 2*a^2 + 3, 4*a + 1],
... [ 2*a^2 + 3, 2*a + 2, 4*a^2 + 4*a + 4],
... [ 4*a + 1, 4*a^2 + 4*a + 4, a^2 + 4*a]])
sage: A.is_symmetric()
True
sage: L = A.cholesky()
sage: L*L.transpose() == A
True
TESTS:
This verifies that trac ticket #11274 is resolved.
sage: E = matrix(QQ, [[2, 1], [1, 1]])
sage: E.is_symmetric()
True
sage: E.eigenvalues()
[0.38...?, 2.61...?]
sage: E.det()
1
sage: E.cholesky()
[ 1.414213562373095? 0]
[0.7071067811865475? 0.7071067811865475?]
AUTHOR:
Return the free module over the base ring spanned by the columns of this matrix.
EXAMPLES:
sage: t = matrix(QQ, 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]
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]
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 nottotallyreal number field:
sage: K.<j> = NumberField(x^2+5)
sage: M = matrix(K, [[1+j,1], [0,2*j]])
sage: M.conjugate()
[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
Returns 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, [[2I, 3+4*I], [96*I, 5*I]])
sage: M.base_ring()
Symbolic Ring
sage: M.conjugate_transpose()
[ I + 2 6*I + 9]
[4*I + 3 5*I]
sage: P = matrix(CC, 3, 2, [0.950.63*I, 0.84+0.13*I, 0.94+0.23*I, 0.23+0.59*I, 0.520.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]
Conjugation does not make sense over rings not containing complex numbers.
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'
AUTHOR:
Rob Beezer (20101213)
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
INPUT:
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 \(k1\).
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: gen = polygen(QQ, 'z')
sage: p, S = A.cyclic_subspace(v, var=gen); p
z^3  9*z^2 + 24*z  16
sage: p.degree() == E.dimension()
True
The polynomial has coefficients that yield a nontrivial 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
TESTS:
A small case.
sage: A = matrix(QQ, 5, range(25))
sage: u = zero_vector(QQ, 5)
sage: A.cyclic_subspace(u)
Vector space of degree 5 and dimension 0 over Rational Field
Basis matrix:
[]
Various problem inputs. Notice the vector must have entries that coerce into the base ring of the matrix, and a polynomial ring generator must have a base ring that agrees with the base ring of the matrix.
sage: A = matrix(QQ, 4, range(16))
sage: v = vector(QQ, 4, range(4))
sage: A.cyclic_subspace('junk')
Traceback (most recent call last):
...
TypeError: first input should be a vector, not junk
sage: A.cyclic_subspace(v, var=sin(x))
Traceback (most recent call last):
...
TypeError: polynomial variable must be a string or polynomial ring generator, not sin(x)
sage: t = polygen(GF(7), 't')
sage: A.cyclic_subspace(v, var=t)
Traceback (most recent call last):
...
TypeError: polynomial generator must be over the same ring as the matrix entries
sage: A.cyclic_subspace(v, basis='garbage')
Traceback (most recent call last):
...
ValueError: basis format must be 'echelon' or 'iterates', not garbage
sage: B = matrix(QQ, 4, 5, range(20))
sage: B.cyclic_subspace(v)
Traceback (most recent call last):
...
TypeError: matrix must be square, not 4 x 5
sage: C = matrix(QQ, 5, 5, range(25))
sage: C.cyclic_subspace(v)
Traceback (most recent call last):
...
TypeError: vector must have degree equal to the size of the matrix, not 4
sage: D = matrix(RDF, 4, 4, range(16))
sage: D.cyclic_subspace(v)
Traceback (most recent call last):
...
TypeError: matrix entries must be from an exact ring, not Real Double Field
sage: E = matrix(Integers(6), 4, 4, range(16))
sage: E.cyclic_subspace(v)
Traceback (most recent call last):
...
TypeError: matrix entries must be from an exact field, not Ring of integers modulo 6
sage: F.<a> = GF(2^4)
sage: G = matrix(QQ, 4, range(16))
sage: w = vector(F, 4, [1, a, a^2, a^3])
sage: G.cyclic_subspace(w)
Traceback (most recent call last):
...
TypeError: unable to make vector entries compatible with matrix entries
AUTHOR:
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:
Note
If the base ring is not a field, the kernel algorithm is used.
OUTPUT:
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, 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)
]
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 restrict
with or without check.
kwds – Keywords that will be forwarded to decomposition().
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)
]
TESTS:
sage: t = matrix(QQ, 3, [3, 0, 2, 0, 2, 0, 0, 0, 0]);
sage: t.decomposition_of_subspace(v, check_restrict = False) == t.decomposition_of_subspace(v)
True
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
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.
EXAMPLE:
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()
5211/20000
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 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)
Synonym for self.determinant(...).
EXAMPLES:
sage: A = MatrixSpace(Integers(8),3)([1,7,3, 1,1,1, 3,4,5])
sage: A.det()
6
Returns the determinant of self.
ALGORITHM:
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 nonprime, 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 divisionfree algorithm (specified by "df"). When the base ring is an exact field, the default choice is "hessenberg", otherwise it is "df". Note that for matrices over most rings, more sophisticated algorithms can be used. (Type A.determinant? to see what is done for a specific matrix A.)
INPUT:
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
TESTS:
sage: A = matrix(5, 5, [next_prime(i^2) for i in range(25)])
sage: B = MatrixSpace(ZZ['x'], 5, 5)(A)
sage: A.det()  B.det()
0
We verify that trac ticket #5569 is resolved (otherwise the following would hang for hours):
sage: d = random_matrix(GF(next_prime(10^20)),50).det()
sage: d = random_matrix(Integers(10^50),50).det()
We verify that trac 7704 is resolved:
sage: matrix(ZZ, {(0,0):1,(1,1):2,(2,2):3,(3,3):4}).det()
24
sage: matrix(QQ, {(0,0):1,(1,1):2,(2,2):3,(3,3):4}).det()
24
We verify that trac 10063 is resolved:
sage: A = GF(2)['x,y,z']
sage: A.inject_variables()
Defining x, y, z
sage: R = A.quotient(x^2 + 1).quotient(y^2 + 1).quotient(z^2 + 1)
sage: R.inject_variables()
Defining xbarbarbar, ybarbarbar, zbarbarbar
sage: M = matrix([[1,1,1,1],[xbarbarbar,ybarbarbar,1,1],[0,1,zbarbarbar,1],[xbarbarbar,zbarbarbar,1,1]])
sage: M.determinant()
xbarbarbar*ybarbarbar*zbarbarbar + xbarbarbar*ybarbarbar + xbarbarbar*zbarbarbar + ybarbarbar*zbarbarbar + xbarbarbar + ybarbarbar + zbarbarbar + 1
Check that the determinant is computed from a cached charpoly properly:
sage: A = matrix(RR, [ [1, 0, 1/2],
... [0, 1, 0 ],
... [0, 0, 2 ] ])
sage: B = copy(A)
sage: _ = A.charpoly()
sage: A.determinant() == B.determinant()
True
AUTHORS:
 Unknown: No author specified in the file from 20090625
 Sebastian Pancratz (20090625): Use the divisionfree algorithm for charpoly
 Thierry Monteil (20101005): Bugfix for trac ticket #10063, so that the determinant is computed even for rings for which the is_field method is not implemented.
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()
[]
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:
OUTPUT:
The reduced row echelon form of self, as an immutable matrix. Note that self is not changed by this command. Use echelonize() to change self in place.
If the optional parameter transformation=True is specified, the output consists of a pair \((E,T)\) of matrices where \(E\) is the echelon form of self and \(T\) is the transformation matrix.
EXAMPLES:
sage: MS = MatrixSpace(GF(19),2,3)
sage: C = MS.matrix([1,2,3,4,5,6])
sage: C.rank()
2
sage: C.nullity()
0
sage: C.echelon_form()
[ 1 0 18]
[ 0 1 2]
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
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:
OUTPUT:
The matrix self is put into echelon form. Nothing is returned unless the keyword option transformation=True is specified, in which case the transformation matrix is returned.
EXAMPLES:
sage: a = matrix(QQ,3,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,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,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]
Echelon form is not defined over arbitrary rings:
sage: a = matrix(Integers(9),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
Return matrices D and P, where D is a diagonal matrix of eigenvalues and P is the corresponding matrix where the rows are corresponding eigenvectors (or zero vectors) so that P*self = D*P.
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: A = jordan_block(2,3)
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
TESTS:
For matrices with floating point entries, some platforms will return eigenvectors that are negatives of those returned by the majority of platforms. This test accounts for that possibility. Running this test independently, without adjusting the eigenvectors could indicate this situation on your hardware.
sage: A = matrix(QQ, 3, 3, range(9))
sage: em = A.change_ring(RDF).eigenmatrix_left()
sage: evalues = em[0]; evalues.dense_matrix().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 0.0]
[ 0.0 0.0 0.0]
sage: evectors = em[1];
sage: for i in range(3):
... scale = evectors[i,0].sign()
... evectors.rescale_row(i, scale)
sage: evectors
[ 0.440242867... 0.567868371... 0.695493875...]
[ 0.897878732... 0.278434036... 0.341010658...]
[ 0.408248290... 0.816496580... 0.408248290...]
Return matrices D and P, where D is a diagonal matrix of eigenvalues and P is the corresponding matrix where the columns are corresponding eigenvectors (or zero vectors) so that self*P = P*D.
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: A = jordan_block(2,3)
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
TESTS:
For matrices with floating point entries, some platforms will return eigenvectors that are negatives of those returned by the majority of platforms. This test accounts for that possibility. Running this test independently, without adjusting the eigenvectors could indicate this situation on your hardware.
sage: B = matrix(QQ, 3, 3, range(9))
sage: em = B.change_ring(RDF).eigenmatrix_right()
sage: evalues = em[0]; evalues.dense_matrix().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 0.0]
[ 0.0 0.0 0.0]
sage: evectors = em[1];
sage: for i in range(3):
... scale = evectors[0,i].sign()
... evectors.rescale_col(i, scale)
sage: evectors
[ 0.164763817... 0.799699663... 0.408248290...]
[ 0.505774475... 0.104205787... 0.816496580...]
[ 0.846785134... 0.591288087... 0.408248290...]
Compute the left eigenspaces of a matrix.
Note that eigenspaces_left() and left_eigenspaces() are identical methods. Here “left” refers to the eigenvectors being placed to the left of the matrix.
INPUT:
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)) < 1e10
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_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 9 0 0]
[ 0 9 0 1]
[ 0 10 9 2]
[ 0 9 0 2]
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 0 1 5]),
(3, Vector space of degree 4 and dimension 1 over Finite Field of size 11
User basis matrix:
[1 6 0 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 1 0 5*beta2 + 10])
]
This method is only applicable to exact matrices. The “eigenmatrix” routines for matrices with doubleprecision floatingpoint 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().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 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])
]
TESTS:
We make sure that trac ticket #13308 is fixed.
sage: M = ModularSymbols(Gamma1(23), sign=1)
sage: m = M.cuspidal_subspace().hecke_matrix(2)
sage: [j*m==i[0]*j for i in m.eigenspaces_left(format='all') for j in i[1].basis()] # long time (4s)
[True, True, True, True, True, True, True, True, True, True, True, True]
sage: B = matrix(QQ, 2, 3, range(6))
sage: B.eigenspaces_left()
Traceback (most recent call last):
...
TypeError: matrix must be square, not 2 x 3
sage: B = matrix(QQ, 4, 4, range(16))
sage: B.eigenspaces_left(format='junk')
Traceback (most recent call last):
...
ValueError: format keyword must be None, 'all' or 'galois', not junk
sage: B.eigenspaces_left(algebraic_multiplicity='garbage')
Traceback (most recent call last):
...
ValueError: algebraic_multiplicity keyword must be True or False
Compute the right eigenspaces of a matrix.
Note that eigenspaces_right() and right_eigenspaces() are identical methods. Here “right” refers to the eigenvectors being placed to the right of the matrix.
INPUT:
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)) < 1e10
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 doubleprecision floatingpoint 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().zero_at(1e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 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]
TESTS:
sage: B = matrix(QQ, 2, 3, range(6))
sage: B.eigenspaces_right()
Traceback (most recent call last):
...
TypeError: matrix must be square, not 2 x 3
sage: B = matrix(QQ, 4, 4, range(16))
sage: B.eigenspaces_right(format='junk')
Traceback (most recent call last):
...
ValueError: format keyword must be None, 'all' or 'galois', not junk
sage: B.eigenspaces_right(algebraic_multiplicity='garbage')
Traceback (most recent call last):
...
ValueError: algebraic_multiplicity keyword must be True or False
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(QQ, 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]
Compute the left eigenvectors of a matrix.
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)) < 1e10
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)]
Compute the right eigenvectors of a matrix.
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)) < 1e10
True
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
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: OE = EquationOrder(x^2  x + 2, 'w')
sage: w = OE.ring_generators()[0]
sage: m = Matrix([ [1, w],[w,7]])
sage: m.elementary_divisors()
[1, w + 9]
See also
Returns the elementwise product of two matrices of the same size (also known as the Hadamard product).
INPUT:
OUTPUT:
A matrix of the same size as self and right. The entry in location \((i,j)\) of the output is the product of the two entries in location \((i,j)\) of self and right (in that order).
The parent of the result is determined by Sage’s coercion model. If the base rings are identical, then the result is dense or sparse according to this property for the left operand. If the base rings must be adjusted for one, or both, matrices then the result will be sparse only if both operands are sparse. No subdivisions are present in the result.
If the type of the result is not to your liking, or the ring could be “tighter,” adjust the operands with change_ring(). Adjust sparse versus dense inputs with the methods sparse_matrix() and dense_matrix().
EXAMPLES:
sage: A = matrix(ZZ, 2, range(6))
sage: B = matrix(QQ, 2, [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: D = matrix(ZZ[x],2,[1+x^2,2,3,4x])
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
Noncommutative 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: operand must be a matrix, not an element of Ambient free module of rank 4 over the principal ideal domain Integer Ring
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: incompatible sizes for matrices from: 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, 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. Notice how if base rings are unequal, both operands must be sparse to get a sparse result.
sage: A = matrix(ZZ, 5, range(30), sparse=False)
sage: B = matrix(ZZ, 5, range(30), sparse=True)
sage: C = matrix(QQ, 5, range(30), sparse=True)
sage: A.elementwise_product(C).is_dense()
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
TESTS:
Implementation for dense and sparse matrices are different, this will provide a trivial test that they are working identically.
sage: A = random_matrix(ZZ, 10, x=1000, sparse=False)
sage: B = random_matrix(ZZ, 10, x=1000, sparse=False)
sage: C = A.sparse_matrix()
sage: D = B.sparse_matrix()
sage: E = A.elementwise_product(B)
sage: F = C.elementwise_product(D)
sage: E.is_dense() and F.is_sparse() and (E == F)
True
If the ring has zero divisors, the routines for setting entries of a sparse matrix should intercept zero results and not create an entry.
sage: R = Integers(6)
sage: A = matrix(R, 2, [3, 2, 0, 0], sparse=True)
sage: B = matrix(R, 2, [2, 3, 1, 0], sparse=True)
sage: C = A.elementwise_product(B)
sage: len(C.nonzero_positions()) == 0
True
AUTHOR:
Calculate the exponential of this matrix X, which is the matrix
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*(e^sqrt(33)*(sqrt(33)  11)  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*(e^sqrt(33)*(sqrt(33) + 11)  sqrt(33) + 11)*e^(1/2*sqrt(33) + 5/2)]
sage: type(a.exp())
<type 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>
sage: a=matrix([[1/2,2/3],[3/4,4/5]])
sage: a.exp()
[1/418*(e^(1/10*sqrt(209))*(3*sqrt(209)  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*(e^(1/10*sqrt(209))*(3*sqrt(209) + 209)  3*sqrt(209) + 209)*e^(1/20*sqrt(209) + 13/20)]
sage: a=matrix(RR,[[1,pi.n()],[1e2,1e2]])
sage: a.exp()
[ 1/11882424341266*(e^(3/1275529100*sqrt(227345670387496707609))*(11*sqrt(227345670387496707609) + 5941212170633)  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*(e^(3/1275529100*sqrt(227345670387496707609))*(11*sqrt(227345670387496707609)  5941212170633)  11*sqrt(227345670387496707609)  5941212170633)*e^(3/2551058200*sqrt(227345670387496707609) + 101/200)]
sage: a.change_ring(RDF).exp()
[42748127.3153 7368259.24416]
[234538976.138 40426191.4516]
Returns the echelon form of self augmented with an identity matrix.
INPUT:
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 rowreduce this \(m\times(n+m)\) matrix. This matrix is returned as an immutable matrix.
If subdivide is True then the returned matrix has a single division among the columns and a single division among the rows. The column subdivision has \(n\) columns to the left and \(m\) columns to the right. The row division separates the nonzero rows from the zero rows, when restricted to the first \(n\) columns.
For a nonsingular matrix the final \(m\) columns of the extended echelon form are the inverse of self. For a matrix of any size, the final \(m\) columns provide a matrix that transforms self to echelon form when it multiplies self from the left. When the base ring is a field, the uniqueness of reduced rowechelon form implies that this transformation matrix can be taken as the coefficients giving a canonical set of linear combinations of the rows of self that yield reduced rowechelon form.
When subdivided as described above, and again over a field, the parts of the subdivision in the upperleft corner and lowerright 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 atttempt 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 00 4 6]
[0 1 04 2 2]
[0 0 15 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 flag include_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]
TESTS:
The subdivide keyword is checked.
sage: A = matrix(QQ, 2, range(4))
sage: A.extended_echelon_form(subdivide='junk')
Traceback (most recent call last):
...
TypeError: subdivide must be True or False, not junk
AUTHOR:
Return the factorization of the characteristic polynomial of self.
INPUT:
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 elements in this matrix satisfying the constraints in the function \(f\). The function is evaluated on each element of the matrix .
INPUT:
OUTPUT: If indices is not specified, return a matrix with 1 where \(f\) is satisfied and 0 where it is not. If indices 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
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 00 0]
[0 1 00 0]
[+]
[0 0 10 0]
[0 0 01 0]
[0 0 00 1]
Performs GramSchmidt orthogonalization on the rows of the matrix, returning a new matrix and a matrix accomplishing the transformation.
INPUT:
OUTPUT:
A pair of matrices, G and M such that if A represents self, where the parenthetical properties occur when orthonormal = True:
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 GramSchmidt 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 doubleprecision floating point numbers, RDF and CDF. 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 and RDF.
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()).zero_at(10^12)
[1.0 0.0 0.0]
[0.0 1.0 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: (AM*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(1e14)
[0.0 0.0 0.0]
[0.0 0.0 0.0]
[0.0 0.0 0.0]
sage: (G*G.transpose()).zero_at(1e14)
[1.0 0.0 0.0]
[0.0 1.0 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 the orthonormal 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*GA
[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 trac ticket #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 conjugatetranspose 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*GA).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 of A.
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*GA
[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 GramSchmidt 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*GA
[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
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 53bits 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
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
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
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).
Return the image of the homomorphism on rows defined by this matrix.
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
Decomposes a symmetric or Hermitian matrix into a lower triangular matrix and a diagonal matrix.
INPUT:
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
If \(A\) is Hermitian matrix, then the transpose of \(L\) should be replaced by the conjugatetranspose of \(L\).
If any leading principal submatrix (a square submatrix in the upperleft corner) is singular then this method will fail with a ValueError.
ALGORITHM:
The algorithm employed only uses field operations, but the compuation 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.
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
AUTHOR:
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: L.<w> = NumberField(x^2  x + 2)
sage: OL = L.ring_of_integers()
sage: A = matrix(L, 2, [1, w/2])
sage: A.integer_kernel(OL)
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]
Returns the inverse of self, without changing self.
Note that one can use the Python inverse operator to obtain the inverse as well.
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]
sage: m.I
[ 2 1]
[ 3/2 1/2]
TESTS:
sage: matrix().inverse()
[]
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.
INPUT:
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
Determines if the matrix is similar to a diagonal matrix.
INPUT:
OUTPUT:
If self is the matrix \(A\), then it is diagonalizable if there is an invertible matrix \(S\) and a diagonal matrix \(D\) such that
This routine returns True if self is diagonalizable. The diagonal entries of the matrix \(D\) are the eigenvalues of \(A\). It may be necessary to “increase” 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 \(S\) and \(D\) use the jordan_form() method with the transformation=True keyword.
ALGORITHM:
For each eigenvalue, this routine checks that the algebraic multiplicity (number of occurences 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 \(S\).
EXAMPLES:
A matrix that is diagonalizable over the rationals, as evidenced by its Jordan form.
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.jordan_form(subdivide=False)
[ 2 0 0 0 0]
[ 0 3 0 0 0]
[ 0 0 3 0 0]
[ 0 0 0 1 0]
[ 0 0 0 0 1]
sage: A.is_diagonalizable()
True
A matrix that is not diagonalizable over the rationals, as evidenced by 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.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]
sage: A.is_diagonalizable()
False
If any eigenvalue of a matrix is outside the base ring, then this routine raises an error. However, the ring can be “expanded” 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()
Traceback (most recent call last):
...
RuntimeError: an eigenvalue of the matrix is not contained in 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 expand 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.jordan_form()
[ 4 1 0 0]
[ 0 4 0 0]
[+]
[ 0 02*b + 1 1]
[ 0 0 0 2*b + 1]
sage: A.is_diagonalizable()
False
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.jordan_form(subdivide=False)
[ 4 0 0 0]
[ 0 2 0 0]
[ 0 0 c + 3 0]
[ 0 0 0 c + 3]
sage: A.is_diagonalizable()
True
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()
False
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, not Integer Ring
sage: A = matrix(RDF, 4, range(16))
sage: A.is_diagonalizable()
Traceback (most recent call last):
...
ValueError: base field must be exact, not Real Double Field
AUTHOR:
Returns True if the matrix commutes with its conjugatetranspose.
OUTPUT:
True if the matrix is square and commutes with its conjugatetranspose, and False 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 or CC. For matrices with approximate entries, the rings of doubleprecision floatingpoint numbers, RDF and CDF, are a better choice since the sage.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, range(25)) + I*matrix(QQ, 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
Skewsymmetric matrices are normal.
sage: A = matrix(QQ, 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 nontrivial.
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 nonreal diagonal entry.
sage: A = matrix(QQbar, [[ 2, 2I, 1+4*I],
... [ 2+I, 3+I, 26*I],
... [14*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
AUTHOR:
Return True if this matrix is the identity matrix.
EXAMPLES:
sage: m = matrix(QQ,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
Determines if a real or symmetric matrix is positive definite.
A square matrix \(A\) is postive definite if it is symmetric with real entries or Hermitan with complex entries, and for every nonzero vector \(\vec{x}\)
Here \(\vec{x}^\ast\) is the conjugatetranspose, which can be simplified to just the transpose in the real case.
ALGORITHM:
A matrix is positive definite if and only if the diagonal entries from the indefinite factorization are all positive (see indefinite_factorization()). So this algorithm is of order n^3/3 and may be applied to matrices with elements of any ring that has a fraction field contained within the reals or complexes.
INPUT:
Any square matrix.
OUTPUT:
This routine will return True if the matrix is square, symmetric or Hermitian, and meets the condition above for the quadratic form.
The base ring for the elements of the matrix needs to have a fraction field implemented and the computations that result from the indefinite factorization must be convertable to real numbers that are comparable to zero.
EXAMPLES:
A real symmetric matrix that is positive definite, as evidenced by the positive entries for the diagonal matrix of the indefinite factorization and the postive determinants of the 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: _, d = A.indefinite_factorization(algorithm='symmetric')
sage: d
(4, 9, 4, 1)
sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
[4, 36, 144, 144]
A real symmetric matrix which is not positive definite, along with a vector that makes the 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: _, d = A.indefinite_factorization(algorithm='symmetric')
sage: d
(3, 1, 5, 2, 1)
sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
[3, 3, 15, 30, 30]
sage: u = vector(QQ, [2, 2, 0, 1, 0])
sage: u.row()*A*u
(3)
A real symmetric matrix with a singular leading principal submatrix, that is therefore not positive definite. The vector u makes the 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: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
[21, 27, 0, 75]
sage: u = vector(QQ, [1,1,3,0])
sage: u.row()*A*u
(0)
An Hermitian matrix that is positive definite.
sage: C.<I> = NumberField(x^2 + 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_positive_definite()
True
sage: _, d = A.indefinite_factorization(algorithm='hermitian')
sage: d
(23, 576/23, 89885/144, 142130/17977)
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. The vector u makes the 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: _, d = B.indefinite_factorization(algorithm='hermitian')
sage: d
(2, 2, 3)
sage: [B[:i,:i].determinant() for i in range(1,B.nrows()+1)]
[2, 4, 12]
sage: u = vector(C, [5 + 10*I, 4  3*I, 0])
sage: u.row().conjugate()*B*u
(50)
A positive definite matrix over an algebraically closed field.
sage: A = matrix(QQbar, [[ 2, 4 + 2*I, 6  4*I],
... [ 2*I + 4, 11, 10  12*I],
... [ 4*I + 6, 10 + 12*I, 37]])
sage: A.is_positive_definite()
True
sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
[2, 2, 6]
TESTS:
If the base ring lacks a conjugate method, it will be assumed to not be Hermitian and thus symmetric. If the base ring does not make sense as a subfield of the reals, then this routine will fail since comparison to zero is meaningless.
sage: F.<a> = FiniteField(5^3)
sage: a.conjugate()
Traceback (most recent call last):
...
AttributeError: 'sage.rings.finite_rings.element_givaro.FiniteField_givaroElement'
object has no attribute 'conjugate'
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_positive_definite()
Traceback (most recent call last):
...
TypeError: cannot convert computations from
Finite Field in a of size 5^3 into real numbers
AUTHOR:
Return True if this matrix is a scalar matrix.
INPUT
OUTPUT
EXAMPLES:
sage: m = matrix(QQ,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
Returns True if self and other are similar, i.e. related by a changeofbasis matrix.
INPUT:
OUTPUT:
Two matrices, \(A\) and \(B\) are similar if there is an invertible matrix \(S\) such that \(A=S^{1}BS\). \(S\) can be interpreted as a changeofbasis matrix if \(A\) and \(B\) are viewed as matrix representations of the same linear transformation.
When transformation=False this method will return True if such a matrix \(S\) exists, otherwise it will return False. When transformation=True the method returns a pair. The first part of the pair is True or False depending on if the matrices are similar and the second part is the changeofbasis matrix, or None should it not exist.
When the transformation matrix is requested, it will satisfy self = S.inverse()*other*S.
If the base rings for any of the matrices is the integers, the rationals, or the field of algebraic numbers (QQbar), then the matrices are converted to have QQbar as their base ring prior to checking the equality of the base rings.
It is possible for this routine to fail over most fields, even when the matrices are similar. However, since the field of algebraic numbers is algebraically closed, the routine will always produce a result for matrices with rational entries.
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 changeofbasis 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.0000000000000? + 0.?e13*I 0.?e13 + 0.?e13*I 0.?e13 + 0.?e13*I]
[0.6666666666667? + 0.?e13*I 0.16666666666667? + 0.?e14*I 0.8333333333334? + 0.?e13*I]
[ 0.6666666666667? + 0.?e13*I 0.?e13 + 0.?e13*I 0.333333333334? + 0.?e13*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 Jordan form, as provided by jordan_form(). The matrices below have identical eigenvalues (as evidenced by equal characteristic polynomials), but slightly different Jordan 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.jordan_form()
[3 0 0]
[+]
[ 0 4 1]
[ 0 0 4]
sage: B.jordan_form()
[3 0 0]
[++]
[ 0 4 0]
[++]
[ 0 0 4]
sage: A.is_similar(B)
False
Obtaining the Jordan form requires computing the eigenvalues of the matrix, which may not lie in the field used for entries of the matrix. So the routine first checks the characteristic polynomials  if they are unequal, then the matrices cannot be similar. However, when the characteristic polynomials are equal, we must examine the Jordan form. In this case, the method may fail, 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, but the routine fails to compute a result.
sage: F.<a> = FiniteField(7^2)
sage: C = matrix(F,[[ a + 2, 5*a + 4],
... [6*a + 6, 6*a + 4]])
sage: S = matrix(ZZ, [[0, 1],
... [1, 0]])
sage: D = S.inverse()*C*S
sage: C.is_similar(D)
Traceback (most recent call last):
...
ValueError: unable to compute Jordan canonical form for a matrix
sage: C.jordan_form()
Traceback (most recent call last):
...
RuntimeError: Some eigenvalue does not exist in Finite Field in a of size 7^2.
Inexact rings and fields are also not supported.
sage: A = matrix(CDF, 2, 2, range(4))
sage: B = copy(A)
sage: A.is_similar(B)
Traceback (most recent call last):
...
ValueError: unable to compute Jordan canonical form for a matrix
Rectangular matrices and mismatched sizes return quickly.
sage: A = matrix(3, 2, range(6))
sage: B = copy(A)
sage: A.is_similar(B)
False
sage: A = matrix(2, 2, range(4))
sage: B = matrix(3, 3, range(9))
sage: A.is_similar(B, transformation=True)
(False, None)
If the fraction fields of the entries are unequal, it is an error, except in the case when the rationals gets promoted to the algebraic numbers.
sage: A = matrix(ZZ, 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: matrices need to have entries with identical fraction fields, not Algebraic Field and Finite Field of size 2
sage: A = matrix(ZZ, 2, 2, range(4))
sage: B = matrix(QQbar, 2, 2, range(4))
sage: A.is_similar(B)
True
Inputs are checked.
sage: A = matrix(ZZ, 2, 2, range(4))
sage: A.is_similar('garbage')
Traceback (most recent call last):
...
TypeError: similarity requires a matrix as an argument, not garbage
sage: B = copy(A)
sage: A.is_similar(B, transformation='junk')
Traceback (most recent call last):
...
ValueError: transformation keyword must be True or False, not junk
Returns 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 conjugatetranspose is its inverse, and False otherwise. In other words, a matrix is orthogonal or unitary if the product of its conjugatetranspose times the matrix is the identity matrix.
For numerical matrices a specialized routine available over RDF and CDF 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))*(1i), (1/sqrt(55))*(2+2*I), (1/sqrt(22))*(3+I)],
... [ (1/sqrt(5))*I, (1/sqrt(55))*(35*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
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 \(AxI\), 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\).
INPUT:
NOTES:
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()
[20 00]
[++]
[01 10]
[00 10]
[++]
[00 01]
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,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 (trac ticket #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 0zeta3  1]
If you need the transformation matrix as well as the Jordan form of self, then pass the option transformation=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
[200 0]
[++]
[010 0]
[++]
[004 1]
[000 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.
TESTS:
sage: c = matrix(ZZ, 3, [1]*9); c
[1 1 1]
[1 1 1]
[1 1 1]
sage: c.jordan_form(subdivide=False)
[3 0 0]
[0 0 0]
[0 0 0]
sage: evals = [(i,i) for i in range(1,6)]
sage: n = sum(range(1,6))
sage: jf = block_diagonal_matrix([jordan_block(ev,size) for ev,size in evals])
sage: p = random_matrix(ZZ,n,n)
sage: while p.rank() != n: p = random_matrix(ZZ,n,n)
sage: m = p * jf * ~p
sage: mjf, mp = m.jordan_form(transformation=True)
sage: mjf == jf
True
sage: m = diagonal_matrix([1,1,0,0])
sage: jf,P = m.jordan_form(transformation=True)
sage: jf == ~P*m*P
True
We verify that the bug from trac ticket #6942 is fixed:
sage: M = Matrix(GF(2),[[1,0,1,0,0,0,1],[1,0,0,1,1,1,0],[1,1,0,1,1,1,1],[1,1,1,0,1,1,1],[1,1,1,0,0,1,0],[1,1,1,0,1,0,0],[1,1,1,1,1,1,0]])
sage: J, T = M.jordan_form(transformation=True)
sage: J
[1 10 00 00]
[0 10 00 00]
[+++]
[0 01 10 00]
[0 00 10 00]
[+++]
[0 00 01 10]
[0 00 00 10]
[+++]
[0 00 00 01]
sage: M * T == T * J
True
sage: T.rank()
7
sage: M.rank()
7
We verify that the bug from trac ticket #6932 is fixed:
sage: M=Matrix(1,1,[1])
sage: M.jordan_form(transformation=True)
([1], [1])
We now go through three \(10 \times 10\) matrices to exhibit cases where there are multiple blocks of the same size:
sage: A = matrix(QQ, [[15, 37/3, 16, 104/3, 29, 7/3, 0, 2/3, 29/3, 1/3], [2, 9, 1, 6, 6, 0, 0, 0, 2, 0], [24, 74/3, 41, 208/3, 58, 23/3, 0, 4/3, 58/3, 2/3], [6, 19, 3, 21, 19, 0, 0, 0, 6, 0], [2, 6, 3, 6, 3, 1, 0, 0, 2, 0], [96, 296/3, 176, 832/3, 232, 101/3, 0, 16/3, 232/3, 8/3], [4, 2/3, 21, 16/3, 4, 14/3, 3, 1/3, 4/3, 25/3], [20, 26/3, 66, 199/3, 42, 41/3, 0, 13/3, 55/3, 2/3], [18, 57, 9, 54, 57, 0, 0, 0, 15, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 3]]); A
[ 15 37/3 16 104/3 29 7/3 0 2/3 29/3 1/3]
[ 2 9 1 6 6 0 0 0 2 0]
[ 24 74/3 41 208/3 58 23/3 0 4/3 58/3 2/3]
[ 6 19 3 21 19 0 0 0 6 0]
[ 2 6 3 6 3 1 0 0 2 0]
[ 96 296/3 176 832/3 232 101/3 0 16/3 232/3 8/3]
[ 4 2/3 21 16/3 4 14/3 3 1/3 4/3 25/3]
[ 20 26/3 66 199/3 42 41/3 0 13/3 55/3 2/3]
[ 18 57 9 54 57 0 0 0 15 0]
[ 0 0 0 0 0 0 0 0 0 3]
sage: J, T = A.jordan_form(transformation=True); J
[3 1 00 0 00 0 00]
[0 3 10 0 00 0 00]
[0 0 30 0 00 0 00]
[+++]
[0 0 03 1 00 0 00]
[0 0 00 3 10 0 00]
[0 0 00 0 30 0 00]
[+++]
[0 0 00 0 03 1 00]
[0 0 00 0 00 3 10]
[0 0 00 0 00 0 30]
[+++]
[0 0 00 0 00 0 03]
sage: T * J * T**(1) == A
True
sage: T.rank()
10
sage: A = matrix(QQ, [[15, 37/3, 16, 14/3, 29, 7/3, 0, 2/3, 1/3, 44/3], [2, 9, 1, 0, 6, 0, 0, 0, 0, 3], [24, 74/3, 41, 28/3, 58, 23/3, 0, 4/3, 2/3, 88/3], [6, 19, 3, 3, 19, 0, 0, 0, 0, 9], [2, 6, 3, 0, 3, 1, 0, 0, 0, 3], [96, 296/3, 176, 112/3, 232, 101/3, 0, 16/3, 8/3, 352/3], [4, 2/3, 21, 16/3, 4, 14/3, 3, 1/3, 4/3, 25/3], [20, 26/3, 66, 28/3, 42, 41/3, 0, 13/3, 2/3, 82/3], [18, 57, 9, 0, 57, 0, 0, 0, 3, 28], [0, 0, 0, 0, 0, 0, 0, 0, 0, 3]]); A
[ 15 37/3 16 14/3 29 7/3 0 2/3 1/3 44/3]
[ 2 9 1 0 6 0 0 0 0 3]
[ 24 74/3 41 28/3 58 23/3 0 4/3 2/3 88/3]
[ 6 19 3 3 19 0 0 0 0 9]
[ 2 6 3 0 3 1 0 0 0 3]
[ 96 296/3 176 112/3 232 101/3 0 16/3 8/3 352/3]
[ 4 2/3 21 16/3 4 14/3 3 1/3 4/3 25/3]
[ 20 26/3 66 28/3 42 41/3 0 13/3 2/3 82/3]
[ 18 57 9 0 57 0 0 0 3 28]
[ 0 0 0 0 0 0 0 0 0 3]
sage: J, T = A.jordan_form(transformation=True); J
[3 1 00 0 00 00 0]
[0 3 10 0 00 00 0]
[0 0 30 0 00 00 0]
[+++]
[0 0 03 1 00 00 0]
[0 0 00 3 10 00 0]
[0 0 00 0 30 00 0]
[+++]
[0 0 00 0 03 10 0]
[0 0 00 0 00 30 0]
[+++]
[0 0 00 0 00 03 1]
[0 0 00 0 00 00 3]
sage: T * J * T**(1) == A
True
sage: T.rank()
10
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]]); A
[ 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: J, T = A.jordan_form(transformation=True); J
[3 1 00 00 00 00]
[0 3 10 00 00 00]
[0 0 30 00 00 00]
[++++]
[0 0 03 10 00 00]
[0 0 00 30 00 00]
[++++]
[0 0 00 03 10 00]
[0 0 00 00 30 00]
[++++]
[0 0 00 00 03 10]
[0 0 00 00 00 30]
[++++]
[0 0 00 00 00 03]
sage: T * J * T**(1) == A
True
sage: T.rank()
10
Verify that we smoothly move to QQ from ZZ (trac ticket #12693), i.e. we work in the vector space over the field:
sage: M = matrix(((2,2,2),(0,0,0),(2,2,2)))
sage: J, P = M.jordan_form(transformation=True)
sage: J; P
[0 10]
[0 00]
[+]
[0 00]
[ 2 1 0]
[ 0 0 1]
[2 0 1]
sage: J  ~P * M * P
[0 0 0]
[0 0 0]
[0 0 0]
sage: parent(M)
Full MatrixSpace of 3 by 3 dense matrices over Integer Ring
sage: parent(J) == parent(P) == MatrixSpace(QQ, 3)
True
sage: M.jordan_form(transformation=True) == (M/1).jordan_form(transformation=True)
True
By providing eigenvalues ourselves, we can compute the Jordan form even lacking a polynomial factorization algorithm.
sage: Qx = PolynomialRing(QQ, 'x11, x12, x13, x21, x22, x23, x31, x32, x33')
sage: x11, x12, x13, x21, x22, x23, x31, x32, x33 = Qx.gens()
sage: M = matrix(Qx, [[0, 0, x31], [0, 0, x21], [0, 0, 0]]) # This is a nilpotent matrix.
sage: M.jordan_form(eigenvalues=[(0, 3)])
[0 10]
[0 00]
[+]
[0 00]
sage: M.jordan_form(eigenvalues=[(0, 2)])
Traceback (most recent call last):
...
ValueError: The provided list of eigenvalues is not correct.
sage: M.jordan_form(transformation=True, eigenvalues=[(0, 3)])
(
[0 10]
[0 00] [x31 0 1]
[+] [x21 0 0]
[0 00], [ 0 1 0]
)
TESTS:
The base ring for the matrix needs to have a fraction field and it needs to be implemented.
sage: A = matrix(Integers(6), 2, 2, range(4))
sage: A.jordan_form()
Traceback (most recent call last):
...
ValueError: Matrix entries must be from a field, not Ring of integers modulo 6
Returns the left kernel of this matrix, as a vector space or free module. This is the set of vectors x such that x*self = 0.
Note
For the right kernel, use right_kernel(). The method kernel() is exactly equal to left_kernel().
INPUT:
OUTPUT:
A vector space or free module whose degree equals the number of rows in self and contains all the vectors x such that x*self = 0.
If self has 0 rows, the kernel has dimension 0, while if self 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 of self.
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
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:
OUTPUT:
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 = (t39)*(t^22)
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 = t39
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 trac ticket #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
Return matrices D and P, where D is a diagonal matrix of eigenvalues and P is the corresponding matrix where the rows are corresponding eigenvectors (or zero vectors) so that P*self = D*P.
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: A = jordan_block(2,3)
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
TESTS:
For matrices with floating point entries, some platforms will return eigenvectors that are negatives of those returned by the majority of platforms. This test accounts for that possibility. Running this test independently, without adjusting the eigenvectors could indicate this situation on your hardware.
sage: A = matrix(QQ, 3, 3, range(9))
sage: em = A.change_ring(RDF).eigenmatrix_left()
sage: evalues = em[0]; evalues.dense_matrix().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 0.0]
[ 0.0 0.0 0.0]
sage: evectors = em[1];
sage: for i in range(3):
... scale = evectors[i,0].sign()
... evectors.rescale_row(i, scale)
sage: evectors
[ 0.440242867... 0.567868371... 0.695493875...]
[ 0.897878732... 0.278434036... 0.341010658...]
[ 0.408248290... 0.816496580... 0.408248290...]
Compute the left eigenspaces of a matrix.
Note that eigenspaces_left() and left_eigenspaces() are identical methods. Here “left” refers to the eigenvectors being placed to the left of the matrix.
INPUT:
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)) < 1e10
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_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 9 0 0]
[ 0 9 0 1]
[ 0 10 9 2]
[ 0 9 0 2]
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 0 1 5]),
(3, Vector space of degree 4 and dimension 1 over Finite Field of size 11
User basis matrix:
[1 6 0 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 1 0 5*beta2 + 10])
]
This method is only applicable to exact matrices. The “eigenmatrix” routines for matrices with doubleprecision floatingpoint 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().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 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])
]
TESTS:
We make sure that trac ticket #13308 is fixed.
sage: M = ModularSymbols(Gamma1(23), sign=1)
sage: m = M.cuspidal_subspace().hecke_matrix(2)
sage: [j*m==i[0]*j for i in m.eigenspaces_left(format='all') for j in i[1].basis()] # long time (4s)
[True, True, True, True, True, True, True, True, True, True, True, True]
sage: B = matrix(QQ, 2, 3, range(6))
sage: B.eigenspaces_left()
Traceback (most recent call last):
...
TypeError: matrix must be square, not 2 x 3
sage: B = matrix(QQ, 4, 4, range(16))
sage: B.eigenspaces_left(format='junk')
Traceback (most recent call last):
...
ValueError: format keyword must be None, 'all' or 'galois', not junk
sage: B.eigenspaces_left(algebraic_multiplicity='garbage')
Traceback (most recent call last):
...
ValueError: algebraic_multiplicity keyword must be True or False
Compute the left eigenvectors of a matrix.
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)) < 1e10
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)]
Returns the left kernel of this matrix, as a vector space or free module. This is the set of vectors x such that x*self = 0.
Note
For the right kernel, use right_kernel(). The method kernel() is exactly equal to left_kernel().
INPUT:
OUTPUT:
A vector space or free module whose degree equals the number of rows in self and contains all the vectors x such that x*self = 0.
If self has 0 rows, the kernel has dimension 0, while if self 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 of self.
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
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
Return the requested matrix window.
EXAMPLES:
sage: A = matrix(QQ, 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
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:
OUTPUT:
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, 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)]
This is a synonym for self.minpoly
EXAMPLES:
sage: a = matrix(QQ, 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
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:
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]
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
We can take the minimal polynomial of symbolic matrices:
sage: t = var('t')
sage: m = matrix(2,[1,2,4,t])
sage: m.minimal_polynomial()
x^2 + (t  1)*x + t  8
Return a numerical approximation of self as either a real or complex number with at least the requested number of bits or digits of precision.
INPUT:
OUTPUT: A matrix coerced to a real or complex field with prec bits of precision.
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: (eI).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()
[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]
Return the pnorm of this matrix, where \(p\) can be 1, 2, \(\inf\), or the Frobenius norm.
INPUT:
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()
15.0
Norms of numerical matrices over highprecision 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,1I,1,2,0])
sage: A.norm('frob')
5.65685424949
sage: A.norm(2)
5.47068444321
sage: A.norm(1)
6.0
sage: A.norm(Infinity)
8.41421356237
sage: a = matrix([[],[],[],[]])
sage: a.norm()
0.0
sage: a.norm(Infinity) == a.norm(1)
True
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
Return a numerical approximation of self as either a real or complex number with at least the requested number of bits or digits of precision.
INPUT:
OUTPUT: A matrix coerced to a real or complex field with prec bits of precision.
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: (eI).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()
[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]
Calculate and return the permanent of the \(m \times n\) matrix self using Ryser’s algorithm.
Let \(A = (a_{i,j})\) be an \(m \times n\) matrix over any commutative ring, with \(m \le n\). The permanent of \(A\) is
where the summation extends over all onetoone 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\).
INPUT:
OUTPUT:
permanent of the matrix \(A\)
ALGORITHM:
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.
EXAMPLES:
sage: M = MatrixSpace(ZZ,4,4)
sage: A = M([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1])
sage: A.permanent()
24
sage: M = MatrixSpace(QQ,3,6)
sage: A = M([1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1])
sage: A.permanent()
36
sage: M = MatrixSpace(RR,3,6)
sage: A = M([1.0,1.0,1.0,1.0,0,0,0,1.0,1.0,1.0,1.0,0,0,0,1.0,1.0,1.0,1.0])
sage: A.permanent()
36.0000000000000
See Sloane’s sequence OEIS A079908(3) = 36, “The Dancing School Problems”
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: M = MatrixSpace(ZZ,4,5)
sage: A = M([1,1,0,1,1,0,1,1,1,1,1,0,1,0,1,1,1,0,1,0])
sage: A.permanent()
32
See Minc: Permanents, Example 2.1, p. 5.
sage: M = MatrixSpace(QQ,2,2)
sage: A = M([1/5,2/7,3/2,4/5])
sage: A.permanent()
103/175
sage: R.<a> = PolynomialRing(ZZ)
sage: A = MatrixSpace(R,2)([[a,1], [a,a+1]])
sage: A.permanent()
a^2 + 2*a
sage: R.<x,y> = PolynomialRing(ZZ,2)
sage: A = MatrixSpace(R,2)([x, y, x^2, y^2])
sage: A.permanent()
x^2*y + x*y^2
AUTHORS:
Return the permanental \(k\)minor of an \(m \times n\) matrix.
This is the sum of the permanents of all possible \(k\) by \(k\) submatrices of \(A\).
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.
Note that the permanental \(m\)minor equals \(\mathrm{per}(A)\) if \(m = n\).
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.
INPUT:
OUTPUT:
The permanental \(k\)minor of the matrix self.
EXAMPLES:
sage: M = MatrixSpace(ZZ,4,4)
sage: A = M([1,0,1,0,1,0,1,0,1,0,10,10,1,0,1,1])
sage: A.permanental_minor(2)
114
sage: M = MatrixSpace(ZZ,3,6)
sage: A = M([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
sage: A.permanental_minor(5)
0
For \(C\) the “complement” of \(A\):
sage: M = MatrixSpace(ZZ,3,6)
sage: C = M([0,0,0,0,1,1,1,0,0,0,0,1,1,1,0,0,0,0])
sage: m, n = 3, 6
sage: sum([(1)^k * C.permanental_minor(k)*factorial(nk)/factorial(nm) for k in range(m+1)])
36
See Theorem 7.2.1 of Brualdi and Ryser: Combinatorial Matrix Theory: per(A)
AUTHORS:
Return the Pfaffian of self, assuming that self is an alternating matrix.
INPUT:
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 oneline 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
The Pfaffian of \(A\) is commonly denoted by \(\mathrm{Pf}(A)\). It is wellknown 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 [Kn95], [DW95] and [Rote2001], just to name three sources, for further properties of Pfaffians.
ALGORITHM:
The current implementation uses the definition given above. It checks alternatingness of the matrix self only if check is True (this is important because even if self is alternating, a nondiscrete base ring might prevent Sage from being able to check this).
REFERENCES:
[Kn95]  Donald E. Knuth, Overlapping Pfaffians, Arxiv math/9503234v1. 
[Rote2001]  (1, 2) Gunter Rote, DivisionFree Algorithms for the Determinant and the Pfaffian: Algebraic and Combinatorial Approaches, H. Alt (Ed.): Computational Discrete Mathematics, LNCS 2122, pp. 119–135, 2001. http://page.mi.fuberlin.de/rote/Papers/pdf/Divisionfree+algorithms.pdf 
[DW95]  Andreas W.M. Dress, Walter Wenzel, A Simple Proof of an Identity Concerning Pfaffians of Skew Symmetric Matrices, Advances in Mathematics, volume 112, Issue 1, April 1995, pp. 120134. http://www.sciencedirect.com/science/article/pii/S0001870885710298 
Todo
Implement faster algorithms, including a divisionfree one. Does [Rote2001], section 3.3 give one?
Check the implementation of the matchings used here for performance?
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
AUTHORS:
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,)
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 to matrix_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()
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')
Another random plot, but over GF(389):
sage: A = random_matrix(GF(389), 10)
sage: A.plot(cmap='Oranges')
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
AUTHORS:
Randomize density proportion of the entries of this matrix, leaving the rest unchanged.
Note
We actually choose at random density proportion of entries of the matrix and set them to random elements. It’s possible that the same position can be chosen multiple times, especially for a very small matrix.
INPUT:
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
[ 1/2*x^2  x  12 1/2*x^2  1/95*x  1/2 0]
[5/2*x^2 + 2/3*x  1/4 0 0]
[ x^2 + 2/3*x 0 0]
Now we randomize all the entries of the resulting matrix:
sage: a.randomize()
sage: a
[ 1/3*x^2  x + 1 x^2 + 1 x^2  x]
[ 1/14*x^2  x  1/4 4*x  1/5 1/4*x^2  1/2*x + 4]
[ 1/9*x^2 + 5/2*x  3 x^2 + 3/2*x + 1 2/7*x^2  x  1/2]
We create the zero matrix over the integers:
sage: a = matrix(ZZ, 2); a
[0 0]
[0 0]
Then we randomize it; the x and y parameters, which determine the size of the random elements, are passed onto the ZZ random_element method.
sage: a.randomize(x=2^64, y=2^64)
sage: a
[12401200298100116246 1709403521783430739]
[ 4417091203680573707 17094769731745295000]
Returns the rational canonical form, also known as Frobenius form.
INPUT:
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 changeofbasis 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 [STORJOHANNEMAIL].
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 the polys 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 and D below were built as similar matrices, while E was built to be just slightly different. All three matrices have equal characteristic polynomials though E‘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*a + 6, a + 6, 2, 0, 3*a + 6, 5*a + 4, 5*a + 6, 5*a + 2, 3, 4*a + 2],
... [ 3*a, 6*a, 3*a, 4*a, 4*a + 4, 3*a + 6, 6*a, 4, 3*a + 4, 6*a + 2, 4*a],
... [4*a + 5, a + 1, 4*a + 3, 6*a + 5, 5*a + 2, 5*a + 2, 6*a, 4*a + 6, 6*a + 4, 5*a + 3, 3*a + 1],
... [ 3*a, 6*a, 4*a + 1, 6*a + 2, 2*a + 5, 4*a + 6, 2, a + 5, 2*a + 4, 2*a + 1, 2*a + 1],
... [4*a + 5, 3*a + 3, 6, 4*a + 1, 4*a + 3, 6*a + 3, 6, 3*a + 3, 3, a + 3, 0],
... [6*a + 6, a + 4, 2*a + 6, 3*a + 5, 4*a + 3, 2, a, 3*a + 4, 5*a, 2*a + 5, 4*a + 3],
... [3*a + 5, 6*a + 2, 4*a, a + 5, 0, 5*a, 6*a + 5, 2*a + 1, 3*a + 1, 3*a + 5, 4*a + 2],
... [3*a + 2, a + 3, 3*a + 6, a, 3*a + 5, 5*a + 1, 3*a + 2, a + 3, a + 2, 6*a + 1, 3*a + 3],
... [6*a + 6, 5*a + 1, 4*a, 2, 5*a + 5, 3*a + 5, 3*a + 1, 2*a, 2*a, 2*a + 4, 4*a + 2]])
sage: A.rational_form()
[ a + 2 0 0 0 0 0 0 0 0 0 0]
[++]
[ 0 0 0 a + 6 0 0 0 0 0 0 0]
[ 0 1 0 6*a + 4 0 0 0 0 0 0 0]
[ 0 0 1 6*a + 4 0 0 0 0 0 0 0]
[++]
[ 0 0 0 0 0 0 0 0 0 0 2*a]
[ 0 0 0 0 1 0 0 0 0 0 6*a + 3]
[ 0 0 0 0 0 1 0 0 0 0 6*a + 1]
[ 0 0 0 0 0 0 1 0 0 0 a + 2]
[ 0 0 0 0 0 0 0 1 0 0 a + 6]
[ 0 0 0 0 0 0 0 0 1 0 2*a + 1]
[ 0 0 0 0 0 0 0 0 0 1 2*a + 1]
sage: invariants = A.rational_form(format='invariants')
sage: invariants
[[6*a + 5, 1], [6*a + 1, a + 3, a + 3, 1], [5*a, a + 4, a + 6, 6*a + 5, 6*a + 1, 5*a + 6, 5*a + 6, 1]]
sage: R.<x> = F[]
sage: polys = [R(p) for p in invariants]
sage: [p.factor() for p in polys]
[x + 6*a + 5, (x + 6*a + 5) * (x^2 + (2*a + 5)*x + 5*a), (x + 6*a + 5) * (x^2 + (2*a + 5)*x + 5*a)^3]
sage: polys[1] == A.minimal_polynomial()
True
sage: prod(polys) == A.characteristic_polynomial()
True
sage: A.eigenvalues()
Traceback (most recent call last):
...
NotImplementedError: eigenvalues() is not implemented for matrices with eigenvalues that are not in the fraction field of the base ring or in QQbar
Companion matrices may be selected as any one of four different types. See the documentation for the companion matrix constructor, sage.matrix.constructor.companion_matrix(), for more information.
sage: A = matrix(QQ, [[35, 18, 2, 45],
... [22, 22, 12, 16],
... [ 5, 12, 12, 4],
... [16, 6, 4, 23]])
sage: A.rational_form(format='right')
[ 2 0 0 0]
[+]
[ 0 0 0 10]
[ 0 1 0 1]
[ 0 0 1 0]
sage: A.rational_form(format='bottom')
[ 2 0 0 0]
[+]
[ 0 0 1 0]
[ 0 0 0 1]
[ 010 1 0]
sage: A.rational_form(format='left')
[ 2 0 0 0]
[+]
[ 0 0 1 0]
[ 01 0 1]
[ 010 0 0]
sage: A.rational_form(format='top')
[ 2 0 0 0]
[+]
[ 0 0 1 10]
[ 0 1 0 0]
[ 0 0 1 0]
TESTS:
sage: A = matrix(QQ, 2, 3, range(6))
sage: A.rational_form()
Traceback (most recent call last):
...
TypeError: matrix must be square, not 2 x 3
sage: A = matrix(Integers(6), 2, 2, range(4))
sage: A.rational_form()
Traceback (most recent call last):
...
TypeError: matrix entries must come from an exact field, not Ring of integers modulo 6
sage: A = matrix(RDF, 2, 2, range(4))
sage: A.rational_form()
Traceback (most recent call last):
...
TypeError: matrix entries must come from an exact field, not Real Double Field
sage: A = matrix(QQ, 2, range(4))
sage: A.rational_form(format='junk')
Traceback (most recent call last):
...
ValueError: 'format' keyword must be 'right', 'bottom', 'left', 'top' or 'invariants', not junk
sage: A = matrix(QQ, 2, range(4))
sage: A.rational_form(subdivide='garbage')
Traceback (most recent call last):
...
ValueError: 'subdivide' keyword must be True or False, not garbage
Citations
[STORJOHANNEMAIL] 

AUTHOR:
Returns the matrix that defines the action of self on the chosen basis for the invariant subspace V. If V is an ambient, returns self (not a copy of self).
INPUT:
OUTPUT: a matrix
Warning
This function returns an nxn matrix, where V has dimension n. It does not check that V is in fact invariant under self, unless check is True.
EXAMPLES:
sage: V = VectorSpace(QQ, 3)
sage: M = MatrixSpace(QQ, 3)
sage: A = M([1,2,0, 3,4,0, 0,0,0])
sage: W = V.subspace([[1,0,0], [0,1,0]])
sage: A.restrict(W)
[1 2]
[3 4]
sage: A.restrict(W, check=True)
[1 2]
[3 4]
We illustrate the warning about invariance not being checked by default, by giving a noninvariant subspace. With the default check=False this function returns the ‘restriction’ matrix, which is meaningless as check=True reveals.
sage: W2 = V.subspace([[1,0,0], [0,1,1]])
sage: A.restrict(W2, check=False)
[1 2]
[3 4]
sage: A.restrict(W2, check=True)
Traceback (most recent call last):
...
ArithmeticError: subspace is not invariant under matrix
Suppose that self defines a linear map from some domain to a codomain that contains \(V\) and that the image of self is contained in \(V\). This function returns a new matrix \(A\) that represents this linear map but as a map to \(V\), in the sense that if \(x\) is in the domain, then \(xA\) is the linear combination of the elements of the basis of \(V\) that equals v*self.
INPUT:
See also
EXAMPLES:
sage: A = matrix(QQ,3,[1..9])
sage: V = (QQ^3).span([[1,2,3], [7,8,9]]); V
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[ 1 0 1]
[ 0 1 2]
sage: z = vector(QQ,[1,2,5])
sage: B = A.restrict_codomain(V); B
[1 2]
[4 5]
[7 8]
sage: z*B
(44, 52)
sage: z*A
(44, 52, 60)
sage: 44*V.0 + 52*V.1
(44, 52, 60)
Compute the matrix relative to the basis for V on the domain obtained by restricting self to V, but not changing the codomain of the matrix. This is the matrix whose rows are the images of the basis for V.
INPUT:
See also
EXAMPLES:
sage: V = QQ^3
sage: A = matrix(QQ,3,[1,2,0, 3,4,0, 0,0,0])
sage: W = V.subspace([[1,0,0], [1,2,3]])
sage: A.restrict_domain(W)
[1 2 0]
[3 4 0]
sage: W2 = V.subspace_with_basis([[1,0,0], [1,2,3]])
sage: A.restrict_domain(W2)
[ 1 2 0]
[ 7 10 0]
Return matrices D and P, where D is a diagonal matrix of eigenvalues and P is the corresponding matrix where the columns are corresponding eigenvectors (or zero vectors) so that self*P = P*D.
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: A = jordan_block(2,3)
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
TESTS:
For matrices with floating point entries, some platforms will return eigenvectors that are negatives of those returned by the majority of platforms. This test accounts for that possibility. Running this test independently, without adjusting the eigenvectors could indicate this situation on your hardware.
sage: B = matrix(QQ, 3, 3, range(9))
sage: em = B.change_ring(RDF).eigenmatrix_right()
sage: evalues = em[0]; evalues.dense_matrix().zero_at(2e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 0.0]
[ 0.0 0.0 0.0]
sage: evectors = em[1];
sage: for i in range(3):
... scale = evectors[0,i].sign()
... evectors.rescale_col(i, scale)
sage: evectors
[ 0.164763817... 0.799699663... 0.408248290...]
[ 0.505774475... 0.104205787... 0.816496580...]
[ 0.846785134... 0.591288087... 0.408248290...]
Compute the right eigenspaces of a matrix.
Note that eigenspaces_right() and right_eigenspaces() are identical methods. Here “right” refers to the eigenvectors being placed to the right of the matrix.
INPUT:
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)) < 1e10
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 doubleprecision floatingpoint 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().zero_at(1e15)
[ 13.3484692283 0.0 0.0]
[ 0.0 1.34846922835 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]
TESTS:
sage: B = matrix(QQ, 2, 3, range(6))
sage: B.eigenspaces_right()
Traceback (most recent call last):
...
TypeError: matrix must be square, not 2 x 3
sage: B = matrix(QQ, 4, 4, range(16))
sage: B.eigenspaces_right(format='junk')
Traceback (most recent call last):
...
ValueError: format keyword must be None, 'all' or 'galois', not junk
sage: B.eigenspaces_right(algebraic_multiplicity='garbage')
Traceback (most recent call last):
...
ValueError: algebraic_multiplicity keyword must be True or False
Compute the right eigenvectors of a matrix.
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)) < 1e10
True
Returns the right kernel of this matrix, as a vector space or free module. This is the set of vectors x such that self*x = 0.
Note
For the left kernel, use left_kernel(). The method kernel() is exactly equal to left_kernel().
INPUT:
OUTPUT:
A vector space or free module whose degree equals the number of columns in self and contains all the vectors x such that self*x = 0.
If self has 0 columns, the kernel has dimension 0, while if self has 0 rows the kernel is the entire ambient vector space.
The result is cached. Requesting the right kernel a second time, but with a different basis format, will return the cached result with the format from the first computation.
Note
For more detailed documentation on the selection of algorithms used and a more flexible method for computing a basis matrix for a right kernel (rather than computing a vector space), see right_kernel_matrix(), which powers the computations for this method.
EXAMPLES:
sage: A = matrix(QQ, [[0, 0, 1, 2, 2, 5, 3],
... [1, 5, 2, 2, 1, 7, 5],
... [0, 0, 2, 3, 3, 8, 5],
... [1, 5, 0, 1, 2, 1, 0]])
sage: K = A.right_kernel(); K
Vector space of degree 7 and dimension 4 over Rational Field
Basis matrix:
[ 1 0 0 0 1 1 1]
[ 0 1 0 0 5 5 5]
[ 0 0 1 0 1 2 3]
[ 0 0 0 1 0 1 1]
sage: A*K.basis_matrix().transpose() == zero_matrix(QQ, 4, 4)
True
The default is basis vectors that form a matrix in echelon form. A “pivot basis” instead has a basis matrix where the columns of an identity matrix are in the locations of the nonpivot columns of the original matrix. This alternate format is available whenever the base ring is a field.
sage: A = matrix(QQ, [[0, 0, 1, 2, 2, 5, 3],
... [1, 5, 2, 2, 1, 7, 5],
... [0, 0, 2, 3, 3, 8, 5],
... [1, 5, 0, 1, 2, 1, 0]])
sage: A.rref()
[ 1 5 0 0 1 1 1]
[ 0 0 1 0 0 1 1]
[ 0 0 0 1 1 2 1]
[ 0 0 0 0 0 0 0]
sage: A.nonpivots()
(1, 4, 5, 6)
sage: K = A.right_kernel(basis='pivot'); K
Vector space of degree 7 and dimension 4 over Rational Field
User basis matrix:
[ 5 1 0 0 0 0 0]
[1 0 0 1 1 0 0]
[1 0 1 2 0 1 0]
[ 1 0 1 1 0 0 1]
sage: A*K.basis_matrix().transpose() == zero_matrix(QQ, 4, 4)
True
Matrices may have any field as a base ring. Number fields are computed by PARI library code, matrices over \(GF(2)\) are computed by the M4RI library, and matrices over the rationals are computed by the IML library. For any of these specialized cases, generalpurpose code can be called instead with the keyword setting algorithm='generic'.
Over an arbitrary field, with two basis formats. Same vector space, different bases.
sage: F.<a> = FiniteField(5^2)
sage: A = matrix(F, 3, 4, [[ 1, a, 1+a, a^3+a^5],
... [ a, a^4, a+a^4, a^4+a^8],
... [a^2, a^6, a^2+a^6, a^5+a^10]])
sage: K = A.right_kernel(); K
Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^2
Basis matrix:
[ 1 0 3*a + 4 2*a + 2]
[ 0 1 2*a 3*a + 3]
sage: A*K.basis_matrix().transpose() == zero_matrix(F, 3, 2)
True
sage: B = copy(A)
sage: P = B.right_kernel(basis = 'pivot'); P
Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^2
User basis matrix:
[ 4 4 1 0]
[ a + 2 3*a + 3 0 1]
sage: B*P.basis_matrix().transpose() == zero_matrix(F, 3, 2)
True
sage: K == P
True
Over number fields, PARI is used by default, but generalpurpose code can be requested. Same vector space, same bases, different code.:
sage: Q = QuadraticField(7)
sage: a = Q.gen(0)
sage: A = matrix(Q, [[ 2, 5a, 15a, 16+4*a],
... [2+a, a, 7 + 5*a, 3+3*a]])
sage: K = A.right_kernel(algorithm='default'); K
Vector space of degree 4 and dimension 2 over Number Field in a with defining polynomial x^2 + 7
Basis matrix:
[ 1 0 7/88*a + 3/88 3/176*a  39/176]
[ 0 1 1/88*a  13/88 13/176*a  7/176]
sage: A*K.basis_matrix().transpose() == zero_matrix(Q, 2, 2)
True
sage: B = copy(A)
sage: G = A.right_kernel(algorithm='generic'); G
Vector space of degree 4 and dimension 2 over Number Field in a with defining polynomial x^2 + 7
Basis matrix:
[ 1 0 7/88*a + 3/88 3/176*a  39/176]
[ 0 1 1/88*a  13/88 13/176*a  7/176]
sage: B*G.basis_matrix().transpose() == zero_matrix(Q, 2, 2)
True
sage: K == G
True
For matrices over the integers, several options are possible. The basis can be an LLLreduced basis or an echelon basis. The pivot basis isnot available. A heuristic will decide whether to use a padic algorithm from the IML library or an algorithm from the PARI library. Note how specifying the algorithm can mildly influence the LLL basis.
sage: A = matrix(ZZ, [[0, 1, 1, 2, 9, 4, 4],
... [1, 1, 0, 2, 7, 1, 6],
... [2, 0, 1, 0, 1, 5, 2],
... [1, 1, 1, 3, 10, 10, 9],
... [1, 2, 0, 3, 7, 1, 6]])
sage: A.right_kernel(basis='echelon')
Free module of degree 7 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1 5 8 3 1 1 1]
[ 0 11 19 5 2 3 3]
sage: B = copy(A)
sage: B.right_kernel(basis='LLL')
Free module of degree 7 and rank 2 over Integer Ring
User basis matrix:
[ 2 1 3 1 0 1 1]
[5 3 2 5 1 1 1]
sage: C = copy(A)
sage: C.right_kernel(basis='pivot')
Traceback (most recent call last):
...
ValueError: pivot basis only available over a field, not over Integer Ring
sage: D = copy(A)
sage: D.right_kernel(algorithm='pari')
Free module of degree 7 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1 5 8 3 1 1 1]
[ 0 11 19 5 2 3 3]
sage: E = copy(A)
sage: E.right_kernel(algorithm='padic', basis='LLL')
Free module of degree 7 and rank 2 over Integer Ring
User basis matrix:
[2 1 3 1 0 1 1]
[ 5 3 2 5 1 1 1]
Besides the integers, rings may be as general as principal ideal domains. Results are then free modules.
sage: R.<y> = QQ[]
sage: A = matrix(R, [[ 1, y, 1+y^2],
... [y^3, y^2, 2*y^3]])
sage: K = A.right_kernel(algorithm='default', basis='echelon'); K
Free module of degree 3 and rank 1 over Univariate Polynomial Ring in y over Rational Field
Echelon basis matrix:
[1 y 1]
sage: A*K.basis_matrix().transpose() == zero_matrix(ZZ, 2, 1)
True
It is possible to compute a kernel for a matrix over an integral domain which is not a PID, but usually this will fail.
sage: D.<x> = ZZ[]
sage: A = matrix(D, 2, 2, [[x^2  x, x + 5],
... [x^2  8, x + 2]])
sage: A.right_kernel()
Traceback (most recent call last):
...
ArithmeticError: Ideal Ideal (x^2  x, x^2  8) of Univariate Polynomial Ring in x over Integer Ring not principal
Matrices over noncommutative rings are not a good idea either. These are the “usual” quaternions.
sage: Q.<i,j,k> = QuaternionAlgebra(1,1)
sage: A = matrix(Q, 2, [i,j,1,k])
sage: A.right_kernel()
Traceback (most recent call last):
...
NotImplementedError: Cannot compute a matrix kernel over Quaternion Algebra (1, 1) with base ring Rational Field
Sparse matrices, over the rationals and the integers, use the same routines as the dense versions.
sage: A = matrix(ZZ, [[0, 1, 1, 1, 2],
... [1, 2, 0, 1, 3],
... [1, 2, 0, 1, 3]],
... sparse=True)
sage: A.right_kernel()
Free module of degree 5 and rank 3 over Integer Ring
Echelon basis matrix:
[ 1 0 0 2 1]
[ 0 1 0 1 1]
[ 0 0 1 3 1]
sage: B = A.change_ring(QQ)
sage: B.is_sparse()
True
sage: B.right_kernel()
Vector space of degree 5 and dimension 3 over Rational Field
Basis matrix:
[ 1 0 0 2 1]
[ 0 1 0 1 1]
[ 0 0 1 3 1]
With no columns, the kernel can only have dimension zero. With no rows, every possible vector is in the kernel.
sage: A = matrix(QQ, 2, 0)
sage: A.right_kernel()
Vector space of degree 0 and dimension 0 over Rational Field
Basis matrix:
[]
sage: A = matrix(QQ, 0, 2)
sage: A.right_kernel()
Vector space of degree 2 and dimension 2 over Rational Field
Basis matrix:
[1 0]
[0 1]
Every vector is in the kernel of a zero matrix, the dimension is the number of columns.
sage: A = zero_matrix(QQ, 10, 20)
sage: A.right_kernel()
Vector space of degree 20 and dimension 20 over Rational Field
Basis matrix:
20 x 20 dense matrix over Rational Field
Results are cached as the right kernel of the matrix. Subsequent requests for the right kernel will return the cached result, without regard for new values of the algorithm or format keyword. Work with a copy if you need a new right kernel, or perhaps investigate the right_kernel_matrix() method, which does not cache its results and is more flexible.
sage: A = matrix(QQ, 3, range(9))
sage: K1 = A.right_kernel(basis='echelon')
sage: K1
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 2 1]
sage: K2 = A.right_kernel(basis='pivot')
sage: K2
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 2 1]
sage: K1 is K2
True
sage: B = copy(A)
sage: K3 = B.kernel(basis='pivot')
sage: K3
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[ 1 2 1]
sage: K3 is K1
False
sage: K3 == K1
True
Returns a matrix whose rows form a basis for the right kernel of self.
INPUT:
OUTPUT:
A matrix X whose rows are an independent set spanning the right kernel of self. So self*X.transpose() is a zero matrix.
The output varies depending on the choice of algorithm and the format chosen by basis.
The results of this routine are not cached, so you can call it again with different options to get possibly different output (like the basis format). Conversely, repeated calls on the same matrix will always start from scratch.
Note
If you want to get the most basic description of a kernel, with a minimum of overhead, then ask for the right kernel matrix with the basis format requested as ‘computed’. You are then free to work with the output for whatever purpose. For a left kernel, call this method on the transpose of your matrix.
For greater convenience, plus cached results, request an actual vector space or free module with right_kernel() or left_kernel().
EXAMPLES:
Over the Rational Numbers:
Kernels are computed by the IML library in _right_kernel_matrix(). Setting the \(algorithm\) keyword to ‘default’, ‘padic’ or unspecified will yield the same result, as there is no optional behavior. The ‘computed’ format of the basis vectors are exactly the negatives of the vectors in the ‘pivot’ format.
sage: A = matrix(QQ, [[1, 0, 1, 3, 1],
... [5, 1, 0, 7, 3],
... [0, 1, 4, 6, 2],
... [4, 1, 0, 6, 2]])
sage: C = A.right_kernel_matrix(algorithm='default', basis='computed'); C
[1 2 2 1 0]
[ 1 2 0 0 1]
sage: A*C.transpose() == zero_matrix(QQ, 4, 2)
True
sage: P = A.right_kernel_matrix(algorithm='padic', basis='pivot'); P
[ 1 2 2 1 0]
[1 2 0 0 1]
sage: A*P.transpose() == zero_matrix(QQ, 4, 2)
True
sage: C == P
True
sage: E = A.right_kernel_matrix(algorithm='default', basis='echelon'); E
[ 1 0 1 1/2 1/2]
[ 0 1 1/2 1/4 1/4]
sage: A*E.transpose() == zero_matrix(QQ, 4, 2)
True
Since the rationals are a field, we can call the general code available for any field by using the ‘generic’ keyword.
sage: A = matrix(QQ, [[1, 0, 1, 3, 1],
... [5, 1, 0, 7, 3],
... [0, 1, 4, 6, 2],
... [4, 1, 0, 6, 2]])
sage: G = A.right_kernel_matrix(algorithm='generic', basis='echelon'); G
[ 1 0 1 1/2 1/2]
[ 0 1 1/2 1/4 1/4]
sage: A*G.transpose() == zero_matrix(QQ, 4, 2)
True
We verify that the rational matrix code is called for both dense and sparse rational matrices, with equal result.
sage: A = matrix(QQ, [[1, 0, 1, 3, 1],
... [5, 1, 0, 7, 3],
... [0, 1, 4, 6, 2],
... [4, 1, 0, 6, 2]],
... sparse=False)
sage: B = copy(A).sparse_matrix()
sage: set_verbose(1)
sage: D = A.right_kernel(); D
verbose ...
verbose 1 (<module>) computing right kernel matrix over the rationals for 4x5 matrix
...
verbose 1 (<module>) done computing right kernel matrix over the rationals for 4x5 matrix
...
Vector space of degree 5 and dimension 2 over Rational Field
Basis matrix:
[ 1 0 1 1/2 1/2]
[ 0 1 1/2 1/4 1/4]
sage: S = B.right_kernel(); S
verbose ...
verbose 1 (<module>) computing right kernel matrix over the rationals for 4x5 matrix
...
verbose 1 (<module>) done computing right kernel matrix over the rationals for 4x5 matrix
...
Vector space of degree 5 and dimension 2 over Rational Field
Basis matrix:
[ 1 0 1 1/2 1/2]
[ 0 1 1/2 1/4 1/4]
sage: set_verbose(0)
sage: D == S
True
Over Number Fields:
Kernels are by default computed by PARI, (except for exceptions like the rationals themselves). The raw results from PARI are a pivot basis, so the \(basis\) keywords ‘computed’ and ‘pivot’ will return the same results.
sage: Q = QuadraticField(7)
sage: a = Q.gen(0)
sage: A = matrix(Q, [[2, 5a, 15a, 16+4*a],
... [2+a, a, 7 + 5*a, 3+3*a]])
sage: C = A.right_kernel_matrix(algorithm='default', basis='computed'); C
[ a 3 1 0]
[ 2 a  1 0 1]
sage: A*C.transpose() == zero_matrix(Q, 2, 2)
True
sage: P = A.right_kernel_matrix(algorithm='pari', basis='pivot'); P
[ a 3 1 0]
[ 2 a  1 0 1]
sage: A*P.transpose() == zero_matrix(Q, 2, 2)
True
sage: E = A.right_kernel_matrix(algorithm='default', basis='echelon'); E
[ 1 0 7/88*a + 3/88 3/176*a  39/176]
[ 0 1 1/88*a  13/88 13/176*a  7/176]
sage: A*E.transpose() == zero_matrix(Q, 2, 2)
True
We can bypass using PARI for number fields and use Sage’s general code for matrices over any field. The basis vectors as computed are in pivot format.
sage: Q = QuadraticField(7)
sage: a = Q.gen(0)
sage: A = matrix(Q, [[2, 5a, 15a, 16+4*a],[2+a, a, 7 + 5*a, 3+3*a]])
sage: G = A.right_kernel_matrix(algorithm='generic', basis='computed'); G
[ a 3 1 0]
[ 2 a  1 0 1]
sage: A*G.transpose() == zero_matrix(Q, 2, 2)
True
We check that number fields are handled by the right routine as part of typical right kernel computation.
sage: Q = QuadraticField(7)
sage: a = Q.gen(0)
sage: A = matrix(Q, [[2, 5a, 15a, 16+4*a],[2+a, a, 7 + 5*a, 3+3*a]])
sage: set_verbose(1)
sage: A.right_kernel(algorithm='default')
verbose ...
verbose 1 (<module>) computing right kernel matrix over a number field for 2x4 matrix
verbose 1 (<module>) done computing right kernel matrix over a number field for 2x4 matrix
...
Vector space of degree 4 and dimension 2 over Number Field in a with defining polynomial x^2 + 7
Basis matrix:
[ 1 0 7/88*a + 3/88 3/176*a  39/176]
[ 0 1 1/88*a  13/88 13/176*a  7/176]
sage: set_verbose(0)
Over the Finite Field of Order 2:
Kernels are computed by the M4RI library using PLUQ matrix decomposition in the _right_kernel_matrix() method. There are no options for the algorithm used.
sage: A = matrix(GF(2),[[0, 1, 1, 0, 0, 0],
... [1, 0, 0, 0, 1, 1,],
... [1, 0, 0, 0, 1, 1]])
sage: E = A.right_kernel_matrix(algorithm='default', format='echelon'); E
[1 0 0 0 0 1]
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 1]
sage: A*E.transpose() == zero_matrix(GF(2), 3, 4)
True
Since GF(2) is a field we can route this computation to the generic code and obtain the ‘pivot’ form of the basis. The algorithm keywords, ‘pluq’, ‘default’ and unspecified, all have the same effect as there is no optional behavior.
sage: A = matrix(GF(2),[[0, 1, 1, 0, 0, 0],
... [1, 0, 0, 0, 1, 1,],
... [1, 0, 0, 0, 1, 1]])
sage: P = A.right_kernel_matrix(algorithm='generic', basis='pivot'); P
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[1 0 0 0 1 0]
[1 0 0 0 0 1]
sage: A*P.transpose() == zero_matrix(GF(2), 3, 4)
True
sage: DP = A.right_kernel_matrix(algorithm='default', basis='pivot'); DP
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[1 0 0 0 1 0]
[1 0 0 0 0 1]
sage: A*DP.transpose() == zero_matrix(GF(2), 3, 4)
True
sage: A.right_kernel_matrix(algorithm='pluq', basis='echelon')
[1 0 0 0 0 1]
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 1]
We test that the mod 2 code is called for matrices over GF(2).
sage: A = matrix(GF(2),[[0, 1, 1, 0, 0, 0],
... [1, 0, 0, 0, 1, 1,],
... [1, 0, 0, 0, 1, 1]])
sage: set_verbose(1)
sage: A.right_kernel(algorithm='default')
verbose ...
verbose 1 (<module>) computing right kernel matrix over integers mod 2 for 3x6 matrix
verbose 1 (<module>) done computing right kernel matrix over integers mod 2 for 3x6 matrix
...
Vector space of degree 6 and dimension 4 over Finite Field of size 2
Basis matrix:
[1 0 0 0 0 1]
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 1]
sage: set_verbose(0)
Over Arbitrary Fields:
For kernels over fields not listed above, totally general code will compute a set of basis vectors in the pivot format. These could be returned as a basis in echelon form.
sage: F.<a> = FiniteField(5^2)
sage: A = matrix(F, 3, 4, [[ 1, a, 1+a, a^3+a^5],
... [ a, a^4, a+a^4, a^4+a^8],
... [a^2, a^6, a^2+a^6, a^5+a^10]])
sage: P = A.right_kernel_matrix(algorithm='default', basis='pivot'); P
[ 4 4 1 0]
[ a + 2 3*a + 3 0 1]
sage: A*P.transpose() == zero_matrix(F, 3, 2)
True
sage: E = A.right_kernel_matrix(algorithm='default', basis='echelon'); E
[ 1 0 3*a + 4 2*a + 2]
[ 0 1 2*a 3*a + 3]
sage: A*E.transpose() == zero_matrix(F, 3, 2)
True
This general code can be requested for matrices over any field with the algorithm keyword ‘generic’. Normally, matrices over the rationals would be handled by specific routines from the IML library. The default format is an echelon basis, but a pivot basis may be requested, which is identical to the computed basis.
sage: A = matrix(QQ, 3, 4, [[1,3,2,4],
... [2,0,2,2],
... [1,1,2,0]])
sage: G = A.right_kernel_matrix(algorithm='generic'); G
[ 1 0 1/2 1/2]
[ 0 1 1/2 1/2]
sage: A*G.transpose() == zero_matrix(QQ, 3, 2)
True
sage: C = A.right_kernel_matrix(algorithm='generic', basis='computed'); C
[1 1 1 0]
[1 1 0 1]
sage: A*C.transpose() == zero_matrix(QQ, 3, 2)
True
We test that the generic code is called for matrices over fields, lacking any more specific routine.
sage: F.<a> = FiniteField(5^2)
sage: A = matrix(F, 3, 4, [[ 1, a, 1+a, a^3+a^5],
... [ a, a^4, a+a^4, a^4+a^8],
... [a^2, a^6, a^2+a^6, a^5+a^10]])
sage: set_verbose(1)
sage: A.right_kernel(algorithm='default')
verbose ...
verbose 1 (<module>) computing right kernel matrix over an arbitrary field for 3x4 matrix
...
verbose 1 (<module>) done computing right kernel matrix over an arbitrary field for 3x4 matrix
...
Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^2
Basis matrix:
[ 1 0 3*a + 4 2*a + 2]
[ 0 1 2*a 3*a + 3]
sage: set_verbose(0)
Over the Integers:
Either the IML or PARI libraries are used to provide a set of basis vectors. The algorithm keyword can be used to select either, or when set to ‘default’ a heuristic will choose between the two. Results can be returned in the ‘compute’ format, straight out of the libraries. Unique to the integers, the basis vectors can be returned as an LLL basis. Note the similarities and differences in the results. The ‘pivot’ format is not available, since the integers are not a field.
sage: A = matrix(ZZ, [[8, 0, 7, 1, 3, 4, 6],
... [4, 0, 3, 4, 2, 7, 7],
... [1, 4, 6, 1, 2, 8, 5],
... [0, 3, 1, 2, 3, 6, 2]])
sage: X = A.right_kernel_matrix(algorithm='default', basis='echelon'); X
[ 1 12 3 14 3 10 1]
[ 0 35 0 25 1 31 17]
[ 0 0 7 12 3 1 8]
sage: A*X.transpose() == zero_matrix(ZZ, 4, 3)
True
sage: X = A.right_kernel_matrix(algorithm='padic', basis='LLL'); X
[ 3 1 5 7 2 3 2]
[ 3 1 2 5 5 2 6]
[ 4 13 2 7 5 7 3]
sage: A*X.transpose() == zero_matrix(ZZ, 4, 3)
True
sage: X = A.right_kernel_matrix(algorithm='pari', basis='computed'); X
[ 3 1 5 7 2 3 2]
[ 3 1 2 5 5 2 6]
[ 4 13 2 7 5 7 3]
sage: A*X.transpose() == zero_matrix(ZZ, 4, 3)
True
sage: X = A.right_kernel_matrix(algorithm='padic', basis='computed'); X
[ 265 345 178 17 297 0 0]
[242 314 163 14 271 1 0]
[ 36 47 25 1