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 and so we must import the function matrix_power from the subpackage numpy.linalg.
from numpy.linalg import matrix_power as mpow
M = np.array([[3,4],[-1,5]])
print(M)
[[ 3 4]
[-1 5]]
mpow(M,2)
array([[ 5, 32],
[-8, 21]])
mpow(M,5)
array([[-1525, 3236],
[ -809, 93]])
Compare with the matrix multiplcation operator:
M @ M @ M @ M @ M
array([[-1525, 3236],
[ -809, 93]])
mpow(M,3)
array([[-17, 180],
[-45, 73]])
M @ M @ M
array([[-17, 180],
[-45, 73]])
Tranpose
We can take the transpose with .T attribute:
print(M)
[[ 3 4]
[-1 5]]
print(M.T)
[[ 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 the function 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 the function numpy.trace:
np.trace(A)
5
Norm
Under construction
Determinant
We find the determinant using the function scipy.linalg.det:
A = np.array([[1,2],[3,4]])
print(A)
[[1 2]
[3 4]]
la.det(A)
-2.0
Dot Product
Under construction
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.
print(A)
[[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 9]
[4 3]]
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 write 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 1. 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 2. Compute the matrix equation $AB + 2B^2 - I$ for matrices
$$A = \begin{bmatrix} 3 & 4 \\ -1 & 2 \end{bmatrix}\hspace{20mm} B = \begin{bmatrix} 5 & 2 \\ 8 & -3 \end{bmatrix}$$