Estimated time to complete this notebook: 20 minutes

## 3.5.1 Recap#

In the previous section we introduced numpy array that represents a multidimensional matrix $$M_{i,j,k...n}$$. Which, among other things, allows for vectorised versions of common functions.

import numpy as np


This is another really powerful feature of NumPy.

By default, array operations are element-by-element:

np.arange(5) * np.arange(5)

array([ 0,  1,  4,  9, 16])


If we multiply arrays with non-matching shapes we get an error:

np.arange(5) * np.arange(6)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 np.arange(5) * np.arange(6)

ValueError: operands could not be broadcast together with shapes (5,) (6,)

np.zeros([2, 3]) * np.zeros([2, 4])

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 np.zeros([2, 3]) * np.zeros([2, 4])

ValueError: operands could not be broadcast together with shapes (2,3) (2,4)

m1 = np.arange(100).reshape([10, 10])

m2 = np.arange(100).reshape([10, 5, 2])

m1 + m2

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 1
----> 1 m1 + m2

ValueError: operands could not be broadcast together with shapes (10,10) (10,5,2)


Arrays must match in all dimensions in order to be compatible:

np.ones([3, 3]) * np.ones([3, 3])  # Note elementwise multiply, *not* matrix multiply.

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

m3 = np.arange(9).reshape([3, 3])
m3

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

m4 = np.arange(9, 18).reshape([3, 3])
m4

array([[ 9, 10, 11],
[12, 13, 14],
[15, 16, 17]])

m3 * m4  # Note elementwise multiply, *not* matrix multiply.

array([[  0,  10,  22],
[ 36,  52,  70],
[ 90, 112, 136]])


Except, that if one array has any Dimension 1, then the data is REPEATED to match the other.

col = np.arange(10).reshape([10, 1])
col

array([[0],
[1],
[2],
[3],
[4],
[5],
[6],
[7],
[8],
[9]])

row = col.transpose()
row

array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])

col.shape  # "Column Vector"

(10, 1)

row.shape  # "Row Vector"

(1, 10)

row + col

array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
[ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10],
[ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
[ 3,  4,  5,  6,  7,  8,  9, 10, 11, 12],
[ 4,  5,  6,  7,  8,  9, 10, 11, 12, 13],
[ 5,  6,  7,  8,  9, 10, 11, 12, 13, 14],
[ 6,  7,  8,  9, 10, 11, 12, 13, 14, 15],
[ 7,  8,  9, 10, 11, 12, 13, 14, 15, 16],
[ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17],
[ 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]])

10 * row + col

array([[ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90],
[ 1, 11, 21, 31, 41, 51, 61, 71, 81, 91],
[ 2, 12, 22, 32, 42, 52, 62, 72, 82, 92],
[ 3, 13, 23, 33, 43, 53, 63, 73, 83, 93],
[ 4, 14, 24, 34, 44, 54, 64, 74, 84, 94],
[ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95],
[ 6, 16, 26, 36, 46, 56, 66, 76, 86, 96],
[ 7, 17, 27, 37, 47, 57, 67, 77, 87, 97],
[ 8, 18, 28, 38, 48, 58, 68, 78, 88, 98],
[ 9, 19, 29, 39, 49, 59, 69, 79, 89, 99]])


This works for arrays with more than one unit dimension.

## 3.5.3 Another example#

x = np.array([1, 2]).reshape(1, 2)
x

array([[1, 2]])

y = np.array([3, 4, 5]).reshape(3, 1)
y

array([[3],
[4],
[5]])

result = x + y
result.shape

(3, 2)

result

array([[4, 5],
[5, 6],
[6, 7]])


What numpy is doing:

## 3.5.4 Newaxis#

Broadcasting is very powerful, and numpy allows indexing with np.newaxis to temporarily create new one-long dimensions on the fly.

import numpy as np

x = np.arange(10).reshape(2, 5)
y = np.arange(8).reshape(2, 2, 2)

x

array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])

y

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

[[4, 5],
[6, 7]]])

x_dash = x[:, :, np.newaxis, np.newaxis]
x_dash.shape

(2, 5, 1, 1)

y_dash = y[:, np.newaxis, :, :]
y_dash.shape

(2, 1, 2, 2)

y_dash

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

[[[4, 5],
[6, 7]]]])

res = x_dash * y_dash

res.shape

(2, 5, 2, 2)

np.sum(res)

830


Note that newaxis works because a $$3 \times 1 \times 3$$ array and a $$3 \times 3$$ array contain the same data, differently shaped:

threebythree = np.arange(9).reshape(3, 3)
threebythree

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

threebythree[:, np.newaxis, :]

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

[[3, 4, 5]],

[[6, 7, 8]]])


## 3.5.5 Dot Products using broadcasting#

NumPy multiply is element-by-element, not a dot-product:

a = np.arange(9).reshape(3, 3)
a

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

b = np.arange(3, 12).reshape(3, 3)
b

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

a * b

array([[ 0,  4, 10],
[18, 28, 40],
[54, 70, 88]])


We can use what we’ve learned about the algebra of broadcasting and newaxis to get a dot-product, (matrix inner product).

First we add new axes to $$A$$ and $$B$$:

a[:, :, np.newaxis].shape

(3, 3, 1)

b[np.newaxis, :, :].shape

(1, 3, 3)


Now we use broadcasting to generate $$A_{ij}B_{jk}$$ as a 3-d matrix:

a[:, :, np.newaxis] * b[np.newaxis, :, :]

array([[[ 0,  0,  0],
[ 6,  7,  8],
[18, 20, 22]],

[[ 9, 12, 15],
[24, 28, 32],
[45, 50, 55]],

[[18, 24, 30],
[42, 49, 56],
[72, 80, 88]]])


Then we sum over the middle, $$j$$ axis, [which is the 1-axis of three axes numbered (0,1,2)] of this 3-d matrix. Thus we generate $$\Sigma_j A_{ij}B_{jk}$$.

(a[:, :, np.newaxis] * b[np.newaxis, :, :]).sum(1)

array([[ 24,  27,  30],
[ 78,  90, 102],
[132, 153, 174]])


Or if you prefer:

(a.reshape(3, 3, 1) * b.reshape(1, 3, 3)).sum(1)

array([[ 24,  27,  30],
[ 78,  90, 102],
[132, 153, 174]])


We can see that the broadcasting concept gives us a powerful and efficient way to express many linear algebra operations computationally.

## 3.5.6 Dot Products using numpy functions#

However, as the dot-product is a common operation, numpy has a built in function:

np.dot(a, b)

array([[ 24,  27,  30],
[ 78,  90, 102],
[132, 153, 174]])


This can also be written as:

a.dot(b)

array([[ 24,  27,  30],
[ 78,  90, 102],
[132, 153, 174]])


If you are using Python 3.5 or later, a dedicated matrix multiplication operator has been added, allowing you to do the following:

a @ b

array([[ 24,  27,  30],
[ 78,  90, 102],
[132, 153, 174]])


## 3.5.7 Record Arrays#

These are a special array structure designed to match the CSV “Record and Field” model. It’s a very different structure from the normal NumPy array, and different fields can contain different datatypes. We saw this when we looked at CSV files:

x = np.arange(50).reshape([10, 5])

record_x = x.view(
dtype={"names": ["col1", "col2", "another", "more", "last"], "formats": [int] * 5}
)

record_x

array([[( 0,  1,  2,  3,  4)],
[( 5,  6,  7,  8,  9)],
[(10, 11, 12, 13, 14)],
[(15, 16, 17, 18, 19)],
[(20, 21, 22, 23, 24)],
[(25, 26, 27, 28, 29)],
[(30, 31, 32, 33, 34)],
[(35, 36, 37, 38, 39)],
[(40, 41, 42, 43, 44)],
[(45, 46, 47, 48, 49)]],
dtype=[('col1', '<i8'), ('col2', '<i8'), ('another', '<i8'), ('more', '<i8'), ('last', '<i8')])


Record arrays can be addressed with field names like they were a dictionary:

record_x["col1"]

array([[ 0],
[ 5],
[10],
[15],
[20],
[25],
[30],
[35],
[40],
[45]])


Indeed we can use these methods when parsing CSV files instead of using Pandas.

## 3.5.8 Logical arrays, masking, and selection#

Numpy defines operators like == and < to apply to arrays element by element:

x = np.zeros([3, 4])
x

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

y = np.arange(-1, 2)[:, np.newaxis] * np.arange(-2, 2)[np.newaxis, :]
y

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

y_is_one = y == 1
y_is_one

array([[False,  True, False, False],
[False, False, False, False],
[False, False, False,  True]])

aresame = x == y
aresame

array([[False, False,  True, False],
[ True,  True,  True,  True],
[False, False,  True, False]])


A logical array can be used to select elements from an array:

y[np.logical_not(aresame)]

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


Although when printed, this comes out as a flat list, if assigned to, the selected elements of the array are changed!

y[aresame] = 5

y

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


## 3.5.9 Numpy memory#

NumPy manages memory differently from lists. Changing an element in a copy of a list does not change the original list.

x = list(range(5))
y = x[:]

y[2] = 0
x

[0, 1, 2, 3, 4]


But in NumPy, changing the copy does change the original array!

x = np.arange(5)
y = x[:]

y[2] = 0
x

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


We must use np.copy to force separate memory. Otherwise NumPy tries its hardest to make slices be views on data.