Linear algebra¶
Choosing the base ring: exact vs. numerical¶
Every vector or matrix in Sage is defined over a base ring, and the same computation can behave very differently depending on which one you choose. The most common choices are:
Base ring |
When to use it |
|---|---|
|
Integers: exact. Determinants, Smith/Hermite normal forms, lattice problems. |
|
Rationals: exact. The default for structural questions: rank, kernel, echelon form, inverse, exact solving. |
|
Finite fields: exact modular arithmetic, e.g. for coding theory and cryptography. |
|
Floating-point with adjustable precision (53 bits by
default, any precision via |
|
Hardware double-precision floats: fast numerical linear algebra via LAPACK (SVD, QR, Cholesky, Schur, …). |
|
Ball arithmetic: floating-point intervals with certified error bounds. |
Exact rings give exact answers:
sage: A = matrix(QQ, [[1, 1/2], [1/3, 1/4]])
sage: A.det()
1/12
sage: A.inverse()
[ 3 -6]
[-4 12]
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(1), Integer(1)/Integer(2)], [Integer(1)/Integer(3), Integer(1)/Integer(4)]])
>>> A.det()
1/12
>>> A.inverse()
[ 3 -6]
[-4 12]
A = matrix(QQ, [[1, 1/2], [1/3, 1/4]]) A.det() A.inverse()
The same matrix over RDF computes with floating-point rounding:
sage: A = matrix(QQ, [[1, 1/2], [1/3, 1/4]])
sage: Ad = A.change_ring(RDF)
sage: Ad.det() # abs tol 1e-15
0.08333333333333334
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(1), Integer(1)/Integer(2)], [Integer(1)/Integer(3), Integer(1)/Integer(4)]])
>>> Ad = A.change_ring(RDF)
>>> Ad.det() # abs tol 1e-15
0.08333333333333334
A = matrix(QQ, [[1, 1/2], [1/3, 1/4]]) Ad = A.change_ring(RDF) Ad.det() # abs tol 1e-15
so exact equality tests that hold over QQ can fail over RDF:
sage: matrix(QQ, [[1/10]]) * 3 == matrix(QQ, [[3/10]])
True
sage: matrix(RDF, [[0.1]]) * 3 == matrix(RDF, [[0.3]])
False
>>> from sage.all import *
>>> matrix(QQ, [[Integer(1)/Integer(10)]]) * Integer(3) == matrix(QQ, [[Integer(3)/Integer(10)]])
True
>>> matrix(RDF, [[RealNumber('0.1')]]) * Integer(3) == matrix(RDF, [[RealNumber('0.3')]])
False
matrix(QQ, [[1/10]]) * 3 == matrix(QQ, [[3/10]]) matrix(RDF, [[0.1]]) * 3 == matrix(RDF, [[0.3]])
Over an exact ring, eigenvalues are returned as exact algebraic numbers
(the trailing ? marks an element of QQbar, printed to its known
precision); over RDF they are ordinary floats:
sage: B = matrix(QQ, [[0, 1], [2, 0]])
sage: B.eigenvalues()
[-1.414213562373095?, 1.414213562373095?]
sage: sorted(B.change_ring(RDF).eigenvalues()) # abs tol 1e-14
[-1.414213562373095, 1.4142135623730951]
>>> from sage.all import *
>>> B = matrix(QQ, [[Integer(0), Integer(1)], [Integer(2), Integer(0)]])
>>> B.eigenvalues()
[-1.414213562373095?, 1.414213562373095?]
>>> sorted(B.change_ring(RDF).eigenvalues()) # abs tol 1e-14
[-1.414213562373095, 1.4142135623730951]
B = matrix(QQ, [[0, 1], [2, 0]]) B.eigenvalues() sorted(B.change_ring(RDF).eigenvalues()) # abs tol 1e-14
Over RDF and CDF, numerical decompositions such as SVD, QR,
cholesky and schur are computed by optimized LAPACK routines, and
some of them (SVD, schur) are available over no other ring:
sage: Ad = matrix(RDF, [[1, 1/2], [1/3, 1/4]])
sage: U, S, V = Ad.SVD()
sage: S.diagonal() # abs tol 1e-14
[1.1910987498956362, 0.0699634126394936]
>>> from sage.all import *
>>> Ad = matrix(RDF, [[Integer(1), Integer(1)/Integer(2)], [Integer(1)/Integer(3), Integer(1)/Integer(4)]])
>>> U, S, V = Ad.SVD()
>>> S.diagonal() # abs tol 1e-14
[1.1910987498956362, 0.0699634126394936]
Ad = matrix(RDF, [[1, 1/2], [1/3, 1/4]]) U, S, V = Ad.SVD() S.diagonal() # abs tol 1e-14
For floating-point results carrying rigorous error bounds, use the ball
fields RBF/CBF:
sage: M = matrix(RBF, [[1, 1/3], [1/3, 1]])
sage: M.det()
[0.888888888888889 +/- ...]
>>> from sage.all import *
>>> M = matrix(RBF, [[Integer(1), Integer(1)/Integer(3)], [Integer(1)/Integer(3), Integer(1)]])
>>> M.det()
[0.888888888888889 +/- ...]
M = matrix(RBF, [[1, 1/3], [1/3, 1]]) M.det()
As a rule of thumb: start with QQ when you want exact structural
answers; switch to RDF/CDF for speed on large matrices and for
the numerical decompositions; use RR/CC when you need a
nondefault precision; and use RBF/CBF when you need
floating-point answers you can certify.
Vector spaces¶
The VectorSpace command creates a vector space class, from which
one can create a subspace. Note the basis computed by Sage is
“row reduced”.
sage: V = VectorSpace(GF(2),8)
sage: S = V.subspace([V([1,1,0,0,0,0,0,0]),V([1,0,0,0,0,1,1,0])])
sage: S.basis()
[(1, 0, 0, 0, 0, 1, 1, 0), (0, 1, 0, 0, 0, 1, 1, 0)]
sage: S.dimension()
2
>>> from sage.all import *
>>> V = VectorSpace(GF(Integer(2)),Integer(8))
>>> S = V.subspace([V([Integer(1),Integer(1),Integer(0),Integer(0),Integer(0),Integer(0),Integer(0),Integer(0)]),V([Integer(1),Integer(0),Integer(0),Integer(0),Integer(0),Integer(1),Integer(1),Integer(0)])])
>>> S.basis()
[(1, 0, 0, 0, 0, 1, 1, 0), (0, 1, 0, 0, 0, 1, 1, 0)]
>>> S.dimension()
2
V = VectorSpace(GF(2),8) S = V.subspace([V([1,1,0,0,0,0,0,0]),V([1,0,0,0,0,1,1,0])]) S.basis() S.dimension()
Matrix powers¶
How do I compute matrix powers in Sage? The syntax is illustrated by the example below.
sage: R = IntegerModRing(51)
sage: M = MatrixSpace(R,3,3)
sage: A = M([1,2,3, 4,5,6, 7,8,9])
sage: A^1000*A^1007
[ 3 3 3]
[18 0 33]
[33 48 12]
sage: A^2007
[ 3 3 3]
[18 0 33]
[33 48 12]
>>> from sage.all import *
>>> R = IntegerModRing(Integer(51))
>>> M = MatrixSpace(R,Integer(3),Integer(3))
>>> A = M([Integer(1),Integer(2),Integer(3), Integer(4),Integer(5),Integer(6), Integer(7),Integer(8),Integer(9)])
>>> A**Integer(1000)*A**Integer(1007)
<BLANKLINE>
[ 3 3 3]
[18 0 33]
[33 48 12]
>>> A**Integer(2007)
<BLANKLINE>
[ 3 3 3]
[18 0 33]
[33 48 12]
R = IntegerModRing(51) M = MatrixSpace(R,3,3) A = M([1,2,3, 4,5,6, 7,8,9]) A^1000*A^1007 A^2007
Kernels¶
The kernel is computed by applying the kernel method to the matrix object. The following examples illustrate the syntax.
sage: M = MatrixSpace(IntegerRing(),4,2)(range(8))
sage: M.kernel()
Free module of degree 4 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1 0 -3 2]
[ 0 1 -2 1]
>>> from sage.all import *
>>> M = MatrixSpace(IntegerRing(),Integer(4),Integer(2))(range(Integer(8)))
>>> M.kernel()
Free module of degree 4 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1 0 -3 2]
[ 0 1 -2 1]
M = MatrixSpace(IntegerRing(),4,2)(range(8)) M.kernel()
A kernel of dimension one over \(\QQ\):
sage: A = MatrixSpace(RationalField(),3)(range(9))
sage: A.kernel()
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2 1]
>>> from sage.all import *
>>> A = MatrixSpace(RationalField(),Integer(3))(range(Integer(9)))
>>> A.kernel()
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2 1]
A = MatrixSpace(RationalField(),3)(range(9)) A.kernel()
A trivial kernel:
sage: A = MatrixSpace(RationalField(),2)([1,2,3,4])
sage: A.kernel()
Vector space of degree 2 and dimension 0 over Rational Field
Basis matrix:
[]
sage: M = MatrixSpace(RationalField(),0,2)(0)
sage: M
[]
sage: M.kernel()
Vector space of degree 0 and dimension 0 over Rational Field
Basis matrix:
[]
sage: M = MatrixSpace(RationalField(),2,0)(0)
sage: M.kernel()
Vector space of degree 2 and dimension 2 over Rational Field
Basis matrix:
[1 0]
[0 1]
>>> from sage.all import *
>>> A = MatrixSpace(RationalField(),Integer(2))([Integer(1),Integer(2),Integer(3),Integer(4)])
>>> A.kernel()
Vector space of degree 2 and dimension 0 over Rational Field
Basis matrix:
[]
>>> M = MatrixSpace(RationalField(),Integer(0),Integer(2))(Integer(0))
>>> M
[]
>>> M.kernel()
Vector space of degree 0 and dimension 0 over Rational Field
Basis matrix:
[]
>>> M = MatrixSpace(RationalField(),Integer(2),Integer(0))(Integer(0))
>>> M.kernel()
Vector space of degree 2 and dimension 2 over Rational Field
Basis matrix:
[1 0]
[0 1]
A = MatrixSpace(RationalField(),2)([1,2,3,4]) A.kernel() M = MatrixSpace(RationalField(),0,2)(0) M M.kernel() M = MatrixSpace(RationalField(),2,0)(0) M.kernel()
Kernel of a zero matrix:
sage: A = MatrixSpace(RationalField(),2)(0)
sage: A.kernel()
Vector space of degree 2 and dimension 2 over Rational Field
Basis matrix:
[1 0]
[0 1]
>>> from sage.all import *
>>> A = MatrixSpace(RationalField(),Integer(2))(Integer(0))
>>> A.kernel()
Vector space of degree 2 and dimension 2 over Rational Field
Basis matrix:
[1 0]
[0 1]
A = MatrixSpace(RationalField(),2)(0) A.kernel()
Kernel of a non-square matrix:
sage: A = MatrixSpace(RationalField(),3,2)(range(6))
sage: A.kernel()
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2 1]
>>> from sage.all import *
>>> A = MatrixSpace(RationalField(),Integer(3),Integer(2))(range(Integer(6)))
>>> A.kernel()
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2 1]
A = MatrixSpace(RationalField(),3,2)(range(6)) A.kernel()
The 2-dimensional kernel of a matrix over a cyclotomic field:
sage: K = CyclotomicField(12); a = K.gen()
sage: M = MatrixSpace(K,4,2)([1,-1, 0,-2, 0,-a^2-1, 0,a^2-1])
sage: M
[ 1 -1]
[ 0 -2]
[ 0 -zeta12^2 - 1]
[ 0 zeta12^2 - 1]
sage: M.kernel()
Vector space of degree 4 and dimension 2 over Cyclotomic Field of order 12
and degree 4
Basis matrix:
[ 0 1 0 -2*zeta12^2]
[ 0 0 1 -2*zeta12^2 + 1]
>>> from sage.all import *
>>> K = CyclotomicField(Integer(12)); a = K.gen()
>>> M = MatrixSpace(K,Integer(4),Integer(2))([Integer(1),-Integer(1), Integer(0),-Integer(2), Integer(0),-a**Integer(2)-Integer(1), Integer(0),a**Integer(2)-Integer(1)])
>>> M
[ 1 -1]
[ 0 -2]
[ 0 -zeta12^2 - 1]
[ 0 zeta12^2 - 1]
>>> M.kernel()
Vector space of degree 4 and dimension 2 over Cyclotomic Field of order 12
and degree 4
Basis matrix:
[ 0 1 0 -2*zeta12^2]
[ 0 0 1 -2*zeta12^2 + 1]
K = CyclotomicField(12); a = K.gen() M = MatrixSpace(K,4,2)([1,-1, 0,-2, 0,-a^2-1, 0,a^2-1]) M M.kernel()
A nontrivial kernel over a complicated base field.
sage: K = FractionField(PolynomialRing(RationalField(),2,'x'))
sage: M = MatrixSpace(K, 2)([[K.gen(1),K.gen(0)], [K.gen(1), K.gen(0)]])
sage: M
[x1 x0]
[x1 x0]
sage: M.kernel()
Vector space of degree 2 and dimension 1 over Fraction Field of Multivariate
Polynomial Ring in x0, x1 over Rational Field
Basis matrix:
[ 1 -1]
>>> from sage.all import *
>>> K = FractionField(PolynomialRing(RationalField(),Integer(2),'x'))
>>> M = MatrixSpace(K, Integer(2))([[K.gen(Integer(1)),K.gen(Integer(0))], [K.gen(Integer(1)), K.gen(Integer(0))]])
>>> M
[x1 x0]
[x1 x0]
>>> M.kernel()
Vector space of degree 2 and dimension 1 over Fraction Field of Multivariate
Polynomial Ring in x0, x1 over Rational Field
Basis matrix:
[ 1 -1]
K = FractionField(PolynomialRing(RationalField(),2,'x')) M = MatrixSpace(K, 2)([[K.gen(1),K.gen(0)], [K.gen(1), K.gen(0)]]) M M.kernel()
Other methods for integer matrices are elementary_divisors,
smith_form (for the Smith normal form), echelon_form
for the Hermite normal form, frobenius for the
Frobenius normal form (rational canonical form).
There are many methods for matrices over a field such as
\(\QQ\) or a finite field: row_span, nullity,
transpose, swap_rows, matrix_from_columns,
matrix_from_rows, among many others.
See the file matrix.py for further details.
Eigenvectors and eigenvalues¶
How do you compute eigenvalues and eigenvectors using Sage?
Sage has a full range of functions for computing eigenvalues and both
left and right eigenvectors and eigenspaces. If our matrix is \(A\),
then the eigenmatrix_right (resp. eightmatrix_left) command also
gives matrices \(D\) and \(P\) such that \(AP=PD\) (resp.
\(PA=DP\).)
sage: A = matrix(QQ, [[1,1,0],[0,2,0],[0,0,3]])
sage: A
[1 1 0]
[0 2 0]
[0 0 3]
sage: A.eigenvalues()
[3, 2, 1]
sage: A.eigenvectors_right()
[(3, [(0, 0, 1)], 1), (2, [(1, 1, 0)], 1), (1, [(1, 0, 0)], 1)]
sage: A.eigenspaces_right()
[(3,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[0 0 1]),
(2,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[1 1 0]),
(1,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[1 0 0])]
sage: D, P = A.eigenmatrix_right()
sage: D
[3 0 0]
[0 2 0]
[0 0 1]
sage: P
[0 1 1]
[0 1 0]
[1 0 0]
sage: A*P == P*D
True
>>> from sage.all import *
>>> A = matrix(QQ, [[Integer(1),Integer(1),Integer(0)],[Integer(0),Integer(2),Integer(0)],[Integer(0),Integer(0),Integer(3)]])
>>> A
[1 1 0]
[0 2 0]
[0 0 3]
>>> A.eigenvalues()
[3, 2, 1]
>>> A.eigenvectors_right()
[(3, [(0, 0, 1)], 1), (2, [(1, 1, 0)], 1), (1, [(1, 0, 0)], 1)]
>>> A.eigenspaces_right()
[(3,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[0 0 1]),
(2,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[1 1 0]),
(1,
Vector space of degree 3 and dimension 1 over Rational Field
User basis matrix:
[1 0 0])]
>>> D, P = A.eigenmatrix_right()
>>> D
[3 0 0]
[0 2 0]
[0 0 1]
>>> P
[0 1 1]
[0 1 0]
[1 0 0]
>>> A*P == P*D
True
A = matrix(QQ, [[1,1,0],[0,2,0],[0,0,3]]) A A.eigenvalues() A.eigenvectors_right() A.eigenspaces_right() D, P = A.eigenmatrix_right() D P A*P == P*D
For eigenvalues outside the fraction field of the base ring of the matrix,
you can choose to have all the eigenspaces output when the algebraic closure
of the field is implemented, such as the algebraic numbers, QQbar. Or you
may request just a single eigenspace for each irreducible factor of the
characteristic polynomial, since the others may be formed through Galois
conjugation. The eigenvalues of the matrix below are $pmsqrt{3}$ and
we exhibit each possible output.
Also, currently Sage does not implement multiprecision numerical eigenvalues
and eigenvectors, so calling the eigen functions on a matrix from
CC or RR will probably give inaccurate and nonsensical results
(a warning is also printed). Eigenvalues and eigenvectors of matrices with
floating point entries (over CDF and RDF) can be obtained with the
“eigenmatrix” commands.
sage: MS = MatrixSpace(QQ, 2, 2)
sage: A = MS([1,-4,1, -1])
sage: A.eigenspaces_left(format='all')
[(-1.732050807568878?*I,
Vector space of degree 2 and dimension 1 over Algebraic Field
User basis matrix:
[ 1 -1 - 1.732050807568878?*I]),
(1.732050807568878?*I,
Vector space of degree 2 and dimension 1 over Algebraic Field
User basis matrix:
[ 1 -1 + 1.732050807568878?*I])]
sage: A.eigenspaces_left(format='galois')
[(a0,
Vector space of degree 2 and dimension 1 over Number Field in a0 with defining polynomial x^2 + 3
User basis matrix:
[ 1 a0 - 1])]
>>> from sage.all import *
>>> MS = MatrixSpace(QQ, Integer(2), Integer(2))
>>> A = MS([Integer(1),-Integer(4),Integer(1), -Integer(1)])
>>> A.eigenspaces_left(format='all')
[(-1.732050807568878?*I,
Vector space of degree 2 and dimension 1 over Algebraic Field
User basis matrix:
[ 1 -1 - 1.732050807568878?*I]),
(1.732050807568878?*I,
Vector space of degree 2 and dimension 1 over Algebraic Field
User basis matrix:
[ 1 -1 + 1.732050807568878?*I])]
>>> A.eigenspaces_left(format='galois')
[(a0,
Vector space of degree 2 and dimension 1 over Number Field in a0 with defining polynomial x^2 + 3
User basis matrix:
[ 1 a0 - 1])]
MS = MatrixSpace(QQ, 2, 2) A = MS([1,-4,1, -1]) A.eigenspaces_left(format='all') A.eigenspaces_left(format='galois')
Another approach is to use the interface with Maxima:
sage: A = maxima("matrix ([1, -4], [1, -1])")
sage: eig = A.eigenvectors()
sage: eig.sage()
[[[-I*sqrt(3), I*sqrt(3)], [1, 1]], [[[1, 1/4*I*sqrt(3) + 1/4]], [[1, -1/4*I*sqrt(3) + 1/4]]]]
>>> from sage.all import *
>>> A = maxima("matrix ([1, -4], [1, -1])")
>>> eig = A.eigenvectors()
>>> eig.sage()
[[[-I*sqrt(3), I*sqrt(3)], [1, 1]], [[[1, 1/4*I*sqrt(3) + 1/4]], [[1, -1/4*I*sqrt(3) + 1/4]]]]
A = maxima("matrix ([1, -4], [1, -1])")
eig = A.eigenvectors()
eig.sage()
This tells us that \(\vec{v}_1 = [1,(\sqrt{3}i + 1)/4]\) is an eigenvector of \(\lambda_1 = - \sqrt{3}i\) (which occurs with multiplicity one) and \(\vec{v}_2 = [1,(-\sqrt{3}i + 1)/4]\) is an eigenvector of \(\lambda_2 = \sqrt{3}i\) (which also occurs with multiplicity one).
Here are two more examples:
sage: A = maxima("matrix ([11, 0, 0], [1, 11, 0], [1, 3, 2])")
sage: A.eigenvectors()
[[[2,11],[1,2]],[[[0,0,1]],[[0,1,1/3]]]]
sage: A = maxima("matrix ([-1, 0, 0], [1, -1, 0], [1, 3, 2])")
sage: A.eigenvectors()
[[[-1,2],[2,1]],[[[0,1,-1]],[[0,0,1]]]]
>>> from sage.all import *
>>> A = maxima("matrix ([11, 0, 0], [1, 11, 0], [1, 3, 2])")
>>> A.eigenvectors()
[[[2,11],[1,2]],[[[0,0,1]],[[0,1,1/3]]]]
>>> A = maxima("matrix ([-1, 0, 0], [1, -1, 0], [1, 3, 2])")
>>> A.eigenvectors()
[[[-1,2],[2,1]],[[[0,1,-1]],[[0,0,1]]]]
A = maxima("matrix ([11, 0, 0], [1, 11, 0], [1, 3, 2])")
A.eigenvectors()
A = maxima("matrix ([-1, 0, 0], [1, -1, 0], [1, 3, 2])")
A.eigenvectors()
Warning: Notice how the ordering of the output is reversed, though the matrices are almost the same.
Finally, you can use Sage’s GAP interface as well to compute “rational” eigenvalues and eigenvectors:
sage: A = libgap([[1,2,3],[4,5,6],[7,8,9]]); A
[ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
sage: libgap(QQ).Eigenvectors(A)
[ [ 1, -2, 1 ] ]
sage: libgap(QQ).Eigenvalues(A)
[ 0 ]
>>> from sage.all import *
>>> A = libgap([[Integer(1),Integer(2),Integer(3)],[Integer(4),Integer(5),Integer(6)],[Integer(7),Integer(8),Integer(9)]]); A
[ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
>>> libgap(QQ).Eigenvectors(A)
[ [ 1, -2, 1 ] ]
>>> libgap(QQ).Eigenvalues(A)
[ 0 ]
A = libgap([[1,2,3],[4,5,6],[7,8,9]]); A libgap(QQ).Eigenvectors(A) libgap(QQ).Eigenvalues(A)
Row reduction¶
The row reduced echelon form of a matrix is computed as in the following example.
sage: M = MatrixSpace(RationalField(),2,3)
sage: A = M([1,2,3, 4,5,6])
sage: A
[1 2 3]
[4 5 6]
sage: A.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field
sage: A[0,2] = 389
sage: A
[ 1 2 389]
[ 4 5 6]
sage: A.echelon_form()
[ 1 0 -1933/3]
[ 0 1 1550/3]
>>> from sage.all import *
>>> M = MatrixSpace(RationalField(),Integer(2),Integer(3))
>>> A = M([Integer(1),Integer(2),Integer(3), Integer(4),Integer(5),Integer(6)])
>>> A
[1 2 3]
[4 5 6]
>>> A.parent()
Full MatrixSpace of 2 by 3 dense matrices over Rational Field
>>> A[Integer(0),Integer(2)] = Integer(389)
>>> A
[ 1 2 389]
[ 4 5 6]
>>> A.echelon_form()
[ 1 0 -1933/3]
[ 0 1 1550/3]
M = MatrixSpace(RationalField(),2,3) A = M([1,2,3, 4,5,6]) A A.parent() A[0,2] = 389 A A.echelon_form()
Characteristic polynomial¶
The characteristic polynomial is a Sage method for square matrices.
First a matrix over \(\ZZ\):
sage: A = MatrixSpace(IntegerRing(),2)( [[1,2], [3,4]] )
sage: f = A.charpoly()
sage: f
x^2 - 5*x - 2
sage: f.parent()
Univariate Polynomial Ring in x over Integer Ring
>>> from sage.all import *
>>> A = MatrixSpace(IntegerRing(),Integer(2))( [[Integer(1),Integer(2)], [Integer(3),Integer(4)]] )
>>> f = A.charpoly()
>>> f
x^2 - 5*x - 2
>>> f.parent()
Univariate Polynomial Ring in x over Integer Ring
A = MatrixSpace(IntegerRing(),2)( [[1,2], [3,4]] ) f = A.charpoly() f f.parent()
We compute the characteristic polynomial of a matrix over the polynomial ring \(\ZZ[a]\):
sage: R = PolynomialRing(IntegerRing(),'a'); a = R.gen()
sage: M = MatrixSpace(R,2)([[a,1], [a,a+1]])
sage: M
[ a 1]
[ a a + 1]
sage: f = M.charpoly()
sage: f
x^2 + (-2*a - 1)*x + a^2
sage: f.parent()
Univariate Polynomial Ring in x over Univariate Polynomial Ring in a over
Integer Ring
sage: M.trace()
2*a + 1
sage: M.determinant()
a^2
>>> from sage.all import *
>>> R = PolynomialRing(IntegerRing(),'a'); a = R.gen()
>>> M = MatrixSpace(R,Integer(2))([[a,Integer(1)], [a,a+Integer(1)]])
>>> M
[ a 1]
[ a a + 1]
>>> f = M.charpoly()
>>> f
x^2 + (-2*a - 1)*x + a^2
>>> f.parent()
Univariate Polynomial Ring in x over Univariate Polynomial Ring in a over
Integer Ring
>>> M.trace()
2*a + 1
>>> M.determinant()
a^2
R = PolynomialRing(IntegerRing(),'a'); a = R.gen() M = MatrixSpace(R,2)([[a,1], [a,a+1]]) M f = M.charpoly() f f.parent() M.trace() M.determinant()
We compute the characteristic polynomial of a matrix over the multi-variate polynomial ring \(\ZZ[u,v]\):
sage: R.<u,v> = PolynomialRing(ZZ,2)
sage: A = MatrixSpace(R,2)([u,v,u^2,v^2])
sage: f = A.charpoly(); f
x^2 + (-v^2 - u)*x - u^2*v + u*v^2
>>> from sage.all import *
>>> R = PolynomialRing(ZZ,Integer(2), names=('u', 'v',)); (u, v,) = R._first_ngens(2)
>>> A = MatrixSpace(R,Integer(2))([u,v,u**Integer(2),v**Integer(2)])
>>> f = A.charpoly(); f
x^2 + (-v^2 - u)*x - u^2*v + u*v^2
R.<u,v> = PolynomialRing(ZZ,2) A = MatrixSpace(R,2)([u,v,u^2,v^2]) f = A.charpoly(); f
It’s a little difficult to distinguish the variables. To fix this, we might want to rename the indeterminate “Z”, which we can easily do as follows:
sage: f = A.charpoly('Z'); f
Z^2 + (-v^2 - u)*Z - u^2*v + u*v^2
>>> from sage.all import *
>>> f = A.charpoly('Z'); f
Z^2 + (-v^2 - u)*Z - u^2*v + u*v^2
f = A.charpoly('Z'); f
Solving systems of linear equations¶
Using maxima, you can easily solve linear equations:
sage: var('a,b,c')
(a, b, c)
sage: eqn = [a+b*c==1, b-a*c==0, a+b==5]
sage: s = solve(eqn, a,b,c); s
[[a == -1/4*I*sqrt(79) + 11/4, b == 1/4*I*sqrt(79) + 9/4, c == 1/10*I*sqrt(79) + 1/10], [a == 1/4*I*sqrt(79) + 11/4, b == -1/4*I*sqrt(79) + 9/4, c == -1/10*I*sqrt(79) + 1/10]]
>>> from sage.all import *
>>> var('a,b,c')
(a, b, c)
>>> eqn = [a+b*c==Integer(1), b-a*c==Integer(0), a+b==Integer(5)]
>>> s = solve(eqn, a,b,c); s
[[a == -1/4*I*sqrt(79) + 11/4, b == 1/4*I*sqrt(79) + 9/4, c == 1/10*I*sqrt(79) + 1/10], [a == 1/4*I*sqrt(79) + 11/4, b == -1/4*I*sqrt(79) + 9/4, c == -1/10*I*sqrt(79) + 1/10]]
var('a,b,c')
eqn = [a+b*c==1, b-a*c==0, a+b==5]
s = solve(eqn, a,b,c); s
You can even nicely typeset the solution in LaTeX:
sage.: print(latex(s))
...
To have the above appear onscreen via xdvi, type view(s).
You can also solve linear equations symbolically using the
solve command:
sage: var('x,y,z,a')
(x, y, z, a)
sage: eqns = [x + z == y, 2*a*x - y == 2*a^2, y - 2*z == 2]
sage: solve(eqns, x, y, z)
[[x == a + 1, y == 2*a, z == a - 1]]
>>> from sage.all import *
>>> var('x,y,z,a')
(x, y, z, a)
>>> eqns = [x + z == y, Integer(2)*a*x - y == Integer(2)*a**Integer(2), y - Integer(2)*z == Integer(2)]
>>> solve(eqns, x, y, z)
[[x == a + 1, y == 2*a, z == a - 1]]
var('x,y,z,a')
eqns = [x + z == y, 2*a*x - y == 2*a^2, y - 2*z == 2]
solve(eqns, x, y, z)
Here is a numerical Numpy example:
sage: from numpy import arange, eye, linalg
sage: A = eye(10) ## the 10x10 identity matrix
sage: b = arange(1,11)
sage: x = linalg.solve(A,b)
>>> from sage.all import *
>>> from numpy import arange, eye, linalg
>>> A = eye(Integer(10)) ## the 10x10 identity matrix
>>> b = arange(Integer(1),Integer(11))
>>> x = linalg.solve(A,b)
from numpy import arange, eye, linalg A = eye(10) ## the 10x10 identity matrix b = arange(1,11) x = linalg.solve(A,b)
Another way to solve a system numerically is to use Sage’s octave interface:
sage: M33 = MatrixSpace(QQ,3,3)
sage: A = M33([1,2,3,4,5,6,7,8,0])
sage: V3 = VectorSpace(QQ,3)
sage: b = V3([1,2,3])
sage: octave.solve_linear_system(A,b) # optional - octave
[-0.333333, 0.666667, 0]
>>> from sage.all import *
>>> M33 = MatrixSpace(QQ,Integer(3),Integer(3))
>>> A = M33([Integer(1),Integer(2),Integer(3),Integer(4),Integer(5),Integer(6),Integer(7),Integer(8),Integer(0)])
>>> V3 = VectorSpace(QQ,Integer(3))
>>> b = V3([Integer(1),Integer(2),Integer(3)])
>>> octave.solve_linear_system(A,b) # optional - octave
[-0.333333, 0.666667, 0]
M33 = MatrixSpace(QQ,3,3) A = M33([1,2,3,4,5,6,7,8,0]) V3 = VectorSpace(QQ,3) b = V3([1,2,3]) octave.solve_linear_system(A,b) # optional - octave