Matrices (linear algebra)¶
Creating Matrices¶
The linear algebra module is designed to be as simple as possible. First, we
import and declare our first Matrix
object:
>>> from sympy.interactive.printing import init_printing
>>> init_printing(use_unicode=False, wrap_line=False, no_global=True)
>>> from sympy.matrices import Matrix, eye, zeros, ones, diag, GramSchmidt
>>> M = Matrix([[1,0,0], [0,0,0]]); M
[1 0 0]
[ ]
[0 0 0]
>>> Matrix([M, (0, 0, 1)])
[1 0 0 ]
[ ]
[0 0 0 ]
[ ]
[0 0 1]
>>> Matrix([[1, 2, 3]])
[1 2 3]
>>> Matrix([1, 2, 3])
[1]
[ ]
[2]
[ ]
[3]
In addition to creating a matrix from a list of appropriatelysized lists and/or matrices, SymPy also supports more advanced methods of matrix creation including a single list of values and dimension inputs:
>>> Matrix(2, 3, [1, 2, 3, 4, 5, 6])
[1 2 3]
[ ]
[4 5 6]
More interesting (and useful), is the ability to use a 2variable function
(or lambda
) to create a matrix. Here we create an indicator function which
is 1 on the diagonal and then use it to make the identity matrix:
>>> def f(i,j):
... if i == j:
... return 1
... else:
... return 0
...
>>> Matrix(4, 4, f)
[1 0 0 0]
[ ]
[0 1 0 0]
[ ]
[0 0 1 0]
[ ]
[0 0 0 1]
Finally let’s use lambda
to create a 1line matrix with 1’s in the even
permutation entries:
>>> Matrix(3, 4, lambda i,j: 1  (i+j) % 2)
[1 0 1 0]
[ ]
[0 1 0 1]
[ ]
[1 0 1 0]
There are also a couple of special constructors for quick matrix construction:
eye
is the identity matrix, zeros
and ones
for matrices of all
zeros and ones, respectively, and diag
to put matrices or elements along
the diagonal:
>>> eye(4)
[1 0 0 0]
[ ]
[0 1 0 0]
[ ]
[0 0 1 0]
[ ]
[0 0 0 1]
>>> zeros(2)
[0 0]
[ ]
[0 0]
>>> zeros(2, 5)
[0 0 0 0 0]
[ ]
[0 0 0 0 0]
>>> ones(3)
[1 1 1]
[ ]
[1 1 1]
[ ]
[1 1 1]
>>> ones(1, 3)
[1 1 1]
>>> diag(1, Matrix([[1, 2], [3, 4]]))
[1 0 0]
[ ]
[0 1 2]
[ ]
[0 3 4]
Basic Manipulation¶
While learning to work with matrices, let’s choose one where the entries are readily identifiable. One useful thing to know is that while matrices are 2dimensional, the storage is not and so it is allowable  though one should be careful  to access the entries as if they were a 1d list.
>>> M = Matrix(2, 3, [1, 2, 3, 4, 5, 6])
>>> M[4]
5
Now, the more standard entry access is a pair of indices which will always return the value at the corresponding row and column of the matrix:
>>> M[1, 2]
6
>>> M[0, 0]
1
>>> M[1, 1]
5
Since this is Python we’re also able to slice submatrices; slices always give a matrix in return, even if the dimension is 1 x 1:
>>> M[0:2, 0:2]
[1 2]
[ ]
[4 5]
>>> M[2:2, 2]
[]
>>> M[:, 2]
[3]
[ ]
[6]
>>> M[:1, 2]
[3]
In the second example above notice that the slice 2:2 gives an empty range. Note also (in keeping with 0based indexing of Python) the first row/column is 0.
You cannot access rows or columns that are not present unless they are in a slice:
>>> M[:, 10] # the 10th column (not there)
Traceback (most recent call last):
...
IndexError: Index out of range: a[[0, 10]]
>>> M[:, 10:11] # the 10th column (if there)
[]
>>> M[:, :10] # all columns up to the 10th
[1 2 3]
[ ]
[4 5 6]
Slicing an empty matrix works as long as you use a slice for the coordinate that has no size:
>>> Matrix(0, 3, [])[:, 1]
[]
Slicing gives a copy of what is sliced, so modifications of one object do not affect the other:
>>> M2 = M[:, :]
>>> M2[0, 0] = 100
>>> M[0, 0] == 100
False
Notice that changing M2
didn’t change M
. Since we can slice, we can also assign
entries:
>>> M = Matrix(([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]))
>>> M
[1 2 3 4 ]
[ ]
[5 6 7 8 ]
[ ]
[9 10 11 12]
[ ]
[13 14 15 16]
>>> M[2,2] = M[0,3] = 0
>>> M
[1 2 3 0 ]
[ ]
[5 6 7 8 ]
[ ]
[9 10 0 12]
[ ]
[13 14 15 16]
as well as assign slices:
>>> M = Matrix(([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]))
>>> M[2:,2:] = Matrix(2,2,lambda i,j: 0)
>>> M
[1 2 3 4]
[ ]
[5 6 7 8]
[ ]
[9 10 0 0]
[ ]
[13 14 0 0]
All the standard arithmetic operations are supported:
>>> M = Matrix(([1,2,3],[4,5,6],[7,8,9]))
>>> M  M
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
>>> M + M
[2 4 6 ]
[ ]
[8 10 12]
[ ]
[14 16 18]
>>> M * M
[30 36 42 ]
[ ]
[66 81 96 ]
[ ]
[102 126 150]
>>> M2 = Matrix(3,1,[1,5,0])
>>> M*M2
[11]
[ ]
[29]
[ ]
[47]
>>> M**2
[30 36 42 ]
[ ]
[66 81 96 ]
[ ]
[102 126 150]
As well as some useful vector operations:
>>> M.row_del(0)
>>> M
[4 5 6]
[ ]
[7 8 9]
>>> M.col_del(1)
>>> M
[4 6]
[ ]
[7 9]
>>> v1 = Matrix([1,2,3])
>>> v2 = Matrix([4,5,6])
>>> v3 = v1.cross(v2)
>>> v1.dot(v2)
32
>>> v2.dot(v3)
0
>>> v1.dot(v3)
0
Recall that the row_del()
and col_del()
operations don’t return a value  they
simply change the matrix object. We can also ‘’glue’’ together matrices of the
appropriate size:
>>> M1 = eye(3)
>>> M2 = zeros(3, 4)
>>> M1.row_join(M2)
[1 0 0 0 0 0 0]
[ ]
[0 1 0 0 0 0 0]
[ ]
[0 0 1 0 0 0 0]
>>> M3 = zeros(4, 3)
>>> M1.col_join(M3)
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
[ ]
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
Operations on entries¶
We are not restricted to having multiplication between two matrices:
>>> M = eye(3)
>>> 2*M
[2 0 0]
[ ]
[0 2 0]
[ ]
[0 0 2]
>>> 3*M
[3 0 0]
[ ]
[0 3 0]
[ ]
[0 0 3]
but we can also apply functions to our matrix entries using applyfunc()
. Here we’ll declare a function that double any input number. Then we apply it to the 3x3 identity matrix:
>>> f = lambda x: 2*x
>>> eye(3).applyfunc(f)
[2 0 0]
[ ]
[0 2 0]
[ ]
[0 0 2]
One more useful matrixwide entry application function is the substitution function. Let’s declare a matrix with symbolic entries then substitute a value. Remember we can substitute anything  even another symbol!:
>>> from sympy import Symbol
>>> x = Symbol('x')
>>> M = eye(3) * x
>>> M
[x 0 0]
[ ]
[0 x 0]
[ ]
[0 0 x]
>>> M.subs(x, 4)
[4 0 0]
[ ]
[0 4 0]
[ ]
[0 0 4]
>>> y = Symbol('y')
>>> M.subs(x, y)
[y 0 0]
[ ]
[0 y 0]
[ ]
[0 0 y]
Linear algebra¶
Now that we have the basics out of the way, let’s see what we can do with the actual matrices. Of course, one of the first things that comes to mind is the determinant:
>>> M = Matrix(( [1, 2, 3], [3, 6, 2], [2, 0, 1] ))
>>> M.det()
28
>>> M2 = eye(3)
>>> M2.det()
1
>>> M3 = Matrix(( [1, 0, 0], [1, 0, 0], [1, 0, 0] ))
>>> M3.det()
0
Another common operation is the inverse: In SymPy, this is computed by Gaussian elimination by default (for dense matrices) but we can specify it be done by \(LU\) decomposition as well:
>>> M2.inv()
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
>>> M2.inv(method="LU")
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
>>> M.inv(method="LU")
[3/14 1/14 1/2 ]
[ ]
[1/28 5/28 1/4]
[ ]
[ 3/7 1/7 0 ]
>>> M * M.inv(method="LU")
[1 0 0]
[ ]
[0 1 0]
[ ]
[0 0 1]
We can perform a \(QR\) factorization which is handy for solving systems:
>>> A = Matrix([[1,1,1],[1,1,3],[2,3,4]])
>>> Q, R = A.QRdecomposition()
>>> Q
[ ___ ___ ___ ]
[\/ 6 \/ 3 \/ 2 ]
[  ]
[ 6 3 2 ]
[ ]
[ ___ ___ ___ ]
[\/ 6 \/ 3 \/ 2 ]
[   ]
[ 6 3 2 ]
[ ]
[ ___ ___ ]
[\/ 6 \/ 3 ]
[  0 ]
[ 3 3 ]
>>> R
[ ___ ]
[ ___ 4*\/ 6 ___]
[\/ 6  2*\/ 6 ]
[ 3 ]
[ ]
[ ___ ]
[ \/ 3 ]
[ 0  0 ]
[ 3 ]
[ ]
[ ___ ]
[ 0 0 \/ 2 ]
>>> Q*R
[1 1 1]
[ ]
[1 1 3]
[ ]
[2 3 4]
In addition to the solvers in the solver.py
file, we can solve the system Ax=b
by passing the b vector to the matrix A’s LUsolve function. Here we’ll cheat a
little choose A and x then multiply to get b. Then we can solve for x and check
that it’s correct:
>>> A = Matrix([ [2, 3, 5], [3, 6, 2], [8, 3, 6] ])
>>> x = Matrix(3,1,[3,7,5])
>>> b = A*x
>>> soln = A.LUsolve(b)
>>> soln
[3]
[ ]
[7]
[ ]
[5]
There’s also a nice GramSchmidt orthogonalizer which will take a set of
vectors and orthogonalize them with respect to another. There is an
optional argument which specifies whether or not the output should also be
normalized, it defaults to False
. Let’s take some vectors and orthogonalize
them  one normalized and one not:
>>> L = [Matrix([2,3,5]), Matrix([3,6,2]), Matrix([8,3,6])]
>>> out1 = GramSchmidt(L)
>>> out2 = GramSchmidt(L, True)
Let’s take a look at the vectors:
>>> for i in out1:
... print(i)
...
Matrix([[2], [3], [5]])
Matrix([[23/19], [63/19], [47/19]])
Matrix([[1692/353], [1551/706], [423/706]])
>>> for i in out2:
... print(i)
...
Matrix([[sqrt(38)/19], [3*sqrt(38)/38], [5*sqrt(38)/38]])
Matrix([[23*sqrt(6707)/6707], [63*sqrt(6707)/6707], [47*sqrt(6707)/6707]])
Matrix([[12*sqrt(706)/353], [11*sqrt(706)/706], [3*sqrt(706)/706]])
We can spotcheck their orthogonality with dot() and their normality with norm():
>>> out1[0].dot(out1[1])
0
>>> out1[0].dot(out1[2])
0
>>> out1[1].dot(out1[2])
0
>>> out2[0].norm()
1
>>> out2[1].norm()
1
>>> out2[2].norm()
1
So there is quite a bit that can be done with the module including eigenvalues,
eigenvectors, nullspace calculation, cofactor expansion tools, and so on. From
here one might want to look over the matrices.py
file for all functionality.
MatrixBase Class Reference¶

class
sympy.matrices.matrices.
MatrixBase
[source]¶ Base class for matrix objects.
Attributes
cols rows 
D
¶ Return Dirac conjugate (if self.rows == 4).
See also
conjugate
 Byelement conjugation
H
 Hermite conjugation
Examples
>>> from sympy import Matrix, I, eye >>> m = Matrix((0, 1 + I, 2, 3)) >>> m.D Matrix([[0, 1  I, 2, 3]]) >>> m = (eye(4) + I*eye(4)) >>> m[0, 3] = 2 >>> m.D Matrix([ [1  I, 0, 0, 0], [ 0, 1  I, 0, 0], [ 0, 0, 1 + I, 0], [ 2, 0, 0, 1 + I]])
If the matrix does not have 4 rows an AttributeError will be raised because this property is only defined for matrices with 4 rows.
>>> Matrix(eye(2)).D Traceback (most recent call last): ... AttributeError: Matrix has no attribute D.

LDLdecomposition
()[source]¶ Returns the LDL Decomposition (L, D) of matrix A, such that L * D * L.T == A This method eliminates the use of square root. Further this ensures that all the diagonal entries of L are 1. A must be a square, symmetric, positivedefinite and nonsingular matrix.
See also
Examples
>>> from sympy.matrices import Matrix, eye >>> A = Matrix(((25, 15, 5), (15, 18, 0), (5, 0, 11))) >>> L, D = A.LDLdecomposition() >>> L Matrix([ [ 1, 0, 0], [ 3/5, 1, 0], [1/5, 1/3, 1]]) >>> D Matrix([ [25, 0, 0], [ 0, 9, 0], [ 0, 0, 9]]) >>> L * D * L.T * A.inv() == eye(A.rows) True

LDLsolve
(rhs)[source]¶ Solves Ax = B using LDL decomposition, for a general square and nonsingular matrix.
For a nonsquare matrix with rows > cols, the least squares solution is returned.
See also
LDLdecomposition
,lower_triangular_solve
,upper_triangular_solve
,gauss_jordan_solve
,cholesky_solve
,diagonal_solve
,LUsolve
,QRsolve
,pinv_solve
Examples
>>> from sympy.matrices import Matrix, eye >>> A = eye(2)*2 >>> B = Matrix([[1, 2], [3, 4]]) >>> A.LDLsolve(B) == B/2 True

LUdecomposition
(iszerofunc=<function _iszero>, simpfunc=None, rankcheck=False)[source]¶ Returns (L, U, perm) where L is a lower triangular matrix with unit diagonal, U is an upper triangular matrix, and perm is a list of row swap index pairs. If A is the original matrix, then A = (L*U).permuteBkwd(perm), and the row permutation matrix P such that P*A = L*U can be computed by P=eye(A.row).permuteFwd(perm).
See documentation for LUCombined for details about the keyword argument rankcheck, iszerofunc, and simpfunc.
See also
cholesky
,LDLdecomposition
,QRdecomposition
,LUdecomposition_Simple
,LUdecompositionFF
,LUsolve
Examples
>>> from sympy import Matrix >>> a = Matrix([[4, 3], [6, 3]]) >>> L, U, _ = a.LUdecomposition() >>> L Matrix([ [ 1, 0], [3/2, 1]]) >>> U Matrix([ [4, 3], [0, 3/2]])

LUdecompositionFF
()[source]¶ Compute a fractionfree LU decomposition.
Returns 4 matrices P, L, D, U such that PA = L D**1 U. If the elements of the matrix belong to some integral domain I, then all elements of L, D and U are guaranteed to belong to I.
 Reference
 W. Zhou & D.J. Jeffrey, “Fractionfree matrix factors: new forms for LU and QR factors”. Frontiers in Computer Science in China, Vol 2, no. 1, pp. 6780, 2008.
See also

LUdecomposition_Simple
(iszerofunc=<function _iszero>, simpfunc=None, rankcheck=False)[source]¶ Compute an lu decomposition of m x n matrix A, where P*A = L*U
 L is m x m lower triangular with unit diagonal
 U is m x n upper triangular
 P is an m x m permutation matrix
Returns an m x n matrix lu, and an m element list perm where each element of perm is a pair of row exchange indices.
The factors L and U are stored in lu as follows: The subdiagonal elements of L are stored in the subdiagonal elements of lu, that is lu[i, j] = L[i, j] whenever i > j. The elements on the diagonal of L are all 1, and are not explicitly stored. U is stored in the upper triangular portion of lu, that is lu[i ,j] = U[i, j] whenever i <= j. The output matrix can be visualized as:
 Matrix([
 [u, u, u, u], [l, u, u, u], [l, l, u, u], [l, l, l, u]])
where l represents a subdiagonal entry of the L factor, and u represents an entry from the upper triangular entry of the U factor.
perm is a list row swap index pairs such that if A is the original matrix, then A = (L*U).permuteBkwd(perm), and the row permutation matrix P such that P*A = L*U can be computed by soP=eye(A.row).permuteFwd(perm).
The keyword argument rankcheck determines if this function raises a ValueError when passed a matrix whose rank is strictly less than min(num rows, num cols). The default behavior is to decompose a rank deficient matrix. Pass rankcheck=True to raise a ValueError instead. (This mimics the previous behavior of this function).
The keyword arguments iszerofunc and simpfunc are used by the pivot search algorithm. iszerofunc is a callable that returns a boolean indicating if its input is zero, or None if it cannot make the determination. simpfunc is a callable that simplifies its input. The default is simpfunc=None, which indicate that the pivot search algorithm should not attempt to simplify any candidate pivots. If simpfunc fails to simplify its input, then it must return its input instead of a copy.
When a matrix contains symbolic entries, the pivot search algorithm differs from the case where every entry can be categorized as zero or nonzero. The algorithm searches column by column through the submatrix whose top left entry coincides with the pivot position. If it exists, the pivot is the first entry in the current search column that iszerofunc guarantees is nonzero. If no such candidate exists, then each candidate pivot is simplified if simpfunc is not None. The search is repeated, with the difference that a candidate may be the pivot if \(iszerofunc()\) cannot guarantee that it is nonzero. In the second search the pivot is the first candidate that iszerofunc can guarantee is nonzero. If no such candidate exists, then the pivot is the first candidate for which iszerofunc returns None. If no such candidate exists, then the search is repeated in the next column to the right. The pivot search algorithm differs from the one in \(rref()\), which relies on \(_find_reasonable_pivot()\). Future versions of \(LUdecomposition_simple()\) may use \(_find_reasonable_pivot()\).
See also

LUsolve
(rhs, iszerofunc=<function _iszero>)[source]¶ Solve the linear system Ax = rhs for x where A = self.
This is for symbolic matrices, for real or complex ones use mpmath.lu_solve or mpmath.qr_solve.

QRdecomposition
()[source]¶ Return Q, R where A = Q*R, Q is orthogonal and R is upper triangular.
See also
Examples
This is the example from wikipedia:
>>> from sympy import Matrix >>> A = Matrix([[12, 51, 4], [6, 167, 68], [4, 24, 41]]) >>> Q, R = A.QRdecomposition() >>> Q Matrix([ [ 6/7, 69/175, 58/175], [ 3/7, 158/175, 6/175], [2/7, 6/35, 33/35]]) >>> R Matrix([ [14, 21, 14], [ 0, 175, 70], [ 0, 0, 35]]) >>> A == Q*R True
QR factorization of an identity matrix:
>>> A = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) >>> Q, R = A.QRdecomposition() >>> Q Matrix([ [1, 0, 0], [0, 1, 0], [0, 0, 1]]) >>> R Matrix([ [1, 0, 0], [0, 1, 0], [0, 0, 1]])

QRsolve
(b)[source]¶ Solve the linear system ‘Ax = b’.
‘self’ is the matrix ‘A’, the method argument is the vector ‘b’. The method returns the solution vector ‘x’. If ‘b’ is a matrix, the system is solved for each column of ‘b’ and the return value is a matrix of the same shape as ‘b’.
This method is slower (approximately by a factor of 2) but more stable for floatingpoint arithmetic than the LUsolve method. However, LUsolve usually uses an exact arithmetic, so you don’t need to use QRsolve.
This is mainly for educational purposes and symbolic matrices, for real (or complex) matrices use mpmath.qr_solve.

cholesky
()[source]¶ Returns the Cholesky decomposition L of a matrix A such that L * L.T = A
A must be a square, symmetric, positivedefinite and nonsingular matrix.
See also
Examples
>>> from sympy.matrices import Matrix >>> A = Matrix(((25, 15, 5), (15, 18, 0), (5, 0, 11))) >>> A.cholesky() Matrix([ [ 5, 0, 0], [ 3, 3, 0], [1, 1, 3]]) >>> A.cholesky() * A.cholesky().T Matrix([ [25, 15, 5], [15, 18, 0], [5, 0, 11]])

cholesky_solve
(rhs)[source]¶ Solves Ax = B using Cholesky decomposition, for a general square nonsingular matrix. For a nonsquare matrix with rows > cols, the least squares solution is returned.

condition_number
()[source]¶ Returns the condition number of a matrix.
This is the maximum singular value divided by the minimum singular value
See also
singular_values
Examples
>>> from sympy import Matrix, S >>> A = Matrix([[1, 0, 0], [0, 10, 0], [0, 0, S.One/10]]) >>> A.condition_number() 100

copy
()[source]¶ Returns the copy of a matrix.
Examples
>>> from sympy import Matrix >>> A = Matrix(2, 2, [1, 2, 3, 4]) >>> A.copy() Matrix([ [1, 2], [3, 4]])

cross
(b)[source]¶ Return the cross product of \(self\) and \(b\) relaxing the condition of compatible dimensions: if each has 3 elements, a matrix of the same type and shape as \(self\) will be returned. If \(b\) has the same shape as \(self\) then common identities for the cross product (like \(a x b =  b x a\)) will hold.

diagonal_solve
(rhs)[source]¶ Solves Ax = B efficiently, where A is a diagonal Matrix, with nonzero diagonal entries.
See also
lower_triangular_solve
,upper_triangular_solve
,gauss_jordan_solve
,cholesky_solve
,LDLsolve
,LUsolve
,QRsolve
,pinv_solve
Examples
>>> from sympy.matrices import Matrix, eye >>> A = eye(2)*2 >>> B = Matrix([[1, 2], [3, 4]]) >>> A.diagonal_solve(B) == B/2 True

dot
(b)[source]¶ Return the dot product of Matrix self and b relaxing the condition of compatible dimensions: if either the number of rows or columns are the same as the length of b then the dot product is returned. If self is a row or column vector, a scalar is returned. Otherwise, a list of results is returned (and in that case the number of columns in self must match the length of b).
Examples
>>> from sympy import Matrix >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> v = [1, 1, 1] >>> M.row(0).dot(v) 6 >>> M.col(0).dot(v) 12 >>> M.dot(v) [6, 15, 24]

dual
()[source]¶ Returns the dual of a matrix, which is:
\((1/2)*levicivita(i, j, k, l)*M(k, l)\) summed over indices \(k\) and \(l\)
Since the levicivita method is anti_symmetric for any pairwise exchange of indices, the dual of a symmetric matrix is the zero matrix. Strictly speaking the dual defined here assumes that the ‘matrix’ \(M\) is a contravariant anti_symmetric second rank tensor, so that the dual is a covariant second rank tensor.

gauss_jordan_solve
(b, freevar=False)[source]¶ Solves Ax = b using Gauss Jordan elimination.
There may be zero, one, or infinite solutions. If one solution exists, it will be returned. If infinite solutions exist, it will be returned parametrically. If no solutions exist, It will throw ValueError.
Parameters: b : Matrix
The right hand side of the equation to be solved for. Must have the same number of rows as matrix A.
freevar : List
If the system is underdetermined (e.g. A has more columns than rows), infinite solutions are possible, in terms of an arbitrary values of free variables. Then the index of the free variables in the solutions (column Matrix) will be returned by freevar, if the flag \(freevar\) is set to \(True\).
Returns: x : Matrix
The matrix that will satisfy Ax = B. Will have as many rows as matrix A has columns, and as many columns as matrix B.
params : Matrix
If the system is underdetermined (e.g. A has more columns than rows), infinite solutions are possible, in terms of an arbitrary parameters. These arbitrary parameters are returned as params Matrix.
See also
lower_triangular_solve
,upper_triangular_solve
,cholesky_solve
,diagonal_solve
,LDLsolve
,LUsolve
,QRsolve
,pinv
References
[R391] http://en.wikipedia.org/wiki/Gaussian_elimination Examples
>>> from sympy import Matrix >>> A = Matrix([[1, 2, 1, 1], [1, 2, 2, 1], [2, 4, 0, 6]]) >>> b = Matrix([7, 12, 4]) >>> sol, params = A.gauss_jordan_solve(b) >>> sol Matrix([ [2*_tau0  3*_tau1 + 2], [ _tau0], [ 2*_tau1 + 5], [ _tau1]]) >>> params Matrix([ [_tau0], [_tau1]])
>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]]) >>> b = Matrix([3, 6, 9]) >>> sol, params = A.gauss_jordan_solve(b) >>> sol Matrix([ [1], [ 2], [ 0]]) >>> params Matrix(0, 1, [])

inv
(method=None, **kwargs)[source]¶ Return the inverse of a matrix.
CASE 1: If the matrix is a dense matrix.
Return the matrix inverse using the method indicated (default is Gauss elimination).
Parameters: method : (‘GE’, ‘LU’, or ‘ADJ’)
Raises: ValueError
If the determinant of the matrix is zero.
See also
Notes
According to the
method
keyword, it calls the appropriate method:LDL ... inverse_LDL(); default CH .... inverse_CH()Kwargs
method : (‘CH’, ‘LDL’)

inv_mod
(m)[source]¶ Returns the inverse of the matrix \(K\) (mod \(m\)), if it exists.
Method to find the matrix inverse of \(K\) (mod \(m\)) implemented in this function:
 Compute \(\mathrm{adj}(K) = \mathrm{cof}(K)^t\), the adjoint matrix of \(K\).
 Compute \(r = 1/\mathrm{det}(K) \pmod m\).
 \(K^{1} = r\cdot \mathrm{adj}(K) \pmod m\).
Examples
>>> from sympy import Matrix >>> A = Matrix(2, 2, [1, 2, 3, 4]) >>> A.inv_mod(5) Matrix([ [3, 1], [4, 2]]) >>> A.inv_mod(3) Matrix([ [1, 1], [0, 1]])

inverse_ADJ
(iszerofunc=<function _iszero>)[source]¶ Calculates the inverse using the adjugate matrix and a determinant.
See also

inverse_GE
(iszerofunc=<function _iszero>)[source]¶ Calculates the inverse using Gaussian elimination.
See also

inverse_LU
(iszerofunc=<function _iszero>)[source]¶ Calculates the inverse using LU decomposition.
See also

is_nilpotent
()[source]¶ Checks if a matrix is nilpotent.
A matrix B is nilpotent if for some integer k, B**k is a zero matrix.
Examples
>>> from sympy import Matrix >>> a = Matrix([[0, 0, 0], [1, 0, 0], [1, 1, 0]]) >>> a.is_nilpotent() True
>>> a = Matrix([[1, 0, 1], [1, 0, 0], [1, 1, 0]]) >>> a.is_nilpotent() False

key2bounds
(keys)[source]¶ Converts a key with potentially mixed types of keys (integer and slice) into a tuple of ranges and raises an error if any index is out of self’s range.
See also

key2ij
(key)[source]¶ Converts key into canonical form, converting integers or indexable items into valid integers for self’s range or returning slices unchanged.
See also

norm
(ord=None)[source]¶ Return the Norm of a Matrix or Vector. In the simplest case this is the geometric size of the vector Other norms can be specified by the ord parameter
ord norm for matrices norm for vectors None Frobenius norm 2norm ‘fro’ Frobenius norm  does not exist
inf – max(abs(x)) inf – min(abs(x)) 1 – as below 1 – as below 2 2norm (largest sing. value) as below 2 smallest singular value as below other  does not exist
sum(abs(x)**ord)**(1./ord) See also
Examples
>>> from sympy import Matrix, Symbol, trigsimp, cos, sin, oo >>> x = Symbol('x', real=True) >>> v = Matrix([cos(x), sin(x)]) >>> trigsimp( v.norm() ) 1 >>> v.norm(10) (sin(x)**10 + cos(x)**10)**(1/10) >>> A = Matrix([[1, 1], [1, 1]]) >>> A.norm(2)# Spectral norm (max of Ax/x under 2vectornorm) 2 >>> A.norm(2) # Inverse spectral norm (smallest singular value) 0 >>> A.norm() # Frobenius Norm 2 >>> Matrix([1, 2]).norm(oo) 2 >>> Matrix([1, 2]).norm(oo) 1

pinv
()[source]¶ Calculate the MoorePenrose pseudoinverse of the matrix.
The MoorePenrose pseudoinverse exists and is unique for any matrix. If the matrix is invertible, the pseudoinverse is the same as the inverse.
See also
References
[R392] https://en.wikipedia.org/wiki/MoorePenrose_pseudoinverse Examples
>>> from sympy import Matrix >>> Matrix([[1, 2, 3], [4, 5, 6]]).pinv() Matrix([ [17/18, 4/9], [ 1/9, 1/9], [ 13/18, 2/9]])

pinv_solve
(B, arbitrary_matrix=None)[source]¶ Solve Ax = B using the MoorePenrose pseudoinverse.
There may be zero, one, or infinite solutions. If one solution exists, it will be returned. If infinite solutions exist, one will be returned based on the value of arbitrary_matrix. If no solutions exist, the leastsquares solution is returned.
Parameters: B : Matrix
The right hand side of the equation to be solved for. Must have the same number of rows as matrix A.
arbitrary_matrix : Matrix
If the system is underdetermined (e.g. A has more columns than rows), infinite solutions are possible, in terms of an arbitrary matrix. This parameter may be set to a specific matrix to use for that purpose; if so, it must be the same shape as x, with as many rows as matrix A has columns, and as many columns as matrix B. If left as None, an appropriate matrix containing dummy symbols in the form of
wn_m
will be used, with n and m being row and column position of each symbol.Returns: x : Matrix
The matrix that will satisfy Ax = B. Will have as many rows as matrix A has columns, and as many columns as matrix B.
See also
lower_triangular_solve
,upper_triangular_solve
,gauss_jordan_solve
,cholesky_solve
,diagonal_solve
,LDLsolve
,LUsolve
,QRsolve
,pinv
Notes
This may return either exact solutions or least squares solutions. To determine which, check
A * A.pinv() * B == B
. It will be True if exact solutions exist, and False if only a leastsquares solution exists. Be aware that the left hand side of that equation may need to be simplified to correctly compare to the right hand side.References
[R393] https://en.wikipedia.org/wiki/MoorePenrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system Examples
>>> from sympy import Matrix >>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> B = Matrix([7, 8]) >>> A.pinv_solve(B) Matrix([ [ _w0_0/6  _w1_0/3 + _w2_0/6  55/18], [_w0_0/3 + 2*_w1_0/3  _w2_0/3 + 1/9], [ _w0_0/6  _w1_0/3 + _w2_0/6 + 59/18]]) >>> A.pinv_solve(B, arbitrary_matrix=Matrix([0, 0, 0])) Matrix([ [55/18], [ 1/9], [ 59/18]])

print_nonzero
(symb='X')[source]¶ Shows location of nonzero entries for fast shape lookup.
Examples
>>> from sympy.matrices import Matrix, eye >>> m = Matrix(2, 3, lambda i, j: i*3+j) >>> m Matrix([ [0, 1, 2], [3, 4, 5]]) >>> m.print_nonzero() [ XX] [XXX] >>> m = eye(4) >>> m.print_nonzero("x") [x ] [ x ] [ x ] [ x]

project
(v)[source]¶ Return the projection of
self
onto the line containingv
.Examples
>>> from sympy import Matrix, S, sqrt >>> V = Matrix([sqrt(3)/2, S.Half]) >>> x = Matrix([[1, 0]]) >>> V.project(x) Matrix([[sqrt(3)/2, 0]]) >>> V.project(x) Matrix([[sqrt(3)/2, 0]])

solve
(rhs, method='GE')[source]¶ Return solution to self*soln = rhs using given inversion method.
For a list of possible inversion methods, see the .inv() docstring.

solve_least_squares
(rhs, method='CH')[source]¶ Return the leastsquare fit to the data.
By default the cholesky_solve routine is used (method=’CH’); other methods of matrix inversion can be used. To find out which are available, see the docstring of the .inv() method.
Examples
>>> from sympy.matrices import Matrix, ones >>> A = Matrix([1, 2, 3]) >>> B = Matrix([2, 3, 4]) >>> S = Matrix(A.row_join(B)) >>> S Matrix([ [1, 2], [2, 3], [3, 4]])
If each line of S represent coefficients of Ax + By and x and y are [2, 3] then S*xy is:
>>> r = S*Matrix([2, 3]); r Matrix([ [ 8], [13], [18]])
But let’s add 1 to the middle value and then solve for the leastsquares value of xy:
>>> xy = S.solve_least_squares(Matrix([8, 14, 18])); xy Matrix([ [ 5/3], [10/3]])
The error is given by S*xy  r:
>>> S*xy  r Matrix([ [1/3], [1/3], [1/3]]) >>> _.norm().n(2) 0.58
If a different xy is used, the norm will be higher:
>>> xy += ones(2, 1)/10 >>> (S*xy  r).norm().n(2) 1.5

table
(printer, rowstart='[', rowend=']', rowsep='\n', colsep=', ', align='right')[source]¶ String form of Matrix as a table.
printer
is the printer to use for on the elements (generally something like StrPrinter())rowstart
is the string used to start each row (by default ‘[‘).rowend
is the string used to end each row (by default ‘]’).rowsep
is the string used to separate rows (by default a newline).colsep
is the string used to separate columns (by default ‘, ‘).align
defines how the elements are aligned. Must be one of ‘left’, ‘right’, or ‘center’. You can also use ‘<’, ‘>’, and ‘^’ to mean the same thing, respectively.This is used by the string printer for Matrix.
Examples
>>> from sympy import Matrix >>> from sympy.printing.str import StrPrinter >>> M = Matrix([[1, 2], [33, 4]]) >>> printer = StrPrinter() >>> M.table(printer) '[ 1, 2]\n[33, 4]' >>> print(M.table(printer)) [ 1, 2] [33, 4] >>> print(M.table(printer, rowsep=',\n')) [ 1, 2], [33, 4] >>> print('[%s]' % M.table(printer, rowsep=',\n')) [[ 1, 2], [33, 4]] >>> print(M.table(printer, colsep=' ')) [ 1 2] [33 4] >>> print(M.table(printer, align='center')) [ 1 , 2] [33, 4] >>> print(M.table(printer, rowstart='{', rowend='}')) { 1, 2} {33, 4}

vech
(diagonal=True, check_symmetry=True)[source]¶ Return the unique elements of a symmetric Matrix as a one column matrix by stacking the elements in the lower triangle.
Arguments: diagonal – include the diagonal cells of self or not check_symmetry – checks symmetry of self but not completely reliably
See also
vec
Examples
>>> from sympy import Matrix >>> m=Matrix([[1, 2], [2, 3]]) >>> m Matrix([ [1, 2], [2, 3]]) >>> m.vech() Matrix([ [1], [2], [3]]) >>> m.vech(diagonal=False) Matrix([[2]])

Matrix Exceptions Reference¶
Matrix Functions Reference¶

sympy.matrices.matrices.
classof
(A, B)[source]¶ Get the type of the result when combining matrices of different types.
Currently the strategy is that immutability is contagious.
Examples
>>> from sympy import Matrix, ImmutableMatrix >>> from sympy.matrices.matrices import classof >>> M = Matrix([[1, 2], [3, 4]]) # a Mutable Matrix >>> IM = ImmutableMatrix([[1, 2], [3, 4]]) >>> classof(M, IM) <class 'sympy.matrices.immutable.ImmutableDenseMatrix'>

sympy.matrices.dense.
matrix_multiply_elementwise
(A, B)[source]¶ Return the Hadamard product (elementwise product) of A and B
>>> from sympy.matrices import matrix_multiply_elementwise >>> from sympy.matrices import Matrix >>> A = Matrix([[0, 1, 2], [3, 4, 5]]) >>> B = Matrix([[1, 10, 100], [100, 10, 1]]) >>> matrix_multiply_elementwise(A, B) Matrix([ [ 0, 10, 200], [300, 40, 5]])
See also
__mul__

sympy.matrices.dense.
zeros
(*args, **kwargs)[source]¶ Returns a matrix of zeros with
rows
rows andcols
columns; ifcols
is omitted a square matrix will be returned.See also
ones
,eye
,diag

sympy.matrices.dense.
ones
(*args, **kwargs)[source]¶ Returns a matrix of ones with
rows
rows andcols
columns; ifcols
is omitted a square matrix will be returned.See also
zeros
,eye
,diag

sympy.matrices.dense.
eye
(*args, **kwargs)[source]¶ Create square identity matrix n x n
See also
diag
,zeros
,ones

sympy.matrices.dense.
diag
(*values, **kwargs)[source]¶ Create a sparse, diagonal matrix from a list of diagonal values.
See also
eye
Notes
When arguments are matrices they are fitted in resultant matrix.
The returned matrix is a mutable, dense matrix. To make it a different type, send the desired class for keyword
cls
.Examples
>>> from sympy.matrices import diag, Matrix, ones >>> diag(1, 2, 3) Matrix([ [1, 0, 0], [0, 2, 0], [0, 0, 3]]) >>> diag(*[1, 2, 3]) Matrix([ [1, 0, 0], [0, 2, 0], [0, 0, 3]])
The diagonal elements can be matrices; diagonal filling will continue on the diagonal from the last element of the matrix:
>>> from sympy.abc import x, y, z >>> a = Matrix([x, y, z]) >>> b = Matrix([[1, 2], [3, 4]]) >>> c = Matrix([[5, 6]]) >>> diag(a, 7, b, c) Matrix([ [x, 0, 0, 0, 0, 0], [y, 0, 0, 0, 0, 0], [z, 0, 0, 0, 0, 0], [0, 7, 0, 0, 0, 0], [0, 0, 1, 2, 0, 0], [0, 0, 3, 4, 0, 0], [0, 0, 0, 0, 5, 6]])
When diagonal elements are lists, they will be treated as arguments to Matrix:
>>> diag([1, 2, 3], 4) Matrix([ [1, 0], [2, 0], [3, 0], [0, 4]]) >>> diag([[1, 2, 3]], 4) Matrix([ [1, 2, 3, 0], [0, 0, 0, 4]])
A given band off the diagonal can be made by padding with a vertical or horizontal “kerning” vector:
>>> hpad = ones(0, 2) >>> vpad = ones(2, 0) >>> diag(vpad, 1, 2, 3, hpad) + diag(hpad, 4, 5, 6, vpad) Matrix([ [0, 0, 4, 0, 0], [0, 0, 0, 5, 0], [1, 0, 0, 0, 6], [0, 2, 0, 0, 0], [0, 0, 3, 0, 0]])
The type is mutable by default but can be made immutable by setting the
mutable
flag to False:>>> type(diag(1)) <class 'sympy.matrices.dense.MutableDenseMatrix'> >>> from sympy.matrices import ImmutableMatrix >>> type(diag(1, cls=ImmutableMatrix)) <class 'sympy.matrices.immutable.ImmutableDenseMatrix'>

sympy.matrices.dense.
jordan_cell
(eigenval, n)[source]¶ Create a Jordan block:
Examples
>>> from sympy.matrices import jordan_cell >>> from sympy.abc import x >>> jordan_cell(x, 4) Matrix([ [x, 1, 0, 0], [0, x, 1, 0], [0, 0, x, 1], [0, 0, 0, x]])

sympy.matrices.dense.
hessian
(f, varlist, constraints=[])[source]¶ Compute Hessian matrix for a function f wrt parameters in varlist which may be given as a sequence or a row/column vector. A list of constraints may optionally be given.
See also
sympy.matrices.mutable.Matrix.jacobian
,wronskian
References
http://en.wikipedia.org/wiki/Hessian_matrix
Examples
>>> from sympy import Function, hessian, pprint >>> from sympy.abc import x, y >>> f = Function('f')(x, y) >>> g1 = Function('g')(x, y) >>> g2 = x**2 + 3*y >>> pprint(hessian(f, (x, y), [g1, g2])) [ d d ] [ 0 0 (g(x, y)) (g(x, y)) ] [ dx dy ] [ ] [ 0 0 2*x 3 ] [ ] [ 2 2 ] [d d d ] [(g(x, y)) 2*x (f(x, y)) (f(x, y))] [dx 2 dy dx ] [ dx ] [ ] [ 2 2 ] [d d d ] [(g(x, y)) 3 (f(x, y)) (f(x, y)) ] [dy dy dx 2 ] [ dy ]

sympy.matrices.dense.
GramSchmidt
(vlist, orthonormal=False)[source]¶ Apply the GramSchmidt process to a set of vectors.
see: http://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

sympy.matrices.dense.
wronskian
(functions, var, method='bareiss')[source]¶ Compute Wronskian for [] of functions
 f1 f2 ... fn   f1' f2' ... fn'   . . . .  W(f1, ..., fn) =  . . . .   . . . .   (n) (n) (n)   D (f1) D (f2) ... D (fn) 
see: http://en.wikipedia.org/wiki/Wronskian
See also
sympy.matrices.mutable.Matrix.jacobian
,hessian

sympy.matrices.dense.
casoratian
(seqs, n, zero=True)[source]¶ Given linear difference operator L of order ‘k’ and homogeneous equation Ly = 0 we want to compute kernel of L, which is a set of ‘k’ sequences: a(n), b(n), ... z(n).
Solutions of L are linearly independent iff their Casoratian, denoted as C(a, b, ..., z), do not vanish for n = 0.
Casoratian is defined by k x k determinant:
+ a(n) b(n) . . . z(n) +  a(n+1) b(n+1) . . . z(n+1)   . . . .   . . . .   . . . .  + a(n+k1) b(n+k1) . . . z(n+k1) +
It proves very useful in rsolve_hyper() where it is applied to a generating set of a recurrence to factor out linearly dependent solutions and return a basis:
>>> from sympy import Symbol, casoratian, factorial >>> n = Symbol('n', integer=True)
Exponential and factorial are linearly independent:
>>> casoratian([2**n, factorial(n)], n) != 0 True

sympy.matrices.dense.
randMatrix
(r, c=None, min=0, max=99, seed=None, symmetric=False, percent=100, prng=None)[source]¶ Create random matrix with dimensions
r
xc
. Ifc
is omitted the matrix will be square. Ifsymmetric
is True the matrix must be square. Ifpercent
is less than 100 then only approximately the given percentage of elements will be nonzero.The pseudorandom number generator used to generate matrix is chosen in the following way.
 If
prng
is supplied, it will be used as random number generator. It should be an instance ofrandom.Random
, or at least haverandint
andshuffle
methods with same signatures.  if
prng
is not supplied butseed
is supplied, then newrandom.Random
with givenseed
will be created;  otherwise, a new
random.Random
with default seed will be used.
Examples
>>> from sympy.matrices import randMatrix >>> randMatrix(3) [25, 45, 27] [44, 54, 9] [23, 96, 46] >>> randMatrix(3, 2) [87, 29] [23, 37] [90, 26] >>> randMatrix(3, 3, 0, 2) [0, 2, 0] [2, 0, 1] [0, 0, 1] >>> randMatrix(3, symmetric=True) [85, 26, 29] [26, 71, 43] [29, 43, 57] >>> A = randMatrix(3, seed=1) >>> B = randMatrix(3, seed=2) >>> A == B False >>> A == randMatrix(3, seed=1) True >>> randMatrix(3, symmetric=True, percent=50) [0, 68, 43] [0, 68, 0] [0, 91, 34]
 If
Numpy Utility Functions Reference¶

sympy.matrices.dense.
list2numpy
(l, dtype=<class 'object'>)[source]¶ Converts python list of SymPy expressions to a NumPy array.
See also
matrix2numpy

sympy.matrices.dense.
matrix2numpy
(m, dtype=<class 'object'>)[source]¶ Converts SymPy’s matrix to a NumPy array.
See also
list2numpy

sympy.matrices.dense.
symarray
(prefix, shape, **kwargs)[source]¶ Create a numpy ndarray of symbols (as an object array).
The created symbols are named
prefix_i1_i2_
... You should thus provide a nonempty prefix if you want your symbols to be unique for different output arrays, as SymPy symbols with identical names are the same object.Parameters: prefix : string
A prefix prepended to the name of every symbol.
shape : int or tuple
Shape of the created array. If an int, the array is onedimensional; for more than one dimension the shape must be a tuple.
**kwargs : dict
keyword arguments passed on to Symbol
Examples
These doctests require numpy.
>>> from sympy import symarray >>> symarray('', 3) [_0 _1 _2]
If you want multiple symarrays to contain distinct symbols, you must provide unique prefixes:
>>> a = symarray('', 3) >>> b = symarray('', 3) >>> a[0] == b[0] True >>> a = symarray('a', 3) >>> b = symarray('b', 3) >>> a[0] == b[0] False
Creating symarrays with a prefix:
>>> symarray('a', 3) [a_0 a_1 a_2]
For more than one dimension, the shape must be given as a tuple:
>>> symarray('a', (2, 3)) [[a_0_0 a_0_1 a_0_2] [a_1_0 a_1_1 a_1_2]] >>> symarray('a', (2, 3, 2)) [[[a_0_0_0 a_0_0_1] [a_0_1_0 a_0_1_1] [a_0_2_0 a_0_2_1]] [[a_1_0_0 a_1_0_1] [a_1_1_0 a_1_1_1] [a_1_2_0 a_1_2_1]]]
For setting assumptions of the underlying Symbols:
>>> [s.is_real for s in symarray('a', 2, real=True)] [True, True]

sympy.matrices.dense.
rot_axis1
(theta)[source]¶ Returns a rotation matrix for a rotation of theta (in radians) about the 1axis.
See also
rot_axis2
 Returns a rotation matrix for a rotation of theta (in radians) about the 2axis
rot_axis3
 Returns a rotation matrix for a rotation of theta (in radians) about the 3axis
Examples
>>> from sympy import pi >>> from sympy.matrices import rot_axis1
A rotation of pi/3 (60 degrees):
>>> theta = pi/3 >>> rot_axis1(theta) Matrix([ [1, 0, 0], [0, 1/2, sqrt(3)/2], [0, sqrt(3)/2, 1/2]])
If we rotate by pi/2 (90 degrees):
>>> rot_axis1(pi/2) Matrix([ [1, 0, 0], [0, 0, 1], [0, 1, 0]])

sympy.matrices.dense.
rot_axis2
(theta)[source]¶ Returns a rotation matrix for a rotation of theta (in radians) about the 2axis.
See also
rot_axis1
 Returns a rotation matrix for a rotation of theta (in radians) about the 1axis
rot_axis3
 Returns a rotation matrix for a rotation of theta (in radians) about the 3axis
Examples
>>> from sympy import pi >>> from sympy.matrices import rot_axis2
A rotation of pi/3 (60 degrees):
>>> theta = pi/3 >>> rot_axis2(theta) Matrix([ [ 1/2, 0, sqrt(3)/2], [ 0, 1, 0], [sqrt(3)/2, 0, 1/2]])
If we rotate by pi/2 (90 degrees):
>>> rot_axis2(pi/2) Matrix([ [0, 0, 1], [0, 1, 0], [1, 0, 0]])

sympy.matrices.dense.
rot_axis3
(theta)[source]¶ Returns a rotation matrix for a rotation of theta (in radians) about the 3axis.
See also
rot_axis1
 Returns a rotation matrix for a rotation of theta (in radians) about the 1axis
rot_axis2
 Returns a rotation matrix for a rotation of theta (in radians) about the 2axis
Examples
>>> from sympy import pi >>> from sympy.matrices import rot_axis3
A rotation of pi/3 (60 degrees):
>>> theta = pi/3 >>> rot_axis3(theta) Matrix([ [ 1/2, sqrt(3)/2, 0], [sqrt(3)/2, 1/2, 0], [ 0, 0, 1]])
If we rotate by pi/2 (90 degrees):
>>> rot_axis3(pi/2) Matrix([ [ 0, 1, 0], [1, 0, 0], [ 0, 0, 1]])