Matrices and Vectors

This notebook discusses how to work with vectors and matrices in Numpy.

We will loosely follow LL § 2.3. Another good reference is the Numpy Documentation, especially the

In [1]:
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv

Numpy Arrays are Vectors

Vectors are just lists of numbers you can add or subtract. Numpy arrays serve the purpose of vectors.

In [2]:
v = np.array([1, 2, 5])
v
Out[2]:
array([1, 2, 5])

You can get the length of a numpy array, which is the same as the vectors dimension.

In [3]:
len(v)
Out[3]:
3

You can access the entries using square brackets. See below:

In [4]:
for i in range(len(v)):
    print("v[{}] = {}".format(i, v[i]))
v[0] = 1
v[1] = 2
v[2] = 5

You can also access entries from the back of the list. The last entry is v[-1].

In [5]:
for i in range(-1,-len(v)-1,-1):
    print("v[{}] = {}".format(i, v[i]))
v[-1] = 5
v[-2] = 2
v[-3] = 1

To be sure you have a numpy array you can check its type. Technically the arrays we are building are of type numpy.ndarray and array is just a convienient abbreviation.

In [6]:
type(v)
Out[6]:
numpy.ndarray

Some other functions defined by numpy also return numpy arrays. We have also seen:

In [7]:
x = np.linspace(1, 5, 4)
type(x)
Out[7]:
numpy.ndarray

You can also initialize a vector of all zeros.

In [8]:
z = np.zeros(3)
print("z = {}".format(z))
type(z)
z = [0. 0. 0.]
Out[8]:
numpy.ndarray

It is convientient sometimes to create a vector using zeros and then set its individual entries.

In [9]:
z = np.zeros(10)
for i in range(10):
    z[i] = m.sqrt(i)
print(z)
[0.         1.         1.41421356 1.73205081 2.         2.23606798
 2.44948974 2.64575131 2.82842712 3.        ]

The most important operations you can do with vectors are addition and scalar multiplication. Both work as you would expect.

In [10]:
print("v = {}".format(v))
0.7 * v
v = [1 2 5]
Out[10]:
array([0.7, 1.4, 3.5])
In [11]:
w = np.array([2, 0, -7])
print("v = {}".format(v))
print("w = {}".format(w))
v + w
v = [1 2 5]
w = [ 2  0 -7]
Out[11]:
array([ 3,  2, -2])

Another basic operation is the dot product.

In [12]:
v.dot(w)
Out[12]:
-33

Warning on copying

Note that setting one array equal to another one means the two names point to the same object. This can lead to unexpected issues like the following:

In [13]:
v = np.array([1,2,3])
w = v
w[0] = 5
print("v = {}".format(v))
print("w = {}".format(w))
v = [5 2 3]
w = [5 2 3]

The correct way to do this is the following:

In [14]:
v = np.array([1,2,3])
w = v.copy()
w[0] = 5
print("v = {}".format(v))
print("w = {}".format(w))
v = [1 2 3]
w = [5 2 3]

This is discussed in a bit more detail in § 2.3.4 of LL.

Arrays can store other types other than floats.

The type of number stored by an array v can be accessed with v.dtype.

In [15]:
v = np.array([1.1, 2.3])
print(v)
v.dtype
[1.1 2.3]
Out[15]:
dtype('float64')
In [16]:
v = np.array([0,1,2])
print(v)
v.dtype
[0 1 2]
Out[16]:
dtype('int64')

You can override the data type:

In [17]:
v = np.array([0,1,2], dtype=float)
print(v)
v.dtype
[0. 1. 2.]
Out[17]:
dtype('float64')
In [18]:
v = np.array([0,1,2**64])
print(v)
v.dtype
[0 1 18446744073709551616]
Out[18]:
dtype('O')

Storing multiple precision numbers:

In [19]:
mp.dps = 50
v = np.array([mp.tan(1), mp.sin(3)])
print(v)
v.dtype
[mpf('1.5574077246549022305069748074583601730872507723815206')
 mpf('0.1411200080598672221007448028081102798469332642522656')]
Out[19]:
dtype('O')

Storing intervals:

In [20]:
v = np.array([iv.tan(1), iv.pi])
print(v)
v.dtype
[mpi('1.5574077246549021', '1.5574077246549023')
 mpi('3.1415926535897931', '3.1415926535897936')]
Out[20]:
dtype('O')

Matrices

In [21]:
mat = np.array([[1,2,3],[4,5,6]])
print(mat)
mat.dtype
[[1 2 3]
 [4 5 6]]
Out[21]:
dtype('int64')

You can access the rows using m[row].

In [22]:
for row in range(len(mat)):
    print("row {} of mat is {}".format(row,mat[row]))
row 0 of mat is [1 2 3]
row 1 of mat is [4 5 6]

This means you can access the individual entries with m[row][col].

In [23]:
row = 0
col = 1
mat[row][col]
Out[23]:
2

Recall an $m \times n$ matrix is one with $m$ rows and $n$ columns. You can initialize a $m \times n$ matrix of all zeros with the following:

In [24]:
z = np.zeros((3,2))
z
Out[24]:
array([[0., 0.],
       [0., 0.],
       [0., 0.]])

Note that above, we pass a pair $(3,2)$ to zeros, which is why there are double parethesis. The following is the same:

In [25]:
dim = (3,2)
z = np.zeros(dim)
z
Out[25]:
array([[0., 0.],
       [0., 0.],
       [0., 0.]])

You can recover the number of rows and columns using shape:

In [26]:
z.shape
Out[26]:
(3, 2)

This means that we can use the following to iterate over the entries of any matrix.

In [27]:
num_rows, num_cols = mat.shape
for row in range(num_rows):
    for col in range(num_cols):
        print("mat[{}][{}] = {}".format(row, col, mat[row][col]))
mat[0][0] = 1
mat[0][1] = 2
mat[0][2] = 3
mat[1][0] = 4
mat[1][1] = 5
mat[1][2] = 6

Resizing matrices

You can change the shape of a matrix or a vector using the resize method. Below, we take a vector, resize it into a matrix, and then resize it back into a vector.

In [28]:
A = np.array([n for n in range(12)])
print("Initially A = {}.".format(A))
A.resize((3,4))
print("After first resize, A =\n{}.".format(A))
A.resize(12)
print("After second resize, A = {}.".format(A))
Initially A = [ 0  1  2  3  4  5  6  7  8  9 10 11].
After first resize, A =
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]].
After second resize, A = [ 0  1  2  3  4  5  6  7  8  9 10 11].

Alternately, you can use the reshape method to make a matrix with a new shape that shares the same underlying data. For example:

In [29]:
A = np.array([n for n in range(12)])
B = A.reshape((3,4))
print("A = {}.".format(A))
print("B = \n{}.".format(B))
A = [ 0  1  2  3  4  5  6  7  8  9 10 11].
B = 
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]].

Because $B$ shares data with $A$, both objects change when you change a value in either one. Observe:

In [30]:
A[11]=999
B[0][0]=888
print("A = {}.".format(A))
print("B = \n{}.".format(B))
A = [888   1   2   3   4   5   6   7   8   9  10 999].
B = 
[[888   1   2   3]
 [  4   5   6   7]
 [  8   9  10 999]].

So, if you want to preserve the original matrix, you need to make a copy.

In [31]:
A = np.array([n for n in range(12)])
B = A.reshape((3,4)).copy()
A[11]=999
B[0][0]=888
print("A = {}.".format(A))
print("B = \n{}.".format(B))
A = [  0   1   2   3   4   5   6   7   8   9  10 999].
B = 
[[888   1   2   3]
 [  4   5   6   7]
 [  8   9  10  11]].

Other constructors

The expression np.arrange(start, stop, step) is equivalent to np.array(range(start, stop, step) and default arguments work in the same way. You can then use the reshape method to adjust the dimensions.

In [32]:
A = np.arange(1,19,3)
print("First, A = {}.".format(A))
A.resize((2,3))
print("Now A is\n{}.".format(A))
First, A = [ 1  4  7 10 13 16].
Now A is
[[ 1  4  7]
 [10 13 16]].

A random $2 \times 1$ matrix with entries in $[0,1)$:

In [33]:
r = np.random.random((2,1))
print(r)
[[0.59458428]
 [0.35999069]]

The $3 \times 3$ identity matrix.

In [34]:
I = np.identity(3)
print(I)
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Joining matrices (or block matrices)

Arrays representing matrices can be joined with the function np.concatenate, but you need to be sure that the dimensions match. By default arrays are joined vertically. For example:

In [35]:
A = np.arange(4).reshape(2,2)
print("A =\n{}".format(A))
B = np.arange(10,14).reshape(2,2)
print("B =\n{}".format(B))
C = np.concatenate((A,B)) # Notice we are passing a pair
C
A =
[[0 1]
 [2 3]]
B =
[[10 11]
 [12 13]]
Out[35]:
array([[ 0,  1],
       [ 2,  3],
       [10, 11],
       [12, 13]])

Concatenate has an optional axis parameter that allows the matrices to be joined horizontally too. The above is equivalent to:

In [36]:
np.concatenate((A,B), axis=0)
Out[36]:
array([[ 0,  1],
       [ 2,  3],
       [10, 11],
       [12, 13]])

If we change to axis=1 we get horizontal joining.

In [37]:
np.concatenate((A,B), axis=1)
Out[37]:
array([[ 0,  1, 10, 11],
       [ 2,  3, 12, 13]])

Basic operations

Scalar multiplication:

In [38]:
mat = np.array([[1,2,3],[4,5,6]])
-2*mat
Out[38]:
array([[ -2,  -4,  -6],
       [ -8, -10, -12]])

Many numpy functions are applied coordinate-wise:

In [39]:
mat = np.array([[1,2,3],[4,5,6]])
np.sqrt(mat)
Out[39]:
array([[1.        , 1.41421356, 1.73205081],
       [2.        , 2.23606798, 2.44948974]])

This includes addition and multiplication:

In [40]:
mat = np.array([[1,2,3],[4,5,6]])
print(mat+mat)
[[ 2  4  6]
 [ 8 10 12]]
In [41]:
mat = np.array([[1,2,3],[4,5,6]])
print(mat*mat)
[[ 1  4  9]
 [16 25 36]]

Example of Random matrices:

Here is a more elaborate example. Suppose you want a $3 \times 3$ matrix with entries randomly chosen from $\{1, \ldots, 10\}$.

In [42]:
# First get a random matrix with entries in [0,1):
r = np.random.random((3,3))
print(r)
[[0.72104802 0.98360577 0.01470337]
 [0.32497131 0.26159116 0.55787085]
 [0.39488839 0.59512844 0.58600399]]
In [43]:
# The map x -> 10*(1-x) will be applied coordinatewise
# It puts entries into the interval $(0,10]$.
r = 10*(1-r)
print(r)
[[2.7895198  0.16394227 9.85296629]
 [6.75028688 7.38408839 4.42129147]
 [6.05111612 4.0487156  4.13996013]]
In [44]:
# The ceil function takes a number to the smallest integer 
# greater than that number. We apply it coordinatewise.
r = np.ceil(r)
r
Out[44]:
array([[ 3.,  1., 10.],
       [ 7.,  8.,  5.],
       [ 7.,  5.,  5.]])
In [45]:
# Observe r has data type float:
r.dtype
Out[45]:
dtype('float64')
In [46]:
# if we want integer we can do 
r = np.array(r, dtype=int)
r
Out[46]:
array([[ 3,  1, 10],
       [ 7,  8,  5],
       [ 7,  5,  5]])

Note that even operators behave coordinate wise. So, below we get an array off all booleans:

In [47]:
b = (r > 2)
print(b)
[[ True False  True]
 [ True  True  True]
 [ True  True  True]]

You can check whether all the entries are true with the all method. (This is a logical and of all the entries.)

In [48]:
b.all()
Out[48]:
False

You can check whether any of the entries are true with the any method. (This is a logical or of all the entries.)

In [49]:
b.any()
Out[49]:
True

We can get the total number of trues using the following:

In [50]:
b.sum()
Out[50]:
8

(When we add, True becomes $1$ and False becomes $0$.)

Matrix multiplication

Matrix multiplication can be carried out with the @ operator.

In [51]:
m1 = np.array([[3,0],[1,1]])
print(m1, end="\n\n")
m2 = np.array([[1,2],[0,1]])
print(m2, end="\n\n")
m3 = m1 @ m2
print(m3)
[[3 0]
 [1 1]]

[[1 2]
 [0 1]]

[[3 6]
 [1 3]]

It is equivalent to use the dot method:

In [52]:
m1 = np.array([[3,0],[1,1]])
print(m1, end="\n\n")
m2 = np.array([[1,2],[0,1]])
print(m2, end="\n\n")
m3 = m1.dot(m2)
print(m3)
[[3 0]
 [1 1]]

[[1 2]
 [0 1]]

[[3 6]
 [1 3]]

Transpose

Recall that you get the transpose of a matrix by switching rows and columns.

In [53]:
mat = np.array([n for n in range(12)])
mat.resize((3,4))
mat
Out[53]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [54]:
mat.transpose()
Out[54]:
array([[ 0,  4,  8],
       [ 1,  5,  9],
       [ 2,  6, 10],
       [ 3,  7, 11]])

Rows and columns

In [55]:
mat = np.array([n for n in range(12)])
mat.resize((3,4))
mat
Out[55]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

We have seen that a row can be accessed as vectors with mat[row]. We can also use mat[row,:]. Here the : indicates "all."

In [56]:
# First row:
mat[0, :]
Out[56]:
array([0, 1, 2, 3])

So, to get a column, we can do mat[:, col].

In [57]:
# Second column:
mat[:, 1]
Out[57]:
array([1, 5, 9])

Note that the above give them as a vectors (1-dimensional arrays). You can also access them as submatrices.

In [58]:
# First row:
mat[[0],:]
Out[58]:
array([[0, 1, 2, 3]])
In [59]:
# First column:
mat[:,[0]]
Out[59]:
array([[0],
       [4],
       [8]])

Submatrices

In [60]:
mat = np.array([n for n in range(16)])
mat.resize((4,4))
print(mat)
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

The second method of getting rows and columns is suggestive of a notation for submatrices. The following gets the first two columns as a matrix.

In [61]:
mat[:,[0,1]]
Out[61]:
array([[ 0,  1],
       [ 4,  5],
       [ 8,  9],
       [12, 13]])

Note that they don't have to be in order. The first, third and fourth rows in reverse order:

In [62]:
mat[[3,2,0],:]
Out[62]:
array([[12, 13, 14, 15],
       [ 8,  9, 10, 11],
       [ 0,  1,  2,  3]])

To get the submatrix consiting of those entries in rows $3$, $2$, or $0$ and in columns $0$ or $1$, you can combine the expressions above:

In [63]:
mat[[3,2,0],:][:,[0,1]]
Out[63]:
array([[12, 13],
       [ 8,  9],
       [ 0,  1]])

Slice notation for submatrices

The notation for accessing rows and columns is somewhat similar to the range command, in that you can give a pair of values start : end. Recall that in range this gives a list of values including start up to but not including end.

Here is a simple matrix for playing with these things.

In [64]:
mat = np.array([n for n in range(16)])
mat.resize((4,4))
print(mat)
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

Here we et the submatrix of elements in positions $(m,n)$ where $m,n \in \{1,2\}$. This is the central square in our matrix.

In [65]:
mat[1:3,1:3]
Out[65]:
array([[ 5,  6],
       [ 9, 10]])

If you omit the first element of a pair, it defaults to zero. So, the following gives the submatrix of entries lying in the first three rows and the first two columns.

In [66]:
mat[:3,:2]
Out[66]:
array([[0, 1],
       [4, 5],
       [8, 9]])

Defaults for the second entry are the total number of rows and the total number of columns, respectively. The following gives elements in rows $2$ and $3$ and in columns $1$, $2$ and $3$.

In [67]:
mat[2:,1:]
Out[67]:
array([[ 9, 10, 11],
       [13, 14, 15]])

You can specify just a single number with no colon to restrict attention to a row or column. This gives a vector consisting of the entries in row $2$ and in columns $1,\ldots,3$.

In [68]:
mat[2, 1:]
Out[68]:
array([ 9, 10, 11])

You can also use a triple of values to specify start : end : step for the rows and the columns. The following gives the matrix cosisting of entries in positions $(m,n)$ with $m$ and $n$ even.

In [69]:
print(mat, end="\n\n")
mat[::2,::2]
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

Out[69]:
array([[ 0,  2],
       [ 8, 10]])

The entries in odd positions:

In [70]:
mat[1::2,1::2]
Out[70]:
array([[ 5,  7],
       [13, 15]])

Here we reverse the order of rows and columns:

In [71]:
mat[::-1, ::-1]
Out[71]:
array([[15, 14, 13, 12],
       [11, 10,  9,  8],
       [ 7,  6,  5,  4],
       [ 3,  2,  1,  0]])

Assignments to submatrices

When you get submatrices, the data is passed by reference. This means assignments to submatrices change the original matrix.

In [72]:
mat = np.array([n for n in range(16)]).reshape((4,4))
print(mat)
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

Here we set all entries in even positions to $-1$.

In [73]:
mat[::2,::2]=-1
print(mat)
[[-1  1 -1  3]
 [ 4  5  6  7]
 [-1  9 -1 11]
 [12 13 14 15]]

Here we double all entries in odd positions:

In [74]:
mat[1::2,1::2] *= 2
print(mat)
[[-1  1 -1  3]
 [ 4 10  6 14]
 [-1  9 -1 11]
 [12 26 14 30]]

You can also assign a matrix to assign different values to each entry.

In [75]:
mat[::2,::2] = np.array([[-2, -3], [-4, -5]])
print(mat)
[[-2  1 -3  3]
 [ 4 10  6 14]
 [-4  9 -5 11]
 [12 26 14 30]]

This is particularly useful for row reduction. As a first step in row reduction, we would clear the entries in the first column below the $-2$ in the first row by adding multiples of the first row. We can do this with:

In [76]:
mat[1,:] += 2 * mat[0,:]
mat[2,:] += -2 * mat[0,:]
mat[3,:] += 6 * mat[0,:]
print(mat)
[[-2  1 -3  3]
 [ 0 12  0 20]
 [ 0  7  1  5]
 [ 0 32 -4 48]]

Swapping rows and columns

Recall that to get columns $1$ and $2$ we can do:

In [77]:
mat[:,[1,2]]
Out[77]:
array([[ 1, -3],
       [12,  0],
       [ 7,  1],
       [32, -4]])

To get them in the opposite order we can do:

In [78]:
mat[:,[2,1]]
Out[78]:
array([[-3,  1],
       [ 0, 12],
       [ 1,  7],
       [-4, 32]])

So, to swap these two columns you can do mat[:,[1,2]] = mat[:,[2,1]]:

In [79]:
print(mat, end="\n\n")
mat[:,[1,2]] = mat[:,[2,1]]
print(mat)
[[-2  1 -3  3]
 [ 0 12  0 20]
 [ 0  7  1  5]
 [ 0 32 -4 48]]

[[-2 -3  1  3]
 [ 0  0 12 20]
 [ 0  1  7  5]
 [ 0 -4 32 48]]