The linear algebra module is designed to be as simple as possible. First, we import and declare our first Matrix object:
>>> from sympy.matrices import Matrix
>>> Matrix([[1,0], [0,1]])
⎡1 0⎤
⎣0 1⎦
This is the standard manner one creates a matrix, i.e. with a list of appropriately-sizes lists. 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 interestingly (and usefully), we can use a 2-variable function (or lambda) to make one. 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:
>>> 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]
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:
>>> M[1,2]
6
>>> M[0,0]
1
>>> M[1,1]
5
Since this is Python we’re also able to slice submatrices:
>>> M[0:2,0:2]
⎡1 2⎤
⎣4 5⎦
>>> M[1:2,2:2]
>>> M[:,2]
⎡3⎤
⎣6⎦
Remember in the 2nd example above that slicing 2:2 gives an empty range and that, as in python, a 4 column list is indexed from 0 to 3. In particular, this mean a quick way to create a copy of the matrix is:
>>> M2 = M[:,:]
>>> M2[0,0] = 100
>>> M
⎡1 2 3⎤
⎣4 5 6⎦
See? 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⎦
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!:
>>> 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⎦
Now that we have the basics out of the way, let’s see what we can do with the actual matrices. Of course the first things that come to mind are the basics like 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
and the inverse. In SymPy the inverse is computed by Gaussian elimination by default 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("LU")
⎡1 0 0⎤
⎢0 1 0⎥
⎣0 0 1⎦
>>> M.inv("LU")
⎡-3/14 1/14 1/2⎤
⎢-1/28 5/28 -1/4⎥
⎣ 3/7 -1/7 0⎦
>>> M * M.inv("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
⎡ ⎤
⎢ ⎽⎽⎽ ⎥
⎢ ⎽⎽⎽⎽⎽ -╲╱ 3 ⎽⎽⎽⎽⎽⎥
⎢ ╲╱ 1/6 ────── -╲╱ 1/2 ⎥
⎢ 3 ⎥
⎢ ⎥
⎢ ⎽⎽⎽ ⎥
⎢ ⎽⎽⎽⎽⎽ -╲╱ 3 ⎽⎽⎽⎽⎽⎥
⎢ ╲╱ 1/6 ────── ╲╱ 1/2 ⎥
⎢ 3 ⎥
⎢ ⎥
⎢ ⎽⎽⎽ ⎥
⎢ ⎽⎽⎽⎽⎽ ╲╱ 3 ⎥
⎢2*╲╱ 1/6 ───── 0⎥
⎣ 3 ⎦
>>> R
⎡ ⎽⎽⎽ ⎽⎽⎽⎽⎽ ⎽⎽⎽⎽⎽⎤
⎢ ╲╱ 6 8*╲╱ 1/6 12*╲╱ 1/6 ⎥
⎢ ⎽⎽⎽⎽⎽ ⎥
⎢ 0 ╲╱ 1/3 0⎥
⎢ ⎽⎽⎽⎥
⎣ 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 then with respect to another 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
...
[2]
[3]
[5]
[ 23/19]
[ 63/19]
[-47/19]
[ 1692/353]
[-1551/706]
[ -423/706]
>>> for i in out2:
... print i
...
[2*(1/38)**(1/2)]
[3*(1/38)**(1/2)]
[5*(1/38)**(1/2)]
[(23/19)*19**(1/2)*(1/353)**(1/2)]
[(63/19)*19**(1/2)*(1/353)**(1/2)]
[ -47/19*19**(1/2)*(1/353)**(1/2)]
[(12/353)*706**(1/2)]
[ -11/706*706**(1/2)]
[ -3/706*706**(1/2)]
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.
Returns 4 matrices P, L, D, U such that PA = L D**-1 U.
From the paper “fraction-free matrix factors...” by Zhou and Jeffrey
Return Q*R where Q is orthogonal and R is upper triangular.
Assumes full-rank square (for now).
Returns the adjugate matrix.
Adjugate matrix is the transpose of the cofactor matrix.
http://en.wikipedia.org/wiki/Adjugate
See also: .cofactorMatrix(), .T
>>> from sympy import *
>>> m = Matrix(2,2,lambda i,j: i*2+j)
>>> m #doctest: +NORMALIZE_WHITESPACE
[0, 1]
[2, 3]
>>> m.applyfunc(lambda i: 2*i) #doctest: +NORMALIZE_WHITESPACE
[0, 2]
[4, 6]
The Berkowitz algorithm.
Given N x N matrix with symbolic content, compute efficiently coefficients of characteristic polynomials of ‘self’ and all its square sub-matrices composed by removing both i-th row and column, without division in the ground domain.
This method is particulary useful for computing determinant, principal minors and characteristic polynomial, when ‘self’ has complicated coefficients eg. polynomials. Semi-direct usage of this algorithm is also important in computing efficiently subresultant PRS.
Assuming that M is a square matrix of dimension N x N and I is N x N identity matrix, then the following following definition of characteristic polynomial is begin used:
charpoly(M) = det(t*I - M)
As a consequence, all polynomials generated by Berkowitz algorithm are monic.
>>> from sympy import *
>>> x,y,z = symbols('xyz')
>>> M = Matrix([ [x,y,z], [1,0,0], [y,z,x] ])
>>> p, q, r = M.berkowitz()
>>> print p # 1 x 1 M's sub-matrix
(1, -x)
>>> print q # 2 x 2 M's sub-matrix
(1, -x, -y)
>>> print r # 3 x 3 M's sub-matrix
(1, -2*x, -y - y*z + x**2, x*y - z**2)
For more information on the implemented algorithm refer to:
>>> import sympy
>>> M = sympy.matrices.eye(3)
>>> M.col_del(1)
>>> M #doctest: +NORMALIZE_WHITESPACE
[1, 0]
[0, 0]
[0, 1]
>>> from sympy import *
>>> M = Matrix(3,3,lambda i,j: i+j)
>>> M
[0, 1, 2]
[1, 2, 3]
[2, 3, 4]
>>> V = zeros((3, 1))
>>> V
[0]
[0]
[0]
>>> M.col_insert(1,V)
[0, 0, 1, 2]
[1, 0, 2, 3]
[2, 0, 3, 4]
Concatenates two matrices along self’s last and bott’s first row
>>> from sympy import *
>>> M = Matrix(3,3,lambda i,j: i+j)
>>> V = Matrix(1,3,lambda i,j: 3+i+j)
>>> M.col_join(V)
[0, 1, 2]
[1, 2, 3]
[2, 3, 4]
[3, 4, 5]
Computes the matrix determinant using the method “method”.
Compute matrix determinant using Bareis’ fraction-free algorithm which is an extension of the well known Gaussian elimination method. This approach is best suited for dense symbolic matrices and will result in a determinant with minimal numer of fractions. It means that less term rewriting is needed on resulting formulae.
TODO: Implement algorithm for sparse matrices (SFF).
Calculates the matrix inverse.
According to the “method” parameter, it calls the appropriate method:
GE .... inverse_GE() LU .... inverse_LU() ADJ ... inverse_ADJ()
Calculates the Jacobian matrix (derivative of a vectorial function).
self is a vector of expression representing functions f_i(x_1, ..., x_n). varlist is the set of x_i’s in order.
>>> from sympy import *
>>> m = Matrix(2,3,lambda i,j: 1)
>>> m #doctest: +NORMALIZE_WHITESPACE
[1, 1, 1]
[1, 1, 1]
>>> m.reshape(1,6) #doctest: +NORMALIZE_WHITESPACE
[1, 1, 1, 1, 1, 1]
>>> m.reshape(3,2) #doctest: +NORMALIZE_WHITESPACE
[1, 1]
[1, 1]
[1, 1]
>>> from sympy import *
>>> M = Matrix(3,3,lambda i,j: i+j)
>>> M
[0, 1, 2]
[1, 2, 3]
[2, 3, 4]
>>> V = zeros((1, 3))
>>> V
[0, 0, 0]
>>> M.row_insert(1,V)
[0, 1, 2]
[0, 0, 0]
[1, 2, 3]
[2, 3, 4]
Concatenates two matrices along self’s last and rhs’s first column
>>> from sympy import *
>>> M = Matrix(3,3,lambda i,j: i+j)
>>> V = Matrix(3,1,lambda i,j: 3+i+j)
>>> M.row_join(V)
[0, 1, 2, 3]
[1, 2, 3, 4]
[2, 3, 4, 5]
Take any matrix and return reduced row-echelon form and indices of pivot vars
To simplify elements before finding nonzero pivots set simplified=True
>>> from sympy import *
>>> m = Matrix(4,4,lambda i,j: i+j)
>>> m #doctest: +NORMALIZE_WHITESPACE
[0, 1, 2, 3]
[1, 2, 3, 4]
[2, 3, 4, 5]
[3, 4, 5, 6]
>>> m[0:1, 1] #doctest: +NORMALIZE_WHITESPACE
[1]
>>> m[0:2, 0:1] #doctest: +NORMALIZE_WHITESPACE
[0]
[1]
>>> m[2:4, 2:4] #doctest: +NORMALIZE_WHITESPACE
[4, 5]
[5, 6]
Return the Matrix converted in a python list.
>>> from sympy import *
>>> m=Matrix(3, 3, range(9))
>>> m
[0, 1, 2]
[3, 4, 5]
[6, 7, 8]
>>> m.tolist()
[[0, 1, 2], [3, 4, 5], [6, 7, 8]]
Matrix transposition.
>>> from sympy import *
>>> m=Matrix(((1,2+I),(3,4)))
>>> m #doctest: +NORMALIZE_WHITESPACE
[1, 2 + I]
[3, 4]
>>> m.transpose() #doctest: +NORMALIZE_WHITESPACE
[ 1, 3]
[2 + I, 4]
>>> m.T == m.transpose()
True