Indexing

Setting up the data

Let's create the structures that will be used later in this notebook


In [1]:
import numpy as np

In [2]:
np.random.seed(42)  # Setting the random seed

In [3]:
# a vector: the argument to the array function is a Python list
v = np.random.rand(10)
v


Out[3]:
array([0.37454012, 0.95071431, 0.73199394, 0.59865848, 0.15601864,
       0.15599452, 0.05808361, 0.86617615, 0.60111501, 0.70807258])

In [4]:
# a matrix: the argument to the array function is a nested Python list
M = np.random.rand(10, 2)
M


Out[4]:
array([[0.02058449, 0.96990985],
       [0.83244264, 0.21233911],
       [0.18182497, 0.18340451],
       [0.30424224, 0.52475643],
       [0.43194502, 0.29122914],
       [0.61185289, 0.13949386],
       [0.29214465, 0.36636184],
       [0.45606998, 0.78517596],
       [0.19967378, 0.51423444],
       [0.59241457, 0.04645041]])

We can index elements in an array using the square bracket and indices:


In [5]:
# v is a vector, and has only one dimension, taking one index
v[0]


Out[5]:
0.3745401188473625

In [6]:
# M is a matrix, or a 2 dimensional array, taking two indices 
M[1,1]


Out[6]:
0.21233911067827616

If we omit an index of a multidimensional array it returns the whole row (or, in general, a N-1 dimensional array)


In [7]:
M[1]


Out[7]:
array([0.83244264, 0.21233911])

The same thing can be achieved with using : instead of an index:


In [8]:
M[1,:] # row 1


Out[8]:
array([0.83244264, 0.21233911])

In [9]:
M[:,1] # column 1


Out[9]:
array([0.96990985, 0.21233911, 0.18340451, 0.52475643, 0.29122914,
       0.13949386, 0.36636184, 0.78517596, 0.51423444, 0.04645041])

We can assign new values to elements in an array using indexing:


In [10]:
M[0,0] = 1

In [11]:
M


Out[11]:
array([[1.        , 0.96990985],
       [0.83244264, 0.21233911],
       [0.18182497, 0.18340451],
       [0.30424224, 0.52475643],
       [0.43194502, 0.29122914],
       [0.61185289, 0.13949386],
       [0.29214465, 0.36636184],
       [0.45606998, 0.78517596],
       [0.19967378, 0.51423444],
       [0.59241457, 0.04645041]])

In [12]:
# also works for rows and columns
M[1,:] = 0
M[:,1] = -1

In [13]:
M


Out[13]:
array([[ 1.        , -1.        ],
       [ 0.        , -1.        ],
       [ 0.18182497, -1.        ],
       [ 0.30424224, -1.        ],
       [ 0.43194502, -1.        ],
       [ 0.61185289, -1.        ],
       [ 0.29214465, -1.        ],
       [ 0.45606998, -1.        ],
       [ 0.19967378, -1.        ],
       [ 0.59241457, -1.        ]])

Index slicing

Index slicing is the technical name for the syntax M[lower:upper:step] to extract part of an array:


In [14]:
a = np.array([1,2,3,4,5])
a


Out[14]:
array([1, 2, 3, 4, 5])

In [15]:
a[1:3]


Out[15]:
array([2, 3])

Array slices are mutable: if they are assigned a new value the original array from which the slice was extracted is modified:


In [16]:
a[1:3] = [-2,-3]

a


Out[16]:
array([ 1, -2, -3,  4,  5])
  • We can omit any of the three parameters in M[lower:upper:step]:

In [17]:
a[::] # lower, upper, step all take the default values


Out[17]:
array([ 1, -2, -3,  4,  5])

In [18]:
a[::2] # step is 2, lower and upper defaults to the beginning and end of the array


Out[18]:
array([ 1, -3,  5])

In [19]:
a[:3] # first three elements


Out[19]:
array([ 1, -2, -3])

In [20]:
a[3:] # elements from index 3


Out[20]:
array([4, 5])
  • Negative indices counts from the end of the array (positive index from the begining):

In [21]:
a = np.array([1,2,3,4,5])

In [22]:
a[-1] # the last element in the array


Out[22]:
5

In [23]:
a[-3:] # the last three elements


Out[23]:
array([3, 4, 5])
  • Index slicing works exactly the same way for multidimensional arrays:

In [24]:
A = np.array([[n+m*10 for n in range(5)] 
              for m in range(5)])
A


Out[24]:
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44]])

In [25]:
# a block from the original array
A[1:4, 1:4]


Out[25]:
array([[11, 12, 13],
       [21, 22, 23],
       [31, 32, 33]])

In [26]:
# strides
A[::2, ::2]


Out[26]:
array([[ 0,  2,  4],
       [20, 22, 24],
       [40, 42, 44]])

Indexing and Array Memory Management

Numpy arrays support two different way of storing data into memory, namely

  • F-Contiguous
    • i.e. column-wise storage, Fortran-like
  • C-Contiguous
    • i.e. row-wise storage, C-like

The storage strategy is controlled by the parameter order of np.array


In [27]:
import numpy as np
FC = np.array([[1, 2, 3],    [4, 5, 6], 
               [7, 8, 9],    [10, 11, 12]], order='F')

In [28]:
CC = np.array([[1, 2, 3],    [4, 5, 6], 
               [7, 8, 9],    [10, 11, 12]], order='C')
  • Note: no changes in meaning for indexing operations

In [29]:
FC[0, 1]


Out[29]:
2

In [30]:
CC[0, 1]


Out[30]:
2

In [31]:
FC.shape


Out[31]:
(4, 3)

In [32]:
CC.shape


Out[32]:
(4, 3)

Fancy Indexing

Fancy indexing is the name for when an array or list is used in-place of an index:


In [33]:
row_indices = [1, 2, 3]
A[row_indices]


Out[33]:
array([[10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34]])

In [34]:
col_indices = [1, 2, -1] # remember, index -1 means the last element
A[row_indices, col_indices]


Out[34]:
array([11, 22, 34])
  • We can also index masks:

    • If the index mask is an Numpy array of with data type bool, then an element is selected (True) or not (False) depending on the value of the index mask at the position each element:

In [35]:
b = np.array([n for n in range(5)])
b


Out[35]:
array([0, 1, 2, 3, 4])

In [36]:
row_mask = np.array([True, False, True, False, False])
b[row_mask]


Out[36]:
array([0, 2])
  • Alternatively:

In [37]:
# same thing
row_mask = np.array([1,0,1,0,0], dtype=bool)
b[row_mask]


Out[37]:
array([0, 2])

This feature is very useful to conditionally select elements from an array, using for example comparison operators:


In [38]:
x = np.arange(0, 10, 0.5)
x


Out[38]:
array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. ,
       6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5])

In [39]:
mask = (5 < x)

mask


Out[39]:
array([False, False, False, False, False, False, False, False, False,
       False, False,  True,  True,  True,  True,  True,  True,  True,
        True,  True])

In [40]:
x[mask]


Out[40]:
array([5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5])

Alternatively, we can use the condition (mask) array directly within brackets to index the array


In [41]:
x[(5 < x)]


Out[41]:
array([5.5, 6. , 6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5])

Exercises on Indexing

Index slicing is the technical name for the syntax M[lower:upper:step] to extract part of an array

Ex 2.1

Generate a three-dimensional array of any size containing random numbers taken from an uniform distribution (guess the numpy function in np.random). Then print out separately the first entry along the three axis (i.e. x, y, z)

  • hint: Slicing with numpy arrays works quite like Python lists

In [ ]:

Ex 2.2

Create a vector and print out elements in reverse order

Hint: Use slicing for this exercise


In [ ]:

Ex 2.3

Generate a $7 \times 7$ matrix and replace all the elements in odd rows and even columns with 1.

Hint: Use slicing to solve this exercise!

Note: Take a look at the original matrix, then.


In [ ]:

Use fancy indexing to get all the elements of the previous matrix that are equals to 1


In [ ]:

Ex 2.4

Generate a 10 x 10 matrix of numbers A. Then, generate a numpy array of integers in range 1-9. Pick 5 random values (with no repetition) from this array and use these values to extract rows from the original matrix A.


In [ ]:

Ex 2.5

Repeat the previous exercise but this time extract columns from A


In [ ]:

Ex 2.6

Generate an array of numbers from 0 to 20 with step 0.5. Extract all the values greater than a randomly generated number in the same range.

  • hint: Try to write the condition as an expression and save it to a variable. Then, use this variable in square brackets to index.... this is when the magic happens!

In [ ]:


Stacking and repeating arrays

Using function repeat, tile, vstack, hstack, and concatenate we can create larger vectors and matrices from smaller ones:

np.tile and np.repeat


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

np.repeat


In [43]:
# repeat each element 3 times
np.repeat(a, 3)


Out[43]:
array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4])

np.tile


In [44]:
# tile the matrix 3 times 
np.tile(a, 3)


Out[44]:
array([[1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4]])

np.concatenate


In [45]:
b = np.array([[5, 6]])

In [46]:
np.concatenate((a, b), axis=0)


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

In [47]:
np.concatenate((a, b.T), axis=1)


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

np.hstack and np.vstack


In [48]:
np.vstack((a,b))


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

In [49]:
np.hstack((a,b.T))


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

Iterating over array elements

Generally, we want to avoid iterating over the elements of arrays whenever we can (at all costs). The reason is that in a interpreted language like Python (or MATLAB), iterations are really slow compared to vectorized operations.

However, sometimes iterations are unavoidable. For such cases, the Python for loop is the most convenient way to iterate over an array:


In [57]:
v = np.array([1,2,3,4])

for element in v:
    print(element)


1
2
3
4

In [58]:
M = np.array([[1,2], [3,4]])

for row in M:
    print("row", row)
    
    for element in row:
        print(element)


row [1 2]
1
2
row [3 4]
3
4
  • When we need to iterate over each element of an array and modify its elements, it is convenient to use the enumerate function to obtain both the element and its index in the for loop:

In [59]:
for row_idx, row in enumerate(M):
    print("row_idx", row_idx, "row", row)
    
    for col_idx, element in enumerate(row):
        print("col_idx", col_idx, "element", element)
       
        # update the matrix M: square each element
        M[row_idx, col_idx] = element ** 2


row_idx 0 row [1 2]
col_idx 0 element 1
col_idx 1 element 2
row_idx 1 row [3 4]
col_idx 0 element 3
col_idx 1 element 4

In [60]:
# each element in M is now squared
M


Out[60]:
array([[ 1,  4],
       [ 9, 16]])

Vectorizing functions

As mentioned several times by now, to get good performance we should try to avoid looping over elements in our vectors and matrices, and instead use vectorized algorithms. The first step in converting a scalar algorithm to a vectorized algorithm is to make sure that the functions we write work with vector inputs.


In [61]:
def Theta(x):
    """
    Scalar implemenation of the Heaviside step function.
    """
    if x >= 0:
        return 1
    else:
        return 0

In [62]:
Theta(array([-3,-2,-1,0,1,2,3]))


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-62-2cb2062a7e18> in <module>
----> 1 Theta(array([-3,-2,-1,0,1,2,3]))

NameError: name 'array' is not defined

OK, that didn't work because we didn't write the Theta function so that it can handle with vector input...

To get a vectorized version of Theta we can use the Numpy function np.vectorize. In many cases it can automatically vectorize a function:

np.vectorize


In [63]:
Theta_vec = np.vectorize(Theta)

In [65]:
Theta_vec(np.array([-3,-2,-1,0,1,2,3]))


Out[65]:
array([0, 0, 0, 1, 1, 1, 1])

np.frompyfunc

Universal functions (Ufuncs) work on arrays, element-by-element, or on scalars.

Ufuncs accept a set of scalars as input, and produce a set of scalars as output.


In [66]:
Theta_ufunc = np.frompyfunc(Theta, 1, 1)
print("Result: ", Theta_ufunc(np.arange(4)))


Result:  [1 1 1 1]

Excercise

Avoiding Vectorize

  • Implement the function to accept vector input from the beginning
    • This requires "more effort" but might give better performance

Vectorisation and Broadcasting

Vectorizing code is the key to writing efficient numerical calculation with Python/Numpy.

That means that as much as possible of a program should be formulated in terms of matrix and vector operations, like matrix-matrix multiplication.

However, sometimes it may happen that we need to perform operations between ndarray objects of different size and shape so we might be tempted of performing this operation explicitly (i.e. explicit for loop).

Although we will get correct results (assumed we iterated arrays in the right way), this will produce very inefficient code.

For this reason, NumPy has a very important concept called Broadcasting.

Broadcasting is Numpy's terminology for performing mathematical operations between arrays with different shapes (assuming scalars being array with 0-dim).

Working example

Assume we have some data in a matrix D:


In [72]:
D = np.array([ [0.3, 2.5, 3.5],
               [2.9, 27.5, 0],
               [0.4, 1.3, 23.9],
               [14.4, 6, 2.3]])

And we also have another vector adj of values that contains some adjusting factors that we might want to apply to each sample (row) of data in D


In [73]:
adj = np.array([9, 4, 4])

Naive Solution


In [79]:
%%timeit

# Create a new array filled with zeros, of the same shape as macros.
result = np.zeros_like(D)

# Now multiply each row of macros by cal_per_macro. In Numpy, `*` is
# element-wise multiplication between two arrays.
for i in range(D.shape[0]):
    result[i, :] = D[i, :] * adj
    
result


12.5 µs ± 375 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

This is a reasonable approach when coding in a low-level programming language: allocate the output, loop over input performing some operation, write result into output.

In Numpy, however, this is fairly bad for performance because the looping is done in (slow) Python code instead of internally by Numpy in (fast) C code

(we've proven this already when we compared ndarray vs list)

Using np.tile

Idea: Leverage on np.tile function to replicate adj over the number of rows D has (i.e. D.shape[0]).


In [80]:
%%timeit 

adj_stretch = np.tile(adj, (D.shape[0], 1))

D * adj_stretch


8.78 µs ± 1.31 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [81]:
adj_stretch


Out[81]:
array([[9, 4, 4],
       [9, 4, 4],
       [9, 4, 4],
       [9, 4, 4]])

Nice, it's shorter too, and slightly faster! To appreciate even more performance gain, of our np.tile solution, we could try increasing the size of D to a bigger structure:


In [84]:
D_large = np.random.rand(10**6, 10)
adj_large = np.random.rand(10)

In [85]:
D_large.shape, adj_large.shape


Out[85]:
((1000000, 10), (10,))

In [86]:
%%timeit

# Create a new array filled with zeros, of the same shape as macros.
result_large = np.zeros_like(D_large)

# Now multiply each row of macros by cal_per_macro. In Numpy, `*` is
# element-wise multiplication between two arrays.
for i in range(D_large.shape[0]):
    result_large[i, :] = D_large[i, :] * adj_large


1.48 s ± 126 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [87]:
%%timeit 

adj_large_stretch = np.tile(adj_large, (D_large.shape[0], 1))

D_large * adj_large_stretch


48.6 ms ± 1.14 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

The loop-in-Python method takes ~1.5 seconds, the stretching method takes ~48 milliseconds, a ~75x speedup.

And now, finally, comes the interesting part.

You see, the operation we just did - stretching one array so that its shape matches that of another and then applying some element-wise operation between them - is actually pretty common.

This often happens when we want to take a lower-dimensional array and use it to perform a computation along some axis of a higher-dimensional array.

In fact, when taken to the extreme this is exactly what happens when we perform an operation between an array and a scalar - the scalar is stretched across the whole array so that the element-wise operation gets the same scalar value for each element it computes.

Numpy generalizes this concept into broadcasting - a set of rules that permit element-wise computations between arrays of different shapes, as long as some constraints apply.

Incidentally, this lets Numpy achieve two separate goals - usefulness as well as more consistent and general semantics.

Binary operators like * are element-wise, but what happens when we apply them between arrays of different shapes? Should it work or should it be rejected? If it works, how should it work?

Broadcasting defines the semantics of these operations.


In [89]:
%%timeit
## Back to our example
D * adj  # Broadcasting


1.79 µs ± 103 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

How Broadcasting works

Element-wise operations on arrays are only valid when the arrays' shapes are either equal or compatible.

The equal shapes case is trivial - this is the stretched array from the example above. What does compatible mean, though?

To determine if two shapes are compatible, Numpy compares their dimensions, starting with the trailing ones and working its way backwards.

For example, for the shape (4, 3, 2) the trailing dimension is 2, and working from 2 "backwards" produces: 2, then 3, then 4.

If two dimensions are equal, or if one of them equals 1, the comparison continues. Otherwise, you'll see a ValueError raised (saying something like "operands could not be broadcast together with shapes ...").

When one of the shapes runs out of dimensions (because it has less dimensions than the other shape), Numpy will use 1 in the comparison process until the other shape's dimensions run out as well.

Once Numpy determines that two shapes are compatible, the shape of the result is simply the maximum of the two shapes' sizes in each dimension.

Scalar-array operations

We can use the usual arithmetic operators to multiply, add, subtract, and divide arrays with scalar numbers.


In [90]:
v1 = np.arange(0, 5)

In [91]:
v1 * 2


Out[91]:
array([0, 2, 4, 6, 8])

In [92]:
v1 + 2


Out[92]:
array([2, 3, 4, 5, 6])

In [93]:
A = np.array([[n+m*10 for n in range(5)] for m in range(5)])

In [94]:
print('A * 2: ', '\n', A * 2)
print('A + 2: ', '\n', A + 2)


A * 2:  
 [[ 0  2  4  6  8]
 [20 22 24 26 28]
 [40 42 44 46 48]
 [60 62 64 66 68]
 [80 82 84 86 88]]
A + 2:  
 [[ 2  3  4  5  6]
 [12 13 14 15 16]
 [22 23 24 25 26]
 [32 33 34 35 36]
 [42 43 44 45 46]]

Element-wise array-array operations

When we add, subtract, multiply and divide arrays with each other, the default behaviour is element-wise operations:


In [95]:
A * A # element-wise multiplication


Out[95]:
array([[   0,    1,    4,    9,   16],
       [ 100,  121,  144,  169,  196],
       [ 400,  441,  484,  529,  576],
       [ 900,  961, 1024, 1089, 1156],
       [1600, 1681, 1764, 1849, 1936]])

In [96]:
v1 * v1


Out[96]:
array([ 0,  1,  4,  9, 16])

If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:


In [97]:
A.shape, v1.shape


Out[97]:
((5, 5), (5,))

In [99]:
A * v1  #Broadcasting


Out[99]:
array([[  0,   1,   4,   9,  16],
       [  0,  11,  24,  39,  56],
       [  0,  21,  44,  69,  96],
       [  0,  31,  64,  99, 136],
       [  0,  41,  84, 129, 176]])

Matrix algebra

What about matrix mutiplication?

There are two ways.

We can either use the np.dot function, which applies a matrix-matrix, matrix-vector, or inner vector multiplication to its two arguments:


In [100]:
np.dot(A, A)


Out[100]:
array([[ 300,  310,  320,  330,  340],
       [1300, 1360, 1420, 1480, 1540],
       [2300, 2410, 2520, 2630, 2740],
       [3300, 3460, 3620, 3780, 3940],
       [4300, 4510, 4720, 4930, 5140]])

In [101]:
np.dot(A, v1)


Out[101]:
array([ 30, 130, 230, 330, 430])

In [102]:
np.dot(v1, v1)


Out[102]:
30

A new dedicated Infix operator for Matrix Multiplication

Python 3.5 (released on the 24th Aug.) introduces a new binary operator to be used for matrix multiplication, called @ . (Mnemonic: @ is * for mATrices.)

Some useful hints:


In [103]:
A @ v1


Out[103]:
array([ 30, 130, 230, 330, 430])