Skip to content

Vectors & Matrices (NumPy)

Let's begin by importing the NumPy package.

import numpy as np

NumPy includes some tools for working with linear algebra in the numpy.linalg module. However, unless you really don’t want to add SciPy as a dependency to your project, it’s typically better to use scipy.linalg for the following reason:

  • As explained in the official documentation, scipy.linalg contains all the functions in numpy.linalg plus some extra advanced functions not included in numpy.linalg.

Therefore, when we need to use linear algebra specific functions we'll load the scipy.linalg library. The first instance where we do this below is when we compute matrix inverses.

NumPy Arrays

Vectors and Matrices are created as instances of a numpy array. We can think of a 1D NumPy array as a list of numbers (or row-vector), and a 2D number array as a matrix.

1D Arrays

Let's create a 1D array of integers.

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

We can verify the dimension, shape and size.

print('Dimension:', a.ndim) # (1)
print('Shape:', a.shape) # (2)
print('Size:', a.size) # (3)
  1. ndim returns the dimension of the array. This is how many indices are needed to specify an entry.
  2. shape returns a tuple of ndim numbers. For a 2D array this is (rows, columns).
  3. size returns the total number of entries in the array.
Dimension: 1
Shape: (4,)
Size: 4

Therefore, since ndim is 1 this is a 1-dimensional array. The shape method returns a python tuple of length 1 - this tells us the array is 1D and has 4 entries. The size method also tells us the array has 4 entries.

2D Arrays

Let's create a 2D array.

A=np.array([[1,2,3],[4,5,6]])

array([[1, 2, 3],
       [4, 5, 6]])

The array was entered as a list of rows.

We can verify the dimension of the array. The following tells us it is 2-dimensional, therefore to select an entry in the array we need 2 indices to locate it (row number and column number). We can also verify the shape and size too. 

print('Dimension:', A.ndim) # (1)
print('Shape:', A.shape) # (2)
print('Size:', A.size) # (3)

  1. ndim returns the dimension of the array. This is how many indices are needed to specify an entry.
  2. shape returns a tuple of ndim numbers. For a 2D array this is (rows, columns).
  3. size returns the total number of entries in the array.
Dimension: 2
Shape: (2,3)
Size: 6

.reshape

Another way to create a 2D array is to first create a 1D array and then use the .reshape method to change its shape as desired.

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

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

A.reshape(4,3) # (1)
  1. This does not change the array stored in the variable A, it creates a new instance that has been reshaped.
array([[ 1,  2,  3],
       [ 4,  5,  6],
       [ 7,  8,  9],
       [10, 11, 12]])

.mat

NumPy has a built in matrix object, which is really just a 2D array. The only difference in using .mat as opposed to .array is that the operator * can be used for the usual matrix product. More on this below.

We can create a matrix directly as follows.

A = np.mat([[1,2,3],[4,5,6]])
A

matrix([[1, 2, 3],
        [4, 5, 6]])

We can also take an existing 2D array and convert it to a matrix with asmatrix.

A=np.array([1,1,1,2,2,2,3,3,3]).reshape(3,3)
A=np.asmatrix(A)
A

matrix([[1, 1, 1],
        [2, 2, 2],
        [3, 3, 3]])

There isn't much of an advantage using .mat over .array for 2D arrays so we'll always use an array. Just so you are aware of some shortcuts available when working with .mat objects, we have .T for transpose, .H for Hermetian, and .I for inverse. We also have * for matrix multiplication. However, we have similar methods for .array objects which we will cover below.

Indexing

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

a[0] # recall a = [1,2,3,4]

1

We can also use index slicing to select a subarray.

a[1:3] # (1)
  1. select entries from index 1 (inclusive) to 3 (non-inclusive). That is, select entries with index from 1 to 2.
array([2,3])

This can be done for 2D arrays by specifying both indices.

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

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

A[3,2] # (1)
  1. This selects the entry in row indexed by 3, and column indexed by 2. Remember indexing starts at 0.
12

We can use slicing to select submatricies.

A[0:2,0:3] # (1)
  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.
array([[1, 2, 3],
       [4, 5, 6]])

An entire column can be selected using slicing.

b = A[:,0] # (1)
b
  1. : means run over all values of the index. Therefore, this selects the first column.
array([ 1,  4,  7, 10])

Notice that b looks to be a 1D array (rather than column vector, which is a 2D array of shape 4x1). We can verify this by checking the dimension.

b.ndim

1

Unfortunately, even though we selected a column of the array this returned a 1D array, which is essentially a row vector. Given the circumstance we may have wanted this to be a column vector. We can reshape it into a 4x1 column vector.

b = b.reshape(4,1) # (1)
b
  1. We change the shape of b and then rewrite the new data to the same variable name.
array([[ 1],
       [ 4],
       [ 7],
       [10]])

print('Dimension:', b.ndim)
print('Shape:', b.shape)
print('Size:', b.size)

Dimension: 2
Shape: (4, 1)
Size: 4

Now we have successfully extracted the first column as a column vector. We could do this in one line as follows:

col = A[:,0].reshape(4,1)
col

array([[ 1],
       [ 4],
       [ 7],
       [10]])

Special Matrices

NumPy 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).

np.eye(3)

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Zeros

An nxm zero matrix is constructed with zeros([n,m]).

np.zeros([3,2])

array([[0., 0.],
       [0., 0.],
       [0., 0.])

Ones

An nxm matrix of all 1's is constructed with ones([n,m]).

np.ones([3,2])

array([[1., 1.],
       [1., 1.],
       [1., 1.])

Example

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

\[ A = \begin{bmatrix} 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{bmatrix} \]
A = np.zeros([5,5])
A[1:,:-1] = np.eye(4) # lower diagonal
print(A)

[[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.]]

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 = np.zeros([5,5])
AL[1:,:-1] = np.eye(4) # lower diagonal
AU = np.zeros([5,5])
AU[:-1,1:] = np.eye(4) # upper diagonal
A = AL + AU
print(A)

[[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.]]

diagonal and block

To create a diagonal matrix use the NumPy diag function.

np.diag([1,2,3])

array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])

To create a block matrix use the NumPy block function.

B2 = np.array([[2,2],[2,2]])
B3 = np.array([[3,3,3],[3,3,3],[3,3,3]])
np.block([[B2,np.zeros([2,3])],
          [np.zeros([3,2]),B3]])

array([[2., 2., 0., 0., 0.],
       [2., 2., 0., 0., 0.],
       [0., 0., 3., 3., 3.],
       [0., 0., 3., 3., 3.],
       [0., 0., 3., 3., 3.]])

random

Use np.random.rand to create an array of the given shape and populate it with random samples from a uniform distribution over [0, 1).

A = np.random.rand(3,4)
print(A)

[[0.40104298 0.25146289 0.94189368 0.320897  ]
 [0.20052053 0.25960576 0.18908785 0.07176892]
 [0.29962823 0.31623458 0.40490173 0.99967804]]

To generate a matrix with random integer values use numpy.random.randint. The first two arguments are the lower and upper bounds on the entries.

np.random.randint(-10,20, size=(2, 4))

array([[-4,  2, -9, 12],
       [ 0,  6, -5,  6]])

Vector Operations

Dot Product

The dot product of two vectors (constructed as 1D arrays) can be done with the numpy.dot function. We also give two alternate methods as well.

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

v = np.array([1,2,3])
w = np.array([1,2,1])
sum(v*w)

v = np.array([1,2,3])
w = np.array([1,2,1])
v@w

8

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} \]

This function isn't built directly into NumPy, but numpy.dot and numpy.sqrt are so we can use those.

v = np.array([3,4])
np.sqrt(np.dot(v,v))

5.0

If you plan on computing the length frequently, you may consider just defining a function.

def norm(v):
    return np.sqrt(np.dot(v,v))

def norm(v):
    return np.sqrt(sum(v*v))


v = np.array([3,4])
norm(v)

5.0

This is how I would do it if I wasn't already importing the SciPy package. However, if we are using the SciPy package already, then we could call scipy.linalg.norm to compute the length of a vector.

import scipy.linalg as la

v = np.array([3,4])
la.norm(v)

5.0

Angle between two vectors

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

\[ \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{|a| |b|} \]
a = np.array([1,0])
b = np.array([1,1])

np.arccos( (np.dot(a,b)) / (np.sqrt(np.dot(a,a)*np.dot(b,b)) ) )

0.7853981633974484

Notice, this is just the decimal value of \(\pi/4\).

np.pi/4

0.7853981633974483

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

def angle(a,b):
    return(np.arccos( (np.dot(a,b)) / (np.sqrt(np.dot(a,a)*np.dot(b,b)) ) ))

a = np.array([1,1])
b = np.array([-1,0])
angle(a,b)

2.356194490192345

This value is \(3\pi/4\).

Matrix Operations

Matrix Scalar Multiplication

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

array([[ 3,  6],
       [ 9, 12]])

Matrix Addition

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

array([[-1,  7],
       [10,  5]])

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 is @. 

A = np.array([1,2,3,4]).reshape(2,2)
B = np.array([-2,5,7,1]).reshape(2,2)
A @ B


array([[12,  7],
       [22, 19]])

Alternatively, we could use numpy.matmul.

np.matmul(A,B)

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 = np.array([1,2,3,4,5,6]).reshape(2,3)
D = np.array([1,2,1,2,1,2,1,2,1,2,1,2]).reshape(3,4)
C @ D

array([[ 6, 12,  6, 12],
       [15, 30, 15, 30]])

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 = np.array([1,2,2,1]).reshape(2,2)
B = np.array([1,1,1,1]).reshape(2,2)
I = np.eye(2)
3*I+2*A-A@B

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

Info

One may think to use * as a multiplication operator, but for NumPy arrays it multiplies entries elementwise.

a = np.array([2,3])
b = np.array(10,20])
a*b

array([20, 60])

This will do elementwise multiplication on any two matrices as long as they have the same shape.

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

array([[-2, 10],
       [21,  4]])

For .mat objects * will do the usual matrix multiplication. This used to be one major advantage for using .mat over .array, but with the @ operator for matrix multiplication on 2D arrays, .mat looses this advantage over .array.

Matrix Powers

To compute a matrix power we use the matrix_power function in the numpy.linalg library.

A = np.array([1,2,-1,1]).reshape(2,2)
np.linalg.matrix_power(A,2)

array([[-1,  4],
       [-2, -1]])

It is customary to import this in with the name mpow.

from numpy.linalg import matrix_power as mpow

A = np.array([1,2,-1,1]).reshape(2,2)
mpow(A,2)

array([[-1,  4],
       [-2, -1]])

Transpose

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

A = np.array([1,2,-1,1]).reshape(2,2)
print(A)

[[1 2]
 [3 4]]

print(A.T)

[[1 3]
 [2 4]]

This can equivalently be done as a function call.

print(np.transpose(A))

Inverse

To find the inverse of a matrix this is where we need to bring in the scipy library. In particular we need scipy.linalg.inv. To do this we'll just bring in the entire scipy.linalg library as la.

import scipy.linalg as la

A = np.array([1,2,3,4]).reshape(2,2)
print(A)

[[1 2]
 [3 4]]

print(la.inv(A))

[[-2.   1. ]
 [ 1.5 -0.5]]

Trace

We can find the trace of a matrix using the .trace function on an array.

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

5

This can also be done by calling the function numpy.trace.

np.trace(A)

Determinant

We can find the determinant using the scipy.linalg.det function.

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

-2

Rank

We can find the rank using the numpy.linalg.matrix_rank function.

A = np.array([[1,1,0],[1,1,1],[-1,-1,2]])
np.linalg.matrix_rank(A)

2

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 = np.trace(A)
    det_A = la.det(A)
    I = np.eye(2)
    return A@A - trace_A*A + det_A*I

Example

A = np.random.randint(0,10,[2,2])
print(A)

[[7 7]
[7 8]]


check_cayley_hamilton(A)

array([[0., 0.],
       [0., 0.]])

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 ((np.dot(v,w)/np.dot(w,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 = np.array([4,2])
w = np.array([1,1])
proj(v,w)

array([3., 3.])

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

perp(v,w)

array([ 1., -1.])

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

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

array([4., 2.])

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 = np.array([1,3,-1])
w = np.array([1,1,2])
proj(v,w)

array([0.33333333, 0.33333333, 0.66666667])

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

perp(v,w)

array([ 0.66666667,  2.66666667, -1.66666667])

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

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

array([ 1.,  3., -1.])