# css styling of the notebook - execute and then ignore
import requests, IPython.core.display
IPython.core.display.HTML(requests.get("https://git.io/fj7u5").text)
SymPy is a Python module for performing symbolic computations. This notebook explain how use it to perform linear algebra computations.
Execute this cell to load SymPy:
#load SymPy module content
from sympy import *
#this makes printouts of matrices and vectors more readeable:
init_printing(use_latex='mathjax')
The simplest way to create a matrix is to lists its elements. Matrices are entered row by row. Each row must be enclosed in square brackets, and then the whole collection of rows is again enclosed in square brackets:
A = Matrix([[1, 2, 3], [4, 5, 6], [7,8,9]])
A #this displays the matrix
Note. In order to create a matrix with one row we need to enclose it entries in double square brackets (one set indicating the beginning and the end of the row, and another enclosing the whole matrix:
A = Matrix([[12, 13, 14]])
A
Using only one set of square brackets creates a column vector, i.e. a matrix with one column:
v = Matrix([12, 13, 14])
v
Say that we want to enter the matrix
\begin{equation*} \begin{bmatrix} \frac{1}{7} & -1 \\ 0 & -\frac{1}{3} \\ \end{bmatrix} \end{equation*}We can try to do it as follows:
B = Matrix([[1/7 , -1],[0, -1/3]])
B
As you see fractions got converted to decimals, and rounded to some number of decimal places. For some applications in this course this would be inconvenient, since we will want to work with fractions and keep their exact values. We can deal with this issue as follows:
B = Matrix([[Rational(1,7) , -1],[0, Rational(-1,3)]])
B
For larger matrices it is tedious to enter all their entries manually. It is often faster to create first a base matrix whose entries are e.g. all zeros or all ones, and then modify entries of this matrix as needed. It is easy to create such base matrices of any size:
C = zeros(3, 4) # 3 rows, 4 columns
C
D = ones(2, 7) # 2 rows, 7 columns
D
A diagonal matrix is a square matrix whose all entries outside its main diagonal are zeros. Such matrix can be created as follows:
D = diag(1, 2, 3, -3)
D
The identity matrix is a diagonal matrix whose all entries on the main diagonal are ones. Since this matrix is used very often, there is a special way to create it:
I = eye(4) # 4x4 identity matrix
I
In some cases it may be convenient to enter all entries of a matrix as a list, and then to reshape this list into a matrix:
#create a list of elements
a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
a
#reshape the list into a 3x4 matrix
A = Matrix(a).reshape(3, 4)
A
#reshape the list into a 2x6 matrix
B = Matrix(a).reshape(2, 6)
B
Note 1. Reshaping works only if the number of elements of the list is exactly equal to the number of entries of a matrix of specified dimensions.
Note 2. Reshaping makes consecutive chunks of a list into rows of a matrix. In order for such chunks to become columns of a matrix, we can apply transpose of a matrix, which makes rows into columns:
A
# B is the transpose of A, rows of A are columns of B
B = A.T
B
Here is a sample matrix we will be working with:
A = Matrix([[1, 2, 3], [4, 5, 6], [7,8,9]])
A
We can access elements of a matrix by specifying the row and column of the element.
Note. In Python indexing starts with 0. Thus, the first row of a matrix has index 0, the second row has index 1 and so on.
A[1,0] #the element in row 1, column 0
We can use this to modify elements of a matrix
A[1,1] = 100
A
Accessing matrix rows:
A[0, :] # get row 0
Accessing columns:
A[:, 2] # get column 2
We can modify the whole row (or column) of a matrix at once:
A
w = Matrix([[-10, -20, -30]])
w
A[2, :] = w # make row 2 of A equal to w
A
Matrices can be modified by inserting into them additional rows and columns. For example, we start with the following matrix:
A = Matrix([[1, 2, 3], [4, 5, 6], [7,8,9]])
A
Assume the we want to insert the following column into this matrix:
v = Matrix([-1, -2, -3]) # matrix with one column
v
This can be done as follows:
B = A.col_insert(2, v) # insert v as the column number 2
B
Note. Column insertion creates a new matrix, the matrix A
remains unchanged:
A
Row insertion works the same way:
C = Matrix([[100, 200, 300]]) # matrix with one row
C
A.row_insert(1, C) #insert C as the row number 1
Row/column insertion works also with matrices consisting of more then one row or column:
A
I = eye(3)
I
Insert the matrix I
into A
, so that the first column of I
becomes column number 3 of the new matrix:
B = A.col_insert(3, eye(3))
B
Matrices can be also modified by deleting their rows and columns. Unlike row/column insertions, deletion modifies the original matrix:
A = Matrix([[1, 2, 3], [4, 5, 6], [7,8,9], [10, 11, 12]])
A
A.row_del(0) #delete row 0
A
A.col_del(1) #delete column 1
A
We will work with the following matrix:
A = Matrix([[0 , 1, 2, 3, 4], [5, 6, 7, 8, 9], [9, 10, 11, 12, 13]])
A
We will denote by R$_i$ the row number $i$ of the matrix (remember that in Python indexing starts with 0).
Multiply R$_1$ by $\frac{1}{5}$:
A[1, :] = A[1, :]*Rational(1, 5)
A
Subtract 9R$_1$ from R$_2$:
A[2, :] = A[2, :] - 9*A[1, :]
A
Interchange R$_0$ and R$_1$:
A.row_swap(0,1)
A
Elementary operations on columns can be done in the same way. For example, here we interchange column 0 and column 1:
A.col_swap(0,1)
A
A = Matrix([[1, 2, 3, 12], [4, 5, 6, 13], [7,8,9, 14]])
A
The function rref()
computes the row reduced echelon form of a matrix:
R = A.rref()
R
Note. A.rref()
returns a tuple consisting of a reduced matrix and a list of indices of pivot columns. To get the reduced matrix only we need to select the first entry of this tuple:
B = A.rref()[0]
B
Here is an example showing that row reduction of matrices with decimal entries can give wrong results:
A = Matrix([[1, 3],[0.1, 0.3]])
A
A.rref() # the result is wrong!
Such issues occur because computer operations on numbers in the decimal form involve rounding errors, and this distorts the end result.
We can avoid such problems by entering all entries of a matrix as rational numbers:
A = Matrix([[1, 3],[Rational(1,10), Rational(3,10)]])
A
Now the reduced echelon form is computed correctly:
A.rref()
We will be working with the following matrices A
and B
:
A = Matrix([[0, 1], [2, 3], [4,5]])
A
B = Matrix([[3, 2], [1, 0], [-1,2]])
B
Addition of matrices:
A+B
Multiplication of a matrix by a scalar:
Rational(1, 2)*A
Transpose of a matrix:
A.T
Matrix multiplication:
C = (A.T)*B
C
We can multiply a matrix by a column vector in the same way:
A
v = Matrix([-1, 2])
v
A*v
Matrix inverse:
C.inv()
Not every matrix is invertible. An attempt to compute an inverse of a non-invertible matrix will give an error:
# this matrix is non-invertible
A = Matrix([[1, 0], [1, 0]])
A
A.inv()
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-95-60cc4e20accc> in <module> ----> 1 A.inv() /anaconda3/lib/python3.7/site-packages/sympy/matrices/matrices.py in inv(self, method, **kwargs) 2916 if method is not None: 2917 kwargs['method'] = method -> 2918 return self._eval_inverse(**kwargs) 2919 2920 def is_nilpotent(self): /anaconda3/lib/python3.7/site-packages/sympy/matrices/dense.py in _eval_inverse(self, **kwargs) 264 M = self.as_mutable() 265 if method == "GE": --> 266 rv = M.inverse_GE(iszerofunc=iszerofunc) 267 elif method == "LU": 268 rv = M.inverse_LU(iszerofunc=iszerofunc) /anaconda3/lib/python3.7/site-packages/sympy/matrices/matrices.py in inverse_GE(self, iszerofunc) 2831 red = big.rref(iszerofunc=iszerofunc, simplify=True)[0] 2832 if any(iszerofunc(red[j, j]) for j in range(red.rows)): -> 2833 raise ValueError("Matrix det == 0; not invertible.") 2834 2835 return self._new(red[:, big.rows:]) ValueError: Matrix det == 0; not invertible.
Matrix powers:
Taking the $n$-th power of a matrix has the same effect as multiplying the matrix by itself $n$ times.
A = Matrix([[1, 2], [3, 4]])
Multiplying the matrix by itself 4 times:
A*A*A*A
Taking the 4-th power of $A$:
A**4
A = Matrix([[1, 1, 2], [3, 4, 5], [6,7,8]])
A
Compute the determinant of the matrix A
:
A.det()
A = Matrix([[1, 2, 3, 4],[5, 6, 7, 8], [4, 4, 4, 4]])
A
A basis of the row space:
A.rowspace()
A basis of the column space:
A.columnspace()
A basis of the null space:
A.nullspace()
We define a few column vectors:
u = Matrix([1, 1, 1])
v = Matrix([1, 2, 3])
w = Matrix([1, -1, 0])
u, v, w
Dot product of vectors v
and w
:
v.dot(w)
Gram-Schmidt orthogonalization of the set of vectors $[$ u
, v
, w
$]$:
GramSchmidt([u, v, w])
Adding the option orthonormal=True
we obtain a set of orthonormal vectors:
GramSchmidt([u, v, w], orthonormal=True)
A = Matrix([[3, -2, 4, -2], [5, 3, -3, -2], [5, -2, 2, -2], [5, -2, -3, 3]])
A
# declare a variable lamda (this is not a misspeling,
#"lambda" already has a different meaning in Python)
lamda = symbols('lamda')
# compute the characteristic polynomial of matrix A with lamda as the variable
p = A.charpoly(lamda)
p
Factorization of the polynomial:
factor(p)
Eigenvalues of the matrix A
with their algebraic multiplicities:
A.eigenvals()
Eigenvalues with their algebraic multiplicities and bases of eigenspaces:
A.eigenvects()
Diagonalization of a matrix (returns a tuple consistsing of matrices $P$ and $D$ such that $A = PDP^{-1}$):
A.diagonalize()
Complex diagonalization:
B = Matrix([[0, 1],[-1, 0]])
B
B.diagonalize()
Not every matrix is diagonalizable. For such matrices we will get an error:
# this matrix is not diagonalizable
C = Matrix([[2, 1], [0, 2]])
C
C.diagonalize()
--------------------------------------------------------------------------- MatrixError Traceback (most recent call last) <ipython-input-97-892acc8b8865> in <module> ----> 1 C.diagonalize() /anaconda3/lib/python3.7/site-packages/sympy/matrices/matrices.py in diagonalize(self, reals_only, sort, normalize) 1085 1086 if not self.is_diagonalizable(reals_only=reals_only, clear_cache=False): -> 1087 raise MatrixError("Matrix is not diagonalizable") 1088 1089 eigenvecs = self._cache_eigenvects MatrixError: Matrix is not diagonalizable
A = Matrix([[0, 1, 2, 3], [4, 6, 7, 0], [1, 1, 0, 0]])
A
sympy
does not include a function for computing singular value decomposition, but such function can be imported from from mpmath
:
from mpmath import svd
U, S, V = svd(A)
U
and V
are orthogonal matrices, S
is the vector of singular values:
U
S
matrix( [['10.3117751931661'], ['3.12657593986521'], ['0.944359707876459']])
V
matrix( [['-0.386544961857121', '-0.595730550459802', '-0.701351582272704', '-0.0615959556733527'], ['0.295892137766083', '0.112289055598665', '-0.176603782603437', '-0.932014009965052'], ['0.536006231242966', '0.395939545766543', '-0.661877715712972', '0.34328863309098']])
S
and V
are not displayed nicely since they are mpmath
matrices, and not sympy
matrices. This can be fixed as follows:
import numpy as np
V = Matrix(np.array(V.tolist(),dtype=np.float64))
V
S = diag(*[Float(x) for x in S])
S