Notebook

Linear Algebra using SymPy¶

Introduction¶

This notebook is a short tutorial of Linear Algebra calculation using SymPy. For further information refer to SymPy official tutorial.

You can also check the SymPy in 10 minutes tutorial.

In [1]:

from sympy import *
init_session()

IPython console for SymPy 1.0 (Python 2.7.13-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

A matrix $A \in \mathbb{R}^{m\times n}$ is a rectangular array of real number with $m$ rows and $n$ columns. To specify a matrix $A$ , we specify the values for its components as a list of lists:

In [2]:

A = Matrix([
        [3, 2, -1, 1],
        [2, -2, 4, -2],
        [-1, S(1)/2, -1, 0]])
display(A)

$\left[\begin{matrix}3 & 2 & -1 & 1\\2 & -2 & 4 & -2\\-1 & \frac{1}{2} & -1 & 0\end{matrix}\right]$

We can access the matrix elements using square brackets, we can also use it for submatrices

In [3]:

A[0, 1] # row 0, column 1

Out[3]:

$2$

In [4]:

A[0:2, 0:3] # top-left 2x3 submatrix

Out[4]:

$\left[\begin{matrix}3 & 2 & -1\\2 & -2 & 4\end{matrix}\right]$

We can also create some common matrices. Let us create an identity matrix

In [5]:

eye(2)

Out[5]:

$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$

In [6]:

zeros(2, 3)

Out[6]:

$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$

We can use algebraic operations like addition $+$ , substraction $-$ , multiplication $*$ , and exponentiation $**$ with Matrix objects.

In [7]:

B = Matrix([
        [2, -3, -8],
        [-2, -1, 2],
        [1, 0, -3]])
C = Matrix([
        [sin(x), exp(x**2), 1],
        [0, cos(x), 1/x],
        [1, 0, 2]])

In [8]:

B + C

Out[8]:

$\left[\begin{matrix}\sin{\left (x \right )} + 2 & e^{x^{2}} - 3 & -7\\-2 & \cos{\left (x \right )} - 1 & 2 + \frac{1}{x}\\2 & 0 & -1\end{matrix}\right]$

In [9]:

B ** 2

Out[9]:

$\left[\begin{matrix}2 & -3 & 2\\0 & 7 & 8\\-1 & -3 & 1\end{matrix}\right]$

In [10]:

C ** 2

Out[10]:

$\left[\begin{matrix}\sin^{2}{\left (x \right )} + 1 & e^{x^{2}} \sin{\left (x \right )} + e^{x^{2}} \cos{\left (x \right )} & \sin{\left (x \right )} + 2 + \frac{e^{x^{2}}}{x}\\\frac{1}{x} & \cos^{2}{\left (x \right )} & \frac{1}{x} \cos{\left (x \right )} + \frac{2}{x}\\\sin{\left (x \right )} + 2 & e^{x^{2}} & 5\end{matrix}\right]$

In [11]:

tan(x) * B ** 5

Out[11]:

$\left[\begin{matrix}52 \tan{\left (x \right )} & 27 \tan{\left (x \right )} & - 28 \tan{\left (x \right )}\\- 2 \tan{\left (x \right )} & - \tan{\left (x \right )} & - 78 \tan{\left (x \right )}\\11 \tan{\left (x \right )} & 30 \tan{\left (x \right )} & 57 \tan{\left (x \right )}\end{matrix}\right]$

And the transpose of the matrix, that flips the matrix through its main diagonal:

In [12]:

A.transpose() # the same as A.T

Out[12]:

$\left[\begin{matrix}3 & 2 & -1\\2 & -2 & \frac{1}{2}\\-1 & 4 & -1\\1 & -2 & 0\end{matrix}\right]$

Row operations¶

In [13]:

M = eye(4)

In [14]:

M[1, :] = M[1, :] + 5*M[0, :]

In [15]:

Out[15]:

$\left[\begin{matrix}1 & 0 & 0 & 0\\5 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$

The notation M[1, :] refers to entire rows of the matrix. The first argument specifies the 0-based row index, for example the first row of M is M[0, :]. The code example above implements the row operation $R_2 \leftarrow R_2 + 5R_1$ . To scale a row by a constant $c$ , use the M[1, :] = c*M[1, :]. To swap rows $1$ and $j$ , we can use the Python tuple-assignment syntax M[1, :], M[j, :] = M[j, :], M[1, :].

Reduced row echelon form¶

The Gauss-Jordan elimination procedure is a sequence of row operations that can be performed on any matrix to bring it to its reduced row echelon form (RREF). In Sympy, matrices have a rref method that compute it:

In [16]:

A.rref()

Out[16]:

$\left ( \left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & -2\\0 & 0 & 1 & -2\end{matrix}\right], \quad \left [ 0, \quad 1, \quad 2\right ]\right )$

It return a tuple, the first value is the RREF of the matrix $A$ , and the second tells the location of the leading ones (pivots). If we just want the RREF, we can just get the first entry of the matrix, i.e.

In [17]:

A.rref()[0]

Out[17]:

$\left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & -2\\0 & 0 & 1 & -2\end{matrix}\right]$

Matrix fundamental spaces¶

Consider the matrix $A \in \mathbb{R}^{m\times n}$ . The fundamental spaces of a matrix are its column space $\mathcal{C}(A)$ , its null space $\mathcal{N}(A)$ , and its row space $\mathcal{R}(A)$ . These vector spaces are importan when we consider the matrix product $A\mathbf{x} = \mathbf{y}$ as a linear transformation $T_A:\mathbb{R}^n\rightarrow \mathbb{R}^n$ of the input vector $\mathbf{x}\in\mathbb{R}^n$ to produce an output vector $\mathbf{y} \in \mathbb{R}^m$ .

Linear transformations $T_A: \mathbb{R}^n \rightarrow \mathbb{R}^m$ can be represented as $m\times n$ matrices. The fundamental spaces of a matrix $A$ gives us information about the domain and image of the linear transformation $T_A$ . The column space $\mathcal{C}(A)$ is the same as the image space $\mathrm{Im}(T_A)$ (the set of all possible outputs). The null space $\mathcal{N}(A)$ is also called kernel $\mathrm{Ker}(T_A)$ , and is the set of all input vectors that are mapped to the zero vector. The row space $\mathcal{R}(A)$ is the orthogonal complement of the null space, i.e., the vectors that are mapped to vectors different from zero. Input vectors in the row space of $A$ are in a one-to-one correspondence with the output vectors in the column space of $A$ .

Let us see how to compute these spaces, or a base for them!

The non-zero rows in the reduced row echelon form $A$ are a basis for its row space, i.e.

In [18]:

[A.rref()[0][row, :] for row in A.rref()[1]]

Out[18]:

$\left [ \left[\begin{matrix}1 & 0 & 0 & 1\end{matrix}\right], \quad \left[\begin{matrix}0 & 1 & 0 & -2\end{matrix}\right], \quad \left[\begin{matrix}0 & 0 & 1 & -2\end{matrix}\right]\right ]$

The column space of $A$ is the span of the columns of $A$ that contain the pivots.

In [19]:

[A[:, col] for col in A.rref()[1]]

Out[19]:

$\left [ \left[\begin{matrix}3\\2\\-1\end{matrix}\right], \quad \left[\begin{matrix}2\\-2\\\frac{1}{2}\end{matrix}\right], \quad \left[\begin{matrix}-1\\4\\-1\end{matrix}\right]\right ]$

We can also use the columnspace method

In [20]:

A.columnspace()

Out[20]:

$\left [ \left[\begin{matrix}3\\2\\-1\end{matrix}\right], \quad \left[\begin{matrix}2\\-2\\\frac{1}{2}\end{matrix}\right], \quad \left[\begin{matrix}-1\\4\\-1\end{matrix}\right]\right ]$

Note that we took columns from the original matrix and not from its RREF.

To find (a base for) the null space of $A$ we use the nullspace method:

In [21]:

A.nullspace()

Out[21]:

$\left [ \left[\begin{matrix}-1\\2\\2\\1\end{matrix}\right]\right ]$

Determinants¶

The determinant of a matrix, denoted by $\det(A)$ or $|A|$ , isis a useful value that can be computed from the elements of a square matrix. It can be viewed as the scaling factor of the transformation described by the matrix.

In [22]:

M = Matrix([
        [1, 2, 2],
        [4, 5, 6],
        [7, 8, 9]])

In [23]:

M.det()

Out[23]:

$3$

Matrix inverse¶

For invertible matrices (those with $\det(A)\neq 0$ ), there is an inverse matrix $A^{-1}$ that have the inverse effect (if we are thinking about linear transformations).

In [24]:

A = Matrix([
        [1, -1, -1],
        [0, 1, 0],
        [1, -2, 1]])

In [25]:

A.inv()

Out[25]:

$\left[\begin{matrix}\frac{1}{2} & \frac{3}{2} & \frac{1}{2}\\0 & 1 & 0\\- \frac{1}{2} & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]$

In [26]:

A.inv() * A

Out[26]:

$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

In [27]:

A * A.inv()

Out[27]:

$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

Eigenvectors and Eigenvalues¶

To find the eigenvalues of a matrix, use eigenvals. eigenvals returns a dictionary of eigenvalue:algebraic multiplicity.

In [28]:

M = Matrix([
        [3, -2,  4, -2],
        [5,  3, -3, -2],
        [5, -2,  2, -2],
        [5, -2, -3,  3]])

M

Out[28]:

$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

In [29]:

M.eigenvals()

Out[29]:

$\left \{ -2 : 1, \quad 3 : 1, \quad 5 : 2\right \}$

This means that M has eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.

To find the eigenvectors of a matrix, use eigenvects. eigenvects returns a list of tuples of the form (eigenvalue:algebraic multiplicity, [eigenvectors]).

In [30]:

M.eigenvects()

Out[30]:

$\left [ \left ( -2, \quad 1, \quad \left [ \left[\begin{matrix}0\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 5, \quad 2, \quad \left [ \left[\begin{matrix}1\\1\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\-1\\0\\1\end{matrix}\right]\right ]\right )\right ]$

This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors. Because the algebraic and geometric multiplicities are the same for all the eigenvalues, M is diagonalizable.

To diagonalize a matrix, use diagonalize. diagonalize returns a tuple $(P,D)$ , where $D$ is diagonal and $M=PDP^{−1}$ .

In [31]:

P, D = M.diagonalize()

In [32]:

Out[32]:

$\left[\begin{matrix}0 & 1 & 1 & 0\\1 & 1 & 1 & -1\\1 & 1 & 1 & 0\\1 & 1 & 0 & 1\end{matrix}\right]$

In [33]:

Out[33]:

$\left[\begin{matrix}-2 & 0 & 0 & 0\\0 & 3 & 0 & 0\\0 & 0 & 5 & 0\\0 & 0 & 0 & 5\end{matrix}\right]$

In [34]:

P * D * P.inv()

Out[34]:

$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

In [35]:

P * D * P.inv() == M

Out[35]:

True

Note that since eigenvects also includes the eigenvalues, you should use it instead of eigenvals if you also want the eigenvectors. However, as computing the eigenvectors may often be costly, eigenvals should be preferred if you only wish to find the eigenvalues.

If all you want is the characteristic polynomial, use charpoly. This is more efficient than eigenvals, because sometimes symbolic roots can be expensive to calculate.

In [36]:

lamda = symbols('lamda')
p = M.charpoly(lamda)
factor(p)

Out[36]:

$\left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)$

Note: lambda is a reserved keyword in Python, so to create a Symbol called λ, while using the same names for SymPy Symbols and Python variables, use lamda (without the b). It will still pretty print as λ.

Non-square matrices don’t have eigenvectors and therefore don’t have an eigendecomposition. Instead, we can use the singular value decomposition to break up a non-square matrix A into left singular vectors, right singular vectors, and a diagonal matrix of singular values. Use the singular_values method on any matrix to find its singular values.

In [37]:

Out[37]:

$\left[\begin{matrix}1 & -1 & -1\\0 & 1 & 0\\1 & -2 & 1\end{matrix}\right]$

In [38]:

A.singular_values()

Out[38]:

$\left [ \sqrt{\sqrt{14} + 4}, \quad \sqrt{- \sqrt{14} + 4}, \quad \sqrt{2}\right ]$

References¶

SymPy Development Team (2016). Sympy Tutorial: Matrices
Ivan Savov (2016). Taming math and physics using SymPy

The following cell change the style of the notebook.

In [39]:

from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()

Out[39]:

In [ ]: