Matrices (linear algebra) ========================= .. module:: sympy.matrices.matrices 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 import pprint >>> import sys >>> sys.displayhook = pprint >>> from sympy.matrices import * >>> Matrix([[1,0], [0,1]]) [1 0] [ ] [0 1] >>> Matrix(( ... Matrix(( ... (1, 0, 0), ... (0, 0, 0) ... )), ... (0, 0, -1) ... )) [1 0 0 ] [ ] [0 0 0 ] [ ] [0 0 -1] This is the standard manner one creates a matrix, i.e. with a list of appropriately-sizes 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 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] 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: >>> 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] [6] >>> 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) None >>> M [4 5 6] [ ] [7 8 9] >>> M.col_del(1) None >>> 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 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 [ ___ ___ ___] [\/ 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 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 ... [ 38**(1/2)/19] [3*38**(1/2)/38] [5*38**(1/2)/38] [ 23*6707**(1/2)/6707] [ 63*6707**(1/2)/6707] [-47*6707**(1/2)/6707] [ 12*706**(1/2)/353] [-11*706**(1/2)/706] [ -3*706**(1/2)/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. Matrix Class Reference ---------------------- .. autoclass:: Matrix :members: