Vectors & Matrices (SymPy)

Let's begin by importing the SymPy package.

from sympy import *

The * means the namespace for all functions in SymPy are loaded, so we can call them directly without having to preface them with the package name.

SymPy Matrices

Vectors and Matrices are created as instances of a SymPy Matrix object. It is constructed by providing a list of row vectors that make up the matrix. Providing just a single list of entries creates a column vector. We'll get up and running with a number of examples below. So learn more check out the full SymPy Matrix documentation.

Vectors

Let's create a vector of integers.

a = Matrix([1,2,3,4])
a

\(\left[ \begin{matrix} 1 \\ 2 \\ 3 \\ 4 \\ \end{matrix} \right]\)

We can verify the shape.

print('Shape:', a.shape) # (1)

shape returns a pair of numbers (rows, columns).

Shape: (4,1)

Matrices

Let's create a matrix.

A = Matrix([[1,2,3],[4,5,6]])
A

\(\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ \end{matrix} \right]\)

The matrix was entered as a list of rows.

Alternatively, we can construct a matrix by first stating its size, then the list of entries.

A = Matrix(3,2,[1,2,3,4,5,6])
A

\(\left[ \begin{matrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \\ \end{matrix} \right]\)

`.reshape`

Another way to create a matrix is to first create column vector and then use the .reshape method to change its shape as desired.

A = Matrix([1,2,3,4,5,6,7,8,9,10,11,12]).reshape(2,6)
A

\(\left[ \begin{matrix} 1 & 2 & 3 & 4 & 5 & 6 \\ 7 & 8 & 9 & 10 & 11 & 12\\ \end{matrix} \right]\)

A.reshape(4,3) # (1)

This does not change the matrix stored in the variable A, it creates a new instance that has been reshaped.

\(\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ 10 & 11 & 12\\ \end{matrix} \right]\)

Variable entries

Variables can be used in SymPy. First we need to declare them.

a,b,c,d = symbols('a:d') # shorthand for multiple symbols

A = Matrix([[a,2,3],[2*a,b+1,b*c+a]])
A

\(\left[ \begin{matrix} a & 2 & 3 \\ 2a & b+1 & bc+a \\ \end{matrix} \right]\)

Indexing

In Python indices begin at 0. Therefore, to select the first entry of vector a we use a[0].

a = Matrix([1,2,3,4])
a[0]

We can also use index slicing to select multiple entries in a matrix.

a[1:3] # (1)

select entries from index 1 (inclusive) to 3 (non-inclusive). That is, select entries with index from 1 to 2.

[2, 3]

This can be done for matrices by specifying both indices.

A = Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])
A

\(\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ 10 & 11 & 12 \\ \end{matrix} \right]\)

A[3,2] # (1)

This selects the entry in row indexed by 3, and column indexed by 2. Remember indexing starts at 0.

We can use slicing to select submatricies.

A[0:2,0:3] # (1)

This selects all entries from row 0 (inclusive) up to row 2 (not inclusive), and column 0 (inclusive) up to column 3 (not inclusive). This will be a 2x3 submatrix.

\(\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ \end{matrix} \right]\)

An entire column can be selected using slicing.

b = A[:,0] # (1)
b

: means run over all values of the index. Therefore, this selects the first column.

\(\left[ \begin{matrix} 1 \\ 4 \\ 7 \\ 10 \\ \end{matrix} \right]\)

Notice that b is a column vector.

Accessing Rows and Columns

Columns and rows can also be extracted using .col and .row.

A = Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])
A

\(\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ 10 & 11 & 12 \\ \end{matrix} \right]\)

A.col(0)

\(\left[ \begin{matrix} 1 \\ 4 \\ 7 \\ 10 \\ \end{matrix} \right]\)

A.row(2)

\(\left[ \begin{matrix} 7 & 8 & 9 \\ \end{matrix} \right]\)

Inserting and Deleting Rows and Columns

Columns and rows can be inserted using .col_insert and .row_insert. This does not change the original matrix.

Here we add a column of \(1\)'s to the end of matrix \(A\).

A.col_insert(A.shape[1]+1,ones(A.shape[0],1))

\(\left[ \begin{matrix} 1 & 2 & 3 & 1\\ 4 & 5 & 6 & 1\\ 7 & 8 & 9 & 1\\ 10 & 11 & 12 & 1\\ \end{matrix} \right]\)

Two matrices can be joined together using row_join or col_join. These names may be at first counterintuitive. The row_join command concatenates rows, therefore making more columns. Similar, col_join concatenates columns.

A = Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])
B = Matrix([[1,2],[1,2],[1,2],[1,2]])
A.row_join(B)

\(\left[ \begin{matrix} 1 & 2 & 3 & 1 & 2\\ 4 & 5 & 6 & 1 & 2\\ 7 & 8 & 9 & 1 & 2\\ 10 & 11 & 12 & 1 & 2\\ \end{matrix} \right]\)

Columns and rows can be deleted using .col_del and .row_del. This changes the original matrix.

A.col_del(2)
A

\(\left[ \begin{matrix} 1 & 2 1\\ 4 & 5 \\ 7 & 8 \\ 10 & 11\\ \end{matrix} \right]\)

Special Matrices

SymPy has some predefined special matrices which can be useful for quickly constructing matrices.

Identity

An \(n \times n\) identity matrix is constructed with eye(n). Can also use eye(m,n) for rectangular matrices.

eye(3)

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

Zeros

An \(n \times m\) zero matrix is constructed with zeros(n,m). Can also use zeros(n) for square matrices.

zeros(3,2)

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

Ones

An nxm matrix of all 1's is constructed with ones(n,m). Can also use ones(n) for square matrices.

ones(3,2)

\(\left[ \begin{matrix} 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ \end{matrix} \right]\)

Example

Construct the \(5 \times 5\) lower diagonal matrix.

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

A = zeros(5)
A[1:,:-1] = eye(4) # lower diagonal
print(A)

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

Example

Construct the \(5 \times 5\) matrix with lower and upper diagonals of \(1\)'s.

\[ A = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ \end{bmatrix} \]

AL = zeros(5)
AL[1:,:-1] = eye(4) # lower diagonal
AU = zeros(5)
AU[:-1,1:] = eye(4) # upper diagonal
A = AL + AU
A

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

diagonal and block

To create a diagonal matrix use the SymPy diag function.

diag(1,2,3)

\(\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{matrix}\right]\)

We can also use diag to create a block matrix.

B2 = Matrix([[2,2],[2,2]])
B3 = Matrix([[3,3,3],[3,3,3],[3,3,3]])
diag(1,B2,B3)

\(\left[\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 2 & 2 & 0 & 0 & 0\\ 0 & 2 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 3 & 3 & 3\\ 0 & 0 & 0 & 3 & 3 & 3\\ 0 & 0 & 0 & 3 & 3 & 3 \end{matrix}\right]\)

For further creation of block matrices lookup the BlockMatrix command in SymPy.

random

To generate a random matrix in SymPy we use the random package to first generate a list of random integers in some range. Then we reshape the list into a matrix.

import random
#Generate 8 random numbers between -10 and 10
randomlist = random.sample(range(-10, 10), 8)
A = Matrix(randomlist).reshape(2,4)
A

\(\left[ \begin{matrix} -2 & 0 & -9 & 4 \\ 8 & -5 & 7 & 3 \\ \end{matrix} \right]\)

Vector Operations

Dot Product

The dot product of two vectors can be done with the .dot function.

v = Matrix([1,2,3])
w = Matrix([1,2,1])
v.dot(w)

Length (Norm)

The length (also called the norm) of an \(n\)-dimensional vector \(\vec{v} = [v_1, v_2, \ldots , v_n]\) is defined by

\[ \begin{align} |\vec{v}| & = \sqrt{v_1^2+v_2^2 + \cdots + v_n^2} \\ & = \sqrt{\vec{v} \cdot \vec{v}} \end{align} \]

The norm of a vector can be done with the .norm function.

v = Matrix([3,4])
v.norm()

5.0

SymPy returns the exact value of the norm.

v = Matrix([1,2,3])
v.norm()

\(\sqrt{14}\)

If we want a decimal approximation we can use .evalf() or .n() or N() to covert to floating point (whichever you prefer). By default numerical evaluation is performed to an accuracy of 15 decimal digits. You can optionally pass a desired accuracy (a positive integer) as an argument to .evalf() or .n() or N(). Below we show all three functions, with both the default number of places, and with \(25\) decimal places.

.evalf().n()N()

v.norm().evalf()
v.norm().evalf(25)

w.norm().n()
w.norm().n(25)

N(w.norm())
N(w.norm(),25)

3.74165738677394
3.741657386773941385583749

Angle between two vectors

The angle \(\theta\) between two vectors satisfies

\[ \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{|a| |b|} \]

a = Matrix([1,0])
b = Matrix([1,1])

acos(a.dot(b)/(a.norm()*b.norm()))

\(\frac{\pi}{4}\)

We can build this into a function to make it easier to call.

def angle(a,b):
    return(acos( (a.dot(b)) / (a.norm()*b.norm()) ) )

a = Matrix([1,1])
b = Matrix([-1,0])
angle(a,b)

\(\frac{3\pi}{4}\)

Let's do an example where the angle isn't so nice.

a = Matrix([1,1/2])
b = Matrix([-3,2])
theta = angle(a,b)
print(theta)
print(theta.evalf())

acos(-0.137604183230756*sqrt(13))
2.08994244104142

Matrix Operations

Matrix Scalar Multiplication

A = Matrix([1,2,3,4]).reshape(2,2)
3*A

\(\left[ \begin{matrix} 3 & 6 \\ 9 & 12 \\ \end{matrix} \right]\)

Matrix Addition

A = Matrix([1,2,3,4]).reshape(2,2)
B = Matrix([-2,5,7,1]).reshape(2,2)
A+B

\(\left[ \begin{matrix} -1 & 7 \\ 10 & 5 \\ \end{matrix} \right]\)

Matrix Multiplication

Let's start with the two matrices.

\(A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\) and \(B = \begin{bmatrix} -2 & 5 \\ 7 & 1 \end{bmatrix}\)

The matix multiplication operator in SymPy is *.

A = Matrix([1,2,3,4]).reshape(2,2)
B = Matrix([-2,5,7,1]).reshape(2,2)
A * B

\(\left[ \begin{matrix} 12 & 7 \\ 22 & 19 \\ \end{matrix} \right]\)

Here is another example,

\(C = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}\) and \(D = \begin{bmatrix} 1 & 2 & 1 & 2 \\ 1 & 2 & 1 & 2 \\ 1 & 2 & 1 & 2 \end{bmatrix}\)

The matrix product is a 2x4 matrix.

C = Matrix([1,2,3,4,5,6]).reshape(2,3)
D = Matrix([1,2,1,2,1,2,1,2,1,2,1,2]).reshape(3,4)
C * D

\(\left[ \begin{matrix} 6 & 12 & 6 & 12 \\ 15 & 30 & 15 & 30 \\ \end{matrix} \right]\)

Example

Putting this all together, let's compute \(3I + 2A -AB\) for

\(A = \begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix}\) and \(B = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}\)

A = Matrix([1,2,2,1]).reshape(2,2)
B = Matrix([1,1,1,1]).reshape(2,2)
I = eye(2)
3*I+2*A-A*B

\(\left[ \begin{matrix} 2 & 1 \\ 1 & 2 \\ \end{matrix} \right]\)

Matrix Powers

To compute a matrix power we use the ** operator.

A = Matrix([1,2,-1,1]).reshape(2,2)
A**2

\(\left[ \begin{matrix} -1 & 4 \\ -2 & -1 \\ \end{matrix} \right]\)

Transpose

The transpose of a matrix can be taken by using either of these methods on an array: .T or .transpose.

A = Matrix([1,2,-1,1]).reshape(2,2)
A

\(\left[ \begin{matrix} 1 & 2 \\ 3 & 4 \\ \end{matrix} \right]\)

A.T

\(\left[ \begin{matrix} 1 & 3 \\ 2 & 4 \\ \end{matrix} \right]\)

Inverse

To find the inverse of a matrix raise it to the power of -1.

A = Matrix([1,2,3,4]).reshape(2,2)
A

\(\left[ \begin{matrix} 1 & 2 \\ 3 & 4 \\ \end{matrix} \right]\)

A**(-1)

\(\left[ \begin{matrix} -2 & 1 \\ 3/2 & -1/2 \\ \end{matrix} \right]\)

Trace

We can find the trace of a matrix by using .trace, either called as a method A.trace() or a function trace(A).

A = Matrix([1,2,3,4]).reshape(2,2)
A.trace()

Determinant

We can find the determinant of a matrix by using .det, either called as a method A.det() or a function det(A).

A = Matrix([1,2,3,4]).reshape(2,2)
A.det()

-2

Rank

We can find the rank using the .rank function.

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

Examples

Characteristic Polynomial and the Cayley-Hamilton Theorem

The characteristic polynomial of a 2x2 matrix \(A\) is

\[ p_A(\lambda) = \operatorname{det}(A-\lambda I) = \lambda^2 - \operatorname{tr}(A)\lambda + \operatorname{det}(A)\]

The Cayley-Hamilton Theorem states that any square matrix satisfies its characteristic polynomial. This means

\[ p_A(A) = A^2 - \operatorname{tr}(A)A + \operatorname{det}(A)I=0\]

Let's check that this is the case for some matrices.

def check_cayley_hamilton(A):
    # input: a 2x2 matrix A
    trace_A = A.trace()
    det_A = A.det()
    I = eye(2)
    return A*A - trace_A*A + det_A*I

Example

A = Matrix(random.sample(range(0,10),4)).reshape(2,2)
A

\(\left[ \begin{matrix} 0 & 1 \\ 7 & 5 \\ \end{matrix} \right]\)

check_cayley_hamilton(A)

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

Projections

The projection of a vector \(\vec{v}\) onto a vector \(\vec{w}\) is

\[\operatorname{proj}_{\vec{w}}(\vec{v}) = \left(\frac{\vec{v}\cdot \vec{w}}{\vec{w}\cdot \vec{w}}\right)\vec{w}\]

The perpendicular of the projection of \(\vec{v}\) onto \(\vec{w}\) is

\[\operatorname{perp}_{\vec{w}}(\vec{v}) = \vec{v} - \operatorname{proj}_{\vec{w}}(\vec{v})\]

This decomposes the vector \(\vec{w}\) into two parts:

\[\vec{v} = \operatorname{proj}_{\vec{w}}(\vec{v}) + \operatorname{perp}_{\vec{w}}(\vec{v})\]

where \(\operatorname{proj}_{\vec{w}}(\vec{v})\) is parallel to \(\vec{w}\) and \(\operatorname{perp}_{\vec{w}}(\vec{v})\) is perpenicular to \(\vec{w}\).

Let's define a function called proj which returns the projection of \(\vec{v}\) onto \(\vec{w}\), and a function called perp.

def proj(v,w):
    # input: v, w are vectors represented as 1D arrays
    # output: projection of v onto w
    return (v.dot(w)/w.dot(w))*w

def perp(v,w):
    # input: v, w are vectors represented as lists
    # output: perpendicular of the projection of v onto w  
    return v-proj(v,w)

Example

Consider the vectors \(\vec{v} = [4,2]\) and \(\vec{w} = [1,1]\). The projection of \(\vec{v}\) onto \(\vec{w}\) is \(3\vec{w} = [3,3]\).

v = Matrix([4,2])
w = Matrix([1,1])
proj(v,w)

\(\left[ \begin{matrix} 3 \\ 3 \\ \end{matrix} \right]\)

The perpendicular is \([1,-1]\):

perp(v,w)

\(\left[ \begin{matrix} 1 \\ -1 \\ \end{matrix} \right]\)

We can verify the sum is the vector \(\vec{v}\).

proj(v,w) + perp(v,w)

\(\left[ \begin{matrix} 4 \\ 2 \\ \end{matrix} \right]\)

Example

Consider the vectors \(\vec{v} = [1,3,-1]\) and \(\vec{w} = [1,1,2]\). The projection of \(\vec{v}\) onto \(\vec{w}\) is \(3\vec{w} = [1/3,1/3,2/3]\).

v = Matrix([1,3,-1])
w = Matrix([1,1,2])
proj(v,w)

\(\left[ \begin{matrix} 1/3 \\ 1/3 \\ 2/3 \\ \end{matrix} \right]\)

The perpendicular is \([2/3, 8/3, -5/3]\):

perp(v,w)

\(\left[ \begin{matrix} 2/3 \\ 8/3 \\ -5/3 \\ \end{matrix} \right]\)

We can check these sum to \(\vec{v}\):

proj(v,w) + perp(v,w)

\(\left[ \begin{matrix} 1 \\ 3 \\ -1 \\ \end{matrix} \right]\)