Linear Algebra with SciPy

The main Python package for linear algebra is the SciPy subpackage scipy.linalg which builds on NumPy. Let's import both packages:

import numpy as np
import scipy.linalg as la

NumPy Arrays

Let's begin with a quick review of NumPy arrays. We can think of a 1D NumPy array as a list of numbers. We can think of a 2D NumPy array as a matrix. And we can think of a 3D array as a cube of numbers.

When we select a row or column from a 2D NumPy array, the result is a 1D NumPy array (called a slice). This is different from MATLAB where when you select a column from a matrix it's returned as a column vector which is a 2D MATLAB matrix.

It can get a bit confusing and so we need to keep track of the shape, size and dimension of our NumPy arrays.

Array Attributes

Create a 1D (one-dimensional) NumPy array and verify its dimensions, shape and size.

a = np.array([1,3,-2,1])
print(a)
[ 1  3 -2  1]

Verify the number of dimensions:

a.ndim
1

Verify the shape of the array:

a.shape
(4,)

The shape of an array is returned as a Python tuple. The output in the cell above is a tuple of length 1. And we verify the size of the array (ie. the total number of entries in the array):

a.size
4

Create a 2D (two-dimensional) NumPy array (ie. matrix):

M = np.array([[1,2],[3,7],[-1,5]])
print(M)
[[ 1  2]
 [ 3  7]
 [-1  5]]

Verify the number of dimensions:

M.ndim
2

Verify the shape of the array:

M.shape
(3, 2)

Finally, verify the total number of entries in the array:

M.size
6

Select a row or column from a 2D NumPy array and we get a 1D array:

col = M[:,1] 
print(col)
[2 7 5]

Verify the number of dimensions of the slice:

col.ndim
1

Verify the shape and size of the slice:

col.shape
(3,)
col.size
3

When we select a row of column from a 2D NumPy array, the result is a 1D NumPy array. However, we may want to select a column as a 2D column vector. This requires us to use the reshape method.

For example, create a 2D column vector from the 1D slice selected from the matrix M above:

print(col)
[2 7 5]
column = np.array([2,7,5]).reshape(3,1)
print(column)
[[2]
 [7]
 [5]]

Verify the dimensions, shape and size of the array:

print('Dimensions:', column.ndim)
print('Shape:', column.shape)
print('Size:', column.size)
Dimensions: 2
Shape: (3, 1)
Size: 3

The variables col and column are different types of objects even though they have the "same" data.

print(col)
[2 7 5]
print('Dimensions:',col.ndim)
print('Shape:',col.shape)
print('Size:',col.size)
Dimensions: 1
Shape: (3,)
Size: 3

Matrix Operations and Functions

Arithmetic Operations

Recall that arithmetic array operations +, -, /, * and ** are performed elementwise on NumPy arrays. Let's create a NumPy array and do some computations:

M = np.array([[3,4],[-1,5]])
print(M)
[[ 3  4]
 [-1  5]]
M * M
array([[ 9, 16],
       [ 1, 25]])

Matrix Multiplication

We use the @ operator to do matrix multiplication with NumPy arrays:

M @ M
array([[ 5, 32],
       [-8, 21]])

Let's compute $2I + 3A - AB$ for

$$ A = \begin{bmatrix} 1 & 3 \\ -1 & 7 \end{bmatrix} \ \ \ \ B = \begin{bmatrix} 5 & 2 \\ 1 & 2 \end{bmatrix} $$ and $I$ is the identity matrix of size 2:

A = np.array([[1,3],[-1,7]])
print(A)
[[ 1  3]
 [-1  7]]
B = np.array([[5,2],[1,2]])
print(B)
[[5 2]
 [1 2]]
I = np.eye(2)
print(I)
[[1. 0.]
 [0. 1.]]
2*I + 3*A - A@B
array([[-3.,  1.],
       [-5., 11.]])

Matrix Powers

There's no symbol for matrix powers:

from numpy.linalg import matrix_power as mpow
mpow(M,2)
array([[2, 3],
       [3, 5]])
mpow(M,5)
array([[34, 55],
       [55, 89]])
M@M@M@M@M
array([[34, 55],
       [55, 89]])
M@M@M
array([[ 5,  8],
       [ 8, 13]])
mpow(M,3)
array([[ 5,  8],
       [ 8, 13]])

Tranpose

We can take the transpose with .T method:

M.T
array([[ 3, -1],
       [ 4,  5]])

Notice that $M M^T$ is a symmetric matrix:

M @ M.T
array([[25, 17],
       [17, 26]])

Inverse

We can find the inverse using scipy.linalg.inv:

A = np.array([[1,2],[3,4]])
print(A)
[[1 2]
 [3 4]]
la.inv(A)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

Trace

We can find the trace of a matrix using numpy.trace:

np.trace(A)
5

Norm

Under construction

Determinant

We find the determinant using scipy.linalg.det:

A = np.array([[1,2],[3,4]])
print(A)
[[1 2]
 [3 4]]
la.det(A)
-2.0

Dot Product

Examples

Characteristic polynomials and Cayley-Hamilton Theorem

The characteristic polynomial of a 2 by 2 square matrix $A$ is

$$ p_A(\lambda) = \det(A - \lambda I) = \lambda^2 - \mathrm{tr}(A) \lambda + \mathrm{det}(A) $$

The Cayley-Hamilton Theorem states that any square matrix satisfies its characteristic polynomial. For a matrix $A$ of size 2, this means that

$$ p_A(A) = A^2 - \mathrm{tr}(A) A + \mathrm{det}(A) I = 0 $$

Let's verify the Cayley-Hamilton Theorem for a few different matrices.

A
array([[1, 2],
       [3, 4]])
trace_A = np.trace(A)
det_A = la.det(A)
I = np.eye(2)
A @ A - trace_A * A + det_A * I
array([[0., 0.],
       [0., 0.]])

Let's do this again for some random matrices:

N = np.random.randint(0,10,[2,2])
print(N)
[[1 8]
 [6 1]]
trace_N = np.trace(N)
det_N = la.det(N)
I = np.eye(2)
N @ N - trace_N * N + det_N * I
array([[0., 0.],
       [0., 0.]])

Projections

The formula to project a vector $v$ onto a vector $w$ is

$$ \mathrm{proj}_w(v) = \frac{v \cdot w}{w \cdot w} w $$

Let's a function called proj which computes the projection $v$ onto $w$.

def proj(v,w):
    '''Project vector v onto w.'''
    v = np.array(v)
    w = np.array(w)
    return np.sum(v * w)/np.sum(w * w) * w   # or (v @ w)/(w @ w) * w

proj([1,2,3],[1,1,1])
array([2., 2., 2.])

Exercises

Exercise. Write a function which takes an input parameter $A$, $i$ and $j$ and returns the dot product of the $i$th and $j$th row (indexing starts at 0).

Exercise. Compute the matrix equation $AB + 2B^2 - I$ for matrices $A = \begin{bmatrix} 3 & 4 \\ -1 & 2 \end{bmatrix}$ and $B = \begin{bmatrix} 5 & 2 \\ 8 & -3 \end{bmatrix}$.