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)
>>> 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 appropriately-sized 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 2-variable 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 1-line 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 2-dimensional, the storage is not and so it is allowable - though one should be careful - to access the entries as if they were a 1-d 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 0-based 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 10-th column (not there)
Traceback (most recent call last):
...
IndexError: Index out of range: a[[0, 10]]
>>> M[:, 10:11] # the 10-th column (if there)
[]
>>> M[:, :10] # all columns up to the 10-th
[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 matrix-wide 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 Gram-Schmidt 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 spot-check 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.

MatrixDeterminant Class Reference

class sympy.matrices.matrices.MatrixDeterminant

Provides basic matrix determinant operations. Should not be instantiated directly.

adjugate(method='berkowitz')

Returns the adjugate, or classical adjoint, of a matrix. That is, the transpose of the matrix of cofactors.

https://en.wikipedia.org/wiki/Adjugate

See also

cofactor_matrix, transpose

charpoly(x='lambda', simplify=<function simplify>)

Computes characteristic polynomial det(x*I - self) where I is the identity matrix.

A PurePoly is returned, so using different variables for x does not affect the comparison or the polynomials:

Examples

>>> from sympy import Matrix
>>> from sympy.abc import x, y
>>> A = Matrix([[1, 3], [2, 0]])
>>> A.charpoly(x) == A.charpoly(y)
True

Specifying x is optional; a symbol named lambda is used by default (which looks good when pretty-printed in unicode):

>>> A.charpoly().as_expr()
lambda**2 - lambda - 6

And if x clashes with an existing symbol, underscores will be preppended to the name to make it unique:

>>> A = Matrix([[1, 2], [x, 0]])
>>> A.charpoly(x).as_expr()
_x**2 - _x - 2*x

Whether you pass a symbol or not, the generator can be obtained with the gen attribute since it may not be the same as the symbol that was passed:

>>> A.charpoly(x).gen
_x
>>> A.charpoly(x).gen == x
False

Notes

The Samuelson-Berkowitz algorithm is used to compute the characteristic polynomial efficiently and without any division operations. Thus the characteristic polynomial over any commutative ring without zero divisors can be computed.

See also

det

cofactor(i, j, method='berkowitz')

Calculate the cofactor of an element.

See also

cofactor_matrix, minor, minor_submatrix

cofactor_matrix(method='berkowitz')

Return a matrix containing the cofactor of each element.

See also

cofactor, minor, minor_submatrix, adjugate

det(method='bareiss', iszerofunc=None)

Computes the determinant of a matrix.

Parameters

method : string, optional

Specifies the algorithm used for computing the matrix determinant.

If the matrix is at most 3x3, a hard-coded formula is used and the specified method is ignored. Otherwise, it defaults to 'bareiss'.

If it is set to 'bareiss', Bareiss’ fraction-free algorithm will be used.

If it is set to 'berkowitz', Berkowitz’ algorithm will be used.

Otherwise, if it is set to 'lu', LU decomposition will be used.

Note

For backward compatibility, legacy keys like “bareis” and “det_lu” can still be used to indicate the corresponding methods. And the keys are also case-insensitive for now. However, it is suggested to use the precise keys for specifying the method.

iszerofunc : FunctionType or None, optional

If it is set to None, it will be defaulted to _iszero if the method is set to 'bareiss', and _is_zero_after_expand_mul if the method is set to 'lu'.

It can also accept any user-specified zero testing function, if it is formatted as a function which accepts a single symbolic argument and returns True if it is tested as zero and False if it tested as non-zero, and also None if it is undecidable.

Returns

det : Basic

Result of determinant.

Raises

ValueError

If unrecognized keys are given for method or iszerofunc.

NonSquareMatrixError

If attempted to calculate determinant from a non-square matrix.

minor(i, j, method='berkowitz')

Return the (i,j) minor of self. That is, return the determinant of the matrix obtained by deleting the \(i\).

See also

minor_submatrix, cofactor, det

minor_submatrix(i, j)

Return the submatrix obtained by removing the \(i\).

See also

minor, cofactor

MatrixReductions Class Reference

class sympy.matrices.matrices.MatrixReductions

Provides basic matrix row/column operations. Should not be instantiated directly.

echelon_form(iszerofunc=<function _iszero>, simplify=False, with_pivots=False)

Returns a matrix row-equivalent to self that is in echelon form. Note that echelon form of a matrix is not unique, however, properties like the row space and the null space are preserved.

elementary_col_op(op='n->kn', col=None, k=None, col1=None, col2=None)

Performs the elementary column operation \(op\).

\(op\) may be one of

• “n->kn” (column n goes to k*n)

• “n<->m” (swap column n and column m)

• “n->n+km” (column n goes to column n + k*column m)

Parameters

op : string; the elementary row operation

col : the column to apply the column operation

k : the multiple to apply in the column operation

col1 : one column of a column swap

col2 : second column of a column swap or column “m” in the column operation

“n->n+km”

elementary_row_op(op='n->kn', row=None, k=None, row1=None, row2=None)

Performs the elementary row operation \(op\).

\(op\) may be one of

• “n->kn” (row n goes to k*n)

• “n<->m” (swap row n and row m)

• “n->n+km” (row n goes to row n + k*row m)

Parameters

op : string; the elementary row operation

row : the row to apply the row operation

k : the multiple to apply in the row operation

row1 : one row of a row swap

row2 : second row of a row swap or row “m” in the row operation

“n->n+km”

is_echelon

Returns \(True\) if the matrix is in echelon form. That is, all rows of zeros are at the bottom, and below each leading non-zero in a row are exclusively zeros.

rank(iszerofunc=<function _iszero>, simplify=False)

Returns the rank of a matrix

>>> from sympy import Matrix
>>> from sympy.abc import x
>>> m = Matrix([[1, 2], [x, 1 - 1/x]])
>>> m.rank()
2
>>> n = Matrix(3, 3, range(1, 10))
>>> n.rank()
2

rref(iszerofunc=<function _iszero>, simplify=False, pivots=True, normalize_last=True)

Return reduced row-echelon form of matrix and indices of pivot vars.

Parameters

iszerofunc : Function

A function used for detecting whether an element can act as a pivot. lambda x: x.is_zero is used by default.

simplify : Function

A function used to simplify elements when looking for a pivot. By default SymPy’s simplify is used.

pivots : True or False

If True, a tuple containing the row-reduced matrix and a tuple of pivot columns is returned. If False just the row-reduced matrix is returned.

normalize_last : True or False

If True, no pivots are normalized to \(1\) until after all entries above and below each pivot are zeroed. This means the row reduction algorithm is fraction free until the very last step. If False, the naive row reduction procedure is used where each pivot is normalized to be \(1\) before row operations are used to zero above and below the pivot.

Notes

The default value of normalize_last=True can provide significant speedup to row reduction, especially on matrices with symbols. However, if you depend on the form row reduction algorithm leaves entries of the matrix, set noramlize_last=False

Examples

>>> from sympy import Matrix
>>> from sympy.abc import x
>>> m = Matrix([[1, 2], [x, 1 - 1/x]])
>>> m.rref()
(Matrix([
[1, 0],
[0, 1]]), (0, 1))
>>> rref_matrix, rref_pivots = m.rref()
>>> rref_matrix
Matrix([
[1, 0],
[0, 1]])
>>> rref_pivots
(0, 1)

MatrixSubspaces Class Reference

class sympy.matrices.matrices.MatrixSubspaces

Provides methods relating to the fundamental subspaces of a matrix. Should not be instantiated directly.

columnspace(simplify=False)

Returns a list of vectors (Matrix objects) that span columnspace of self

Examples

>>> from sympy.matrices import Matrix
>>> m = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
>>> m
Matrix([
[ 1,  3, 0],
[-2, -6, 0],
[ 3,  9, 6]])
>>> m.columnspace()
[Matrix([
[ 1],
[-2],
[ 3]]), Matrix([
[0],
[0],
[6]])]

See also

nullspace, rowspace

nullspace(simplify=False, iszerofunc=<function _iszero>)

Returns list of vectors (Matrix objects) that span nullspace of self

Examples

>>> from sympy.matrices import Matrix
>>> m = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
>>> m
Matrix([
[ 1,  3, 0],
[-2, -6, 0],
[ 3,  9, 6]])
>>> m.nullspace()
[Matrix([
[-3],
[ 1],
[ 0]])]

See also

columnspace, rowspace

classmethod orthogonalize(*vecs, **kwargs)

Apply the Gram-Schmidt orthogonalization procedure to vectors supplied in vecs.

Parameters

vecs

vectors to be made orthogonal

normalize : bool

If true, return an orthonormal basis.

rowspace(simplify=False)

Returns a list of vectors that span the row space of self.

MatrixEigen Class Reference

class sympy.matrices.matrices.MatrixEigen

Provides basic matrix eigenvalue/vector operations. Should not be instantiated directly.

diagonalize(reals_only=False, sort=False, normalize=False)

Return (P, D), where D is diagonal and

D = P^-1 * M * P

where M is current matrix.

Parameters

reals_only : bool. Whether to throw an error if complex numbers are need

to diagonalize. (Default: False)

sort : bool. Sort the eigenvalues along the diagonal. (Default: False)

normalize : bool. If True, normalize the columns of P. (Default: False)

Examples

>>> from sympy import Matrix
>>> m = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
>>> m
Matrix([
[1,  2, 0],
[0,  3, 0],
[2, -4, 2]])
>>> (P, D) = m.diagonalize()
>>> D
Matrix([
[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
>>> P
Matrix([
[-1, 0, -1],
[ 0, 0, -1],
[ 2, 1,  2]])
>>> P.inv() * m * P
Matrix([
[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])

See also

is_diagonal, is_diagonalizable

eigenvals(error_when_incomplete=True, **flags)

Return eigenvalues using the Berkowitz agorithm to compute the characteristic polynomial.

Parameters

error_when_incomplete : bool, optional

If it is set to True, it will raise an error if not all eigenvalues are computed. This is caused by roots not returning a full list of eigenvalues.

simplify : bool or function, optional

If it is set to True, it attempts to return the most simplified form of expressions returned by applying default simplification method in every routine.

If it is set to False, it will skip simplification in this particular routine to save computation resources.

If a function is passed to, it will attempt to apply the particular function as simplification method.

rational : bool, optional

If it is set to True, every floating point numbers would be replaced with rationals before computation. It can solve some issues of roots routine not working well with floats.

multiple : bool, optional

If it is set to True, the result will be in the form of a list.

If it is set to False, the result will be in the form of a dictionary.

Returns

eigs : list or dict

Eigenvalues of a matrix. The return format would be specified by the key multiple.

Raises

MatrixError

If not enough roots had got computed.

NonSquareMatrixError

If attempted to compute eigenvalues from a non-square matrix.

Notes

Eigenvalues of a matrix \(A\) can be computed by solving a matrix equation \(\det(A - \lambda I) = 0\)

See also

MatrixDeterminant.charpoly, eigenvects

eigenvects(error_when_incomplete=True, iszerofunc=<function _iszero>, **flags)

Return list of triples (eigenval, multiplicity, eigenspace).

Parameters

error_when_incomplete : bool, optional

Raise an error when not all eigenvalues are computed. This is caused by roots not returning a full list of eigenvalues.

iszerofunc : function, optional

Specifies a zero testing function to be used in rref.

Default value is _iszero, which uses SymPy’s naive and fast default assumption handler.

It can also accept any user-specified zero testing function, if it is formatted as a function which accepts a single symbolic argument and returns True if it is tested as zero and False if it is tested as non-zero, and None if it is undecidable.

simplify : bool or function, optional

If True, as_content_primitive() will be used to tidy up normalization artifacts.

It will also be used by the nullspace routine.

chop : bool or positive number, optional

If the matrix contains any Floats, they will be changed to Rationals for computation purposes, but the answers will be returned after being evaluated with evalf. The chop flag is passed to evalf. When chop=True a default precision will be used; a number will be interpreted as the desired level of precision.

Returns

ret : [(eigenval, multiplicity, eigenspace), …]

A ragged list containing tuples of data obtained by eigenvals and nullspace.

eigenspace is a list containing the eigenvector for each eigenvalue.

eigenvector is a vector in the form of a Matrix. e.g. a vector of length 3 is returned as Matrix([a_1, a_2, a_3]).

Raises

NotImplementedError

If failed to compute nullspace.

See also

eigenvals, MatrixSubspaces.nullspace

is_diagonalizable(reals_only=False, **kwargs)

Returns true if a matrix is diagonalizable.

Parameters

reals_only : bool. If reals_only=True, determine whether the matrix can be

diagonalized without complex numbers. (Default: False)

Kwargs

clear_cachebool. If True, clear the result of any computations when finished.

(Default: True)

Examples

>>> from sympy import Matrix
>>> m = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
>>> m
Matrix([
[1,  2, 0],
[0,  3, 0],
[2, -4, 2]])
>>> m.is_diagonalizable()
True
>>> m = Matrix(2, 2, [0, 1, 0, 0])
>>> m
Matrix([
[0, 1],
[0, 0]])
>>> m.is_diagonalizable()
False
>>> m = Matrix(2, 2, [0, 1, -1, 0])
>>> m
Matrix([
[ 0, 1],
[-1, 0]])
>>> m.is_diagonalizable()
True
>>> m.is_diagonalizable(reals_only=True)
False

See also

is_diagonal, diagonalize

jordan_form(calc_transform=True, **kwargs)

Return (P, J) where \(J\) is a Jordan block matrix and \(P\) is a matrix such that

self == P*J*P**-1

Parameters

calc_transform : bool

If False, then only \(J\) is returned.

chop : bool

All matrices are convered to exact types when computing eigenvalues and eigenvectors. As a result, there may be approximation errors. If chop==True, these errors will be truncated.

Examples

>>> from sympy import Matrix
>>> m = Matrix([[ 6,  5, -2, -3], [-3, -1,  3,  3], [ 2,  1, -2, -3], [-1,  1,  5,  5]])
>>> P, J = m.jordan_form()
>>> J
Matrix([
[2, 1, 0, 0],
[0, 2, 0, 0],
[0, 0, 2, 1],
[0, 0, 0, 2]])

See also

jordan_block

left_eigenvects(**flags)

Returns left eigenvectors and eigenvalues.

This function returns the list of triples (eigenval, multiplicity, basis) for the left eigenvectors. Options are the same as for eigenvects(), i.e. the **flags arguments gets passed directly to eigenvects().

Examples

>>> from sympy import Matrix
>>> M = Matrix([[0, 1, 1], [1, 0, 0], [1, 1, 1]])
>>> M.eigenvects()
[(-1, 1, [Matrix([
[-1],
[ 1],
[ 0]])]), (0, 1, [Matrix([
[ 0],
[-1],
[ 1]])]), (2, 1, [Matrix([
[2/3],
[1/3],
[  1]])])]
>>> M.left_eigenvects()
[(-1, 1, [Matrix([[-2, 1, 1]])]), (0, 1, [Matrix([[-1, -1, 1]])]), (2,
1, [Matrix([[1, 1, 1]])])]

singular_values()

Compute the singular values of a Matrix

Examples

>>> from sympy import Matrix, Symbol
>>> x = Symbol('x', real=True)
>>> A = Matrix([[0, 1, 0], [0, x, 0], [-1, 0, 0]])
>>> A.singular_values()
[sqrt(x**2 + 1), 1, 0]

See also

condition_number

MatrixCalculus Class Reference

class sympy.matrices.matrices.MatrixCalculus

Provides calculus-related matrix operations.

diff(*args, **kwargs)

Calculate the derivative of each element in the matrix. args will be passed to the integrate function.

Examples

>>> from sympy.matrices import Matrix
>>> from sympy.abc import x, y
>>> M = Matrix([[x, y], [1, 0]])
>>> M.diff(x)
Matrix([
[1, 0],
[0, 0]])

See also

integrate, limit

integrate(*args)

Integrate each element of the matrix. args will be passed to the integrate function.

Examples

>>> from sympy.matrices import Matrix
>>> from sympy.abc import x, y
>>> M = Matrix([[x, y], [1, 0]])
>>> M.integrate((x, ))
Matrix([
[x**2/2, x*y],
[     x,   0]])
>>> M.integrate((x, 0, 2))
Matrix([
[2, 2*y],
[2,   0]])

See also

limit, diff

jacobian(X)

Calculates the Jacobian matrix (derivative of a vector-valued function).

Parameters

``self`` : vector of expressions representing functions f_i(x_1, …, x_n).

X : set of x_i’s in order, it can be a list or a Matrix

Both ``self`` and X can be a row or a column matrix in any order

(i.e., jacobian() should always work).

Examples

>>> from sympy import sin, cos, Matrix
>>> from sympy.abc import rho, phi
>>> X = Matrix([rho*cos(phi), rho*sin(phi), rho**2])
>>> Y = Matrix([rho, phi])
>>> X.jacobian(Y)
Matrix([
[cos(phi), -rho*sin(phi)],
[sin(phi),  rho*cos(phi)],
[   2*rho,             0]])
>>> X = Matrix([rho*cos(phi), rho*sin(phi)])
>>> X.jacobian(Y)
Matrix([
[cos(phi), -rho*sin(phi)],
[sin(phi),  rho*cos(phi)]])

See also

hessian, wronskian

limit(*args)

Calculate the limit of each element in the matrix. args will be passed to the limit function.

Examples

>>> from sympy.matrices import Matrix
>>> from sympy.abc import x, y
>>> M = Matrix([[x, y], [1, 0]])
>>> M.limit(x, 2)
Matrix([
[2, y],
[1, 0]])

See also

integrate, diff

MatrixBase Class Reference

class sympy.matrices.matrices.MatrixBase

Base class for matrix objects.

D

Return Dirac conjugate (if self.rows == 4).

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.

See also

conjugate

By-element conjugation

H

Hermite conjugation

LDLdecomposition(hermitian=True)

Returns the LDL Decomposition (L, D) of matrix A, such that L * D * L.H == A if hermitian flag is True, or L * D * L.T == A if hermitian is False. This method eliminates the use of square root. Further this ensures that all the diagonal entries of L are 1. A must be a Hermitian positive-definite matrix if hermitian is True, or a symmetric matrix otherwise.

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

The matrix can have complex entries:

>>> from sympy import I
>>> A = Matrix(((9, 3*I), (-3*I, 5)))
>>> L, D = A.LDLdecomposition()
>>> L
Matrix([
[   1, 0],
[-I/3, 1]])
>>> D
Matrix([
[9, 0],
[0, 4]])
>>> L*D*L.H == A
True

See also

cholesky, LUdecomposition, QRdecomposition

LDLsolve(rhs)

Solves Ax = B using LDL decomposition, for a general square and non-singular matrix.

For a non-square matrix with rows > cols, the least squares solution is returned.

Examples

>>> from sympy.matrices import Matrix, eye
>>> A = eye(2)*2
>>> B = Matrix([[1, 2], [3, 4]])
>>> A.LDLsolve(B) == B/2
True

See also

LDLdecomposition, lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LUsolve, QRsolve, pinv_solve

LUdecomposition(iszerofunc=<function _iszero>, simpfunc=None, rankcheck=False)

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.

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]])

See also

cholesky, LDLdecomposition, QRdecomposition, LUdecomposition_Simple, LUdecompositionFF, LUsolve

LUdecompositionFF()

Compute a fraction-free 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, “Fraction-free matrix factors: new forms for LU and QR factors”. Frontiers in Computer Science in China, Vol 2, no. 1, pp. 67-80, 2008.

See also

LUdecomposition, LUdecomposition_Simple, LUsolve

LUdecomposition_Simple(iszerofunc=<function _iszero>, simpfunc=None, rankcheck=False)

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 P=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

LUdecomposition, LUdecompositionFF, LUsolve

LUsolve(rhs, iszerofunc=<function _iszero>)

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.

See also

lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LDLsolve, QRsolve, pinv_solve, LUdecomposition

QRdecomposition()

Return Q, R where A = Q*R, Q is orthogonal and R is upper triangular.

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]])

See also

cholesky, LDLdecomposition, LUdecomposition, QRsolve

QRsolve(b)

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 floating-point 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.

See also

lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LDLsolve, LUsolve, pinv_solve, QRdecomposition

add(b)

Return self + b

cholesky(hermitian=True)

Returns the Cholesky-type decomposition L of a matrix A such that L * L.H == A if hermitian flag is True, or L * L.T == A if hermitian is False.

A must be a Hermitian positive-definite matrix if hermitian is True, or a symmetric matrix if it is False.

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]])

The matrix can have complex entries:

>>> from sympy import I
>>> A = Matrix(((9, 3*I), (-3*I, 5)))
>>> A.cholesky()
Matrix([
[ 3, 0],
[-I, 2]])
>>> A.cholesky() * A.cholesky().H
Matrix([
[   9, 3*I],
[-3*I,   5]])

Non-hermitian Cholesky-type decomposition may be useful when the matrix is not positive-definite.

>>> A = Matrix([[1, 2], [2, 1]])
>>> L = A.cholesky(hermitian=False)
>>> L
Matrix([
[1,         0],
[2, sqrt(3)*I]])
>>> L*L.T == A
True

See also

LDLdecomposition, LUdecomposition, QRdecomposition

cholesky_solve(rhs)

Solves Ax = B using Cholesky decomposition, for a general square non-singular matrix. For a non-square matrix with rows > cols, the least squares solution is returned.

See also

lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, diagonal_solve, LDLsolve, LUsolve, QRsolve, pinv_solve

condition_number()

Returns the condition number of a matrix.

This is the maximum singular value divided by the minimum singular value

Examples

>>> from sympy import Matrix, S
>>> A = Matrix([[1, 0, 0], [0, 10, 0], [0, 0, S.One/10]])
>>> A.condition_number()
100

See also

singular_values

copy()

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)

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 \times b = - b \times a\)) will hold.

Parameters

b : 3x1 or 1x3 Matrix

See also

dot, multiply, multiply_elementwise

diagonal_solve(rhs)

Solves Ax = B efficiently, where A is a diagonal Matrix, with non-zero diagonal entries.

Examples

>>> from sympy.matrices import Matrix, eye
>>> A = eye(2)*2
>>> B = Matrix([[1, 2], [3, 4]])
>>> A.diagonal_solve(B) == B/2
True

See also

lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, cholesky_solve, LDLsolve, LUsolve, QRsolve, pinv_solve

dot(b, hermitian=None, conjugate_convention=None)

Return the dot or inner product of two vectors of equal length. Here self must be a Matrix of size 1 x n or n x 1, and b must be either a matrix of size 1 x n, n x 1, or a list/tuple of length n. A scalar is returned.

By default, dot does not conjugate self or b, even if there are complex entries. Set hermitian=True (and optionally a conjugate_convention) to compute the hermitian inner product.

Possible kwargs are hermitian and conjugate_convention.

If conjugate_convention is "left", "math" or "maths", the conjugate of the first vector (self) is used. If "right" or "physics" is specified, the conjugate of the second vector b is used.

Examples

>>> from sympy import Matrix
>>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> v = Matrix([1, 1, 1])
>>> M.row(0).dot(v)
6
>>> M.col(0).dot(v)
12
>>> v = [3, 2, 1]
>>> M.row(0).dot(v)
10

>>> from sympy import I
>>> q = Matrix([1*I, 1*I, 1*I])
>>> q.dot(q, hermitian=False)
-3

>>> q.dot(q, hermitian=True)
3

>>> q1 = Matrix([1, 1, 1*I])
>>> q.dot(q1, hermitian=True, conjugate_convention="maths")
1 - 2*I
>>> q.dot(q1, hermitian=True, conjugate_convention="physics")
1 + 2*I

See also

cross, multiply, multiply_elementwise

dual()

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.

exp()

Return the exponentiation of a square matrix.

gauss_jordan_solve(B, freevar=False)

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 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 arbitrary parameters. These arbitrary parameters are returned as params Matrix.

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]])
>>> taus_zeroes = { tau:0 for tau in params }
>>> sol_unique = sol.xreplace(taus_zeroes)
>>> sol_unique
 Matrix([
[2],
[0],
[5],
[0]])

>>> 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, [])

>>> A = Matrix([[2, -7], [-1, 4]])
>>> B = Matrix([[-21, 3], [12, -2]])
>>> sol, params = A.gauss_jordan_solve(B)
>>> sol
Matrix([
[0, -2],
[3, -1]])
>>> params
Matrix(0, 2, [])

See also

lower_triangular_solve, upper_triangular_solve, cholesky_solve, diagonal_solve, LDLsolve, LUsolve, QRsolve, pinv

References

R458

https://en.wikipedia.org/wiki/Gaussian_elimination

inv(method=None, **kwargs)

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.

Notes

According to the method keyword, it calls the appropriate method:

LDL … inverse_LDL(); default CH …. inverse_CH()

Kwargs

method : (‘CH’, ‘LDL’)

See also

inverse_LU, inverse_GE, inverse_ADJ

inv_mod(m)

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>)

Calculates the inverse using the adjugate matrix and a determinant.

See also

inv, inverse_LU, inverse_GE

inverse_GE(iszerofunc=<function _iszero>)

Calculates the inverse using Gaussian elimination.

See also

inv, inverse_LU, inverse_ADJ

inverse_LU(iszerofunc=<function _iszero>)

Calculates the inverse using LU decomposition.

See also

inv, inverse_GE, inverse_ADJ

is_nilpotent()

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)

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

key2ij(key)

Converts key into canonical form, converting integers or indexable items into valid integers for self’s range or returning slices unchanged.

See also

key2bounds

lower_triangular_solve(rhs)

Solves Ax = B, where A is a lower triangular matrix.

See also

upper_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LDLsolve, LUsolve, QRsolve, pinv_solve

multiply(b)

Returns self*b

See also

dot, cross, multiply_elementwise

norm(ord=None)

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	2-norm
‘fro’	Frobenius norm	does not exist
inf	maximum row sum	max(abs(x))
-inf	–	min(abs(x))
1	maximum column sum	as below
-1	–	as below
2	2-norm (largest sing. value)	as below
-2	smallest singular value	as below
other	does not exist	sum(abs(x)**ord)**(1./ord)

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(1) # maximum sum of absolute values of A is 2
2
>>> A.norm(2) # Spectral norm (max of |Ax|/|x| under 2-vector-norm)
2
>>> A.norm(-2) # Inverse spectral norm (smallest singular value)
0
>>> A.norm() # Frobenius Norm
2
>>> A.norm(oo) # Infinity Norm
2
>>> Matrix([1, -2]).norm(oo)
2
>>> Matrix([-1, 2]).norm(-oo)
1

See also

normalized

normalized(iszerofunc=<function _iszero>)

Return the normalized version of self.

Parameters

iszerofunc : Function, optional

A function to determine whether self is a zero vector. The default _iszero tests to see if each element is exactly zero.

Returns

Matrix

Normalized vector form of self. It has the same length as a unit vector. However, a zero vector will be returned for a vector with norm 0.

Raises

ShapeError

If the matrix is not in a vector form.

See also

norm

pinv()

Calculate the Moore-Penrose pseudoinverse of the matrix.

The Moore-Penrose pseudoinverse exists and is unique for any matrix. If the matrix is invertible, the pseudoinverse is the same as the inverse.

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]])

See also

inv, pinv_solve

References

R459

https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse

pinv_solve(B, arbitrary_matrix=None)

Solve Ax = B using the Moore-Penrose 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 least-squares 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.

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]])

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 least-squares 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.

See also

lower_triangular_solve, upper_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LDLsolve, LUsolve, QRsolve, pinv

References

R460

https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system

print_nonzero(symb='X')

Shows location of non-zero 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)

Return the projection of self onto the line containing v.

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]])

rank_decomposition(iszerofunc=<function _iszero>, simplify=False)

Returns a pair of matrices (\(C\), \(F\)) with matching rank such that \(A = C F\).

Parameters

iszerofunc : Function, optional

A function used for detecting whether an element can act as a pivot. lambda x: x.is_zero is used by default.

simplify : Bool or Function, optional

A function used to simplify elements when looking for a pivot. By default SymPy’s simplify is used.

Returns

(C, F) : Matrices

\(C\) and \(F\) are full-rank matrices with rank as same as \(A\), whose product gives \(A\).

See Notes for additional mathematical details.

Examples

>>> from sympy.matrices import Matrix
>>> A = Matrix([
...     [1, 3, 1, 4],
...     [2, 7, 3, 9],
...     [1, 5, 3, 1],
...     [1, 2, 0, 8]
... ])
>>> C, F = A.rank_decomposition()
>>> C
Matrix([
[1, 3, 4],
[2, 7, 9],
[1, 5, 1],
[1, 2, 8]])
>>> F
Matrix([
[1, 0, -2, 0],
[0, 1,  1, 0],
[0, 0,  0, 1]])
>>> C * F == A
True

Notes

Obtaining \(F\), an RREF of \(A\), is equivalent to creating a product

\[E_n E_{n-1} ... E_1 A = F\]

where \(E_n, E_{n-1}, ... , E_1\) are the elimination matrices or permutation matrices equivalent to each row-reduction step.

The inverse of the same product of elimination matrices gives \(C\):

\[C = (E_n E_{n-1} ... E_1)^{-1}\]

It is not necessary, however, to actually compute the inverse: the columns of \(C\) are those from the original matrix with the same column indices as the indices of the pivot columns of \(F\).

See also

rref

References

R461

https://en.wikipedia.org/wiki/Rank_factorization

R462

Piziak, R.; Odell, P. L. (1 June 1999). “Full Rank Factorization of Matrices”. Mathematics Magazine. 72 (3): 193. doi:10.2307/2690882

solve(rhs, method='GJ')

Solves linear equation where the unique solution exists.

Parameters

rhs : Matrix

Vector representing the right hand side of the linear equation.

method : string, optional

If set to 'GJ', the Gauss-Jordan elimination will be used, which is implemented in the routine gauss_jordan_solve.

If set to 'LU', LUsolve routine will be used.

If set to 'QR', QRsolve routine will be used.

If set to 'PINV', pinv_solve routine will be used.

It also supports the methods available for special linear systems

For positive definite systems:

If set to 'CH', cholesky_solve routine will be used.

If set to 'LDL', LDLsolve routine will be used.

To use a different method and to compute the solution via the inverse, use a method defined in the .inv() docstring.

Returns

solutions : Matrix

Vector representing the solution.

Raises

ValueError

If there is not a unique solution then a ValueError will be raised.

If self is not square, a ValueError and a different routine for solving the system will be suggested.

solve_least_squares(rhs, method='CH')

Return the least-square fit to the data.

Parameters

rhs : Matrix

Vector representing the right hand side of the linear equation.

method : string or boolean, optional

If set to 'CH', cholesky_solve routine will be used.

If set to 'LDL', LDLsolve routine will be used.

If set to 'QR', QRsolve routine will be used.

If set to 'PINV', pinv_solve routine will be used.

Otherwise, the conjugate of self will be used to create a system of equations that is passed to solve along with the hint defined by method.

Returns

solutions : Matrix

Vector representing the solution.

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 least-squares 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')

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}

upper_triangular_solve(rhs)

Solves Ax = B, where A is an upper triangular matrix.

See also

lower_triangular_solve, gauss_jordan_solve, cholesky_solve, diagonal_solve, LDLsolve, LUsolve, QRsolve, pinv_solve

vech(diagonal=True, check_symmetry=True)

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

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]])

See also

vec

Matrix Exceptions Reference

class sympy.matrices.matrices.MatrixError

class sympy.matrices.matrices.ShapeError

Wrong matrix shape

class sympy.matrices.matrices.NonSquareMatrixError

Matrix Functions Reference

sympy.matrices.dense.matrix_multiply_elementwise(A, B)

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)

Returns a matrix of zeros with rows rows and cols columns; if cols is omitted a square matrix will be returned.

See also

ones, eye, diag

sympy.matrices.dense.ones(*args, **kwargs)

Returns a matrix of ones with rows rows and cols columns; if cols is omitted a square matrix will be returned.

See also

zeros, eye, diag

sympy.matrices.dense.eye(*args, **kwargs)

Create square identity matrix n x n

See also

diag, zeros, ones

sympy.matrices.dense.diag(*values, **kwargs)

Create a sparse, diagonal matrix from a list of diagonal values.

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'>

See also

eye

sympy.matrices.dense.jordan_cell(eigenval, n)

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=[])

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.

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           ]

See also

sympy.matrices.mutable.Matrix.jacobian, wronskian

References

https://en.wikipedia.org/wiki/Hessian_matrix

sympy.matrices.dense.GramSchmidt(vlist, orthonormal=False)

Apply the Gram-Schmidt process to a set of vectors.

see: https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

sympy.matrices.dense.wronskian(functions, var, method='bareiss')

Compute Wronskian for [] of functions

                 | f1       f2        ...   fn      |
                 | f1'      f2'       ...   fn'     |
                 |  .        .        .      .      |
W(f1, ..., fn) = |  .        .         .     .      |
                 |  .        .          .    .      |
                 |  (n)      (n)            (n)     |
                 | D   (f1) D   (f2)  ...  D   (fn) |

see: https://en.wikipedia.org/wiki/Wronskian

See also

sympy.matrices.mutable.Matrix.jacobian, hessian

sympy.matrices.dense.casoratian(seqs, n, zero=True)

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+k-1) b(n+k-1) . . . z(n+k-1) +

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)

Create random matrix with dimensions r x c. If c is omitted the matrix will be square. If symmetric is True the matrix must be square. If percent is less than 100 then only approximately the given percentage of elements will be non-zero.

The pseudo-random 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 of random.Random, or at least have randint and shuffle methods with same signatures.

• if prng is not supplied but seed is supplied, then new random.Random with given seed will be created;

• otherwise, a new random.Random with default seed will be used.

Examples

>>> from sympy.matrices import randMatrix
>>> randMatrix(3) # doctest:+SKIP
[25, 45, 27]
[44, 54,  9]
[23, 96, 46]
>>> randMatrix(3, 2) # doctest:+SKIP
[87, 29]
[23, 37]
[90, 26]
>>> randMatrix(3, 3, 0, 2) # doctest:+SKIP
[0, 2, 0]
[2, 0, 1]
[0, 0, 1]
>>> randMatrix(3, symmetric=True) # doctest:+SKIP
[85, 26, 29]
[26, 71, 43]
[29, 43, 57]
>>> A = randMatrix(3, seed=1)
>>> B = randMatrix(3, seed=2)
>>> A == B # doctest:+SKIP
False
>>> A == randMatrix(3, seed=1)
True
>>> randMatrix(3, symmetric=True, percent=50) # doctest:+SKIP
[77, 70,  0],
[70,  0,  0],
[ 0,  0, 88]

Numpy Utility Functions Reference

sympy.matrices.dense.list2numpy(l, dtype=<class 'object'>)

Converts python list of SymPy expressions to a NumPy array.

See also

matrix2numpy

sympy.matrices.dense.matrix2numpy(m, dtype=<class 'object'>)

Converts SymPy’s matrix to a NumPy array.

See also

list2numpy

sympy.matrices.dense.symarray(prefix, shape, **kwargs)

Create a numpy ndarray of symbols (as an object array).

The created symbols are named prefix_i1_i2_… You should thus provide a non-empty 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 one-dimensional; 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]]
<BLANKLINE>
 [[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)

Returns a rotation matrix for a rotation of theta (in radians) about the 1-axis.

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]])

See also

rot_axis2

Returns a rotation matrix for a rotation of theta (in radians) about the 2-axis

rot_axis3

Returns a rotation matrix for a rotation of theta (in radians) about the 3-axis

sympy.matrices.dense.rot_axis2(theta)

Returns a rotation matrix for a rotation of theta (in radians) about the 2-axis.

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]])

See also

rot_axis1

Returns a rotation matrix for a rotation of theta (in radians) about the 1-axis

rot_axis3

Returns a rotation matrix for a rotation of theta (in radians) about the 3-axis

sympy.matrices.dense.rot_axis3(theta)

Returns a rotation matrix for a rotation of theta (in radians) about the 3-axis.

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]])

See also

rot_axis1

Returns a rotation matrix for a rotation of theta (in radians) about the 1-axis

rot_axis2

Returns a rotation matrix for a rotation of theta (in radians) about the 2-axis

sympy.matrices.matrices.a2idx(j, n=None)