Case study 1: a bootstrap routine

Case study 1: a bootstrap routine#

Bootstrapping is a computational approach to resample with replacement \(B\) bootstrap datasets from an original sample. For example, let’s assume you want to construct the distribution of the statistic from a given (sometimes small sample).

From a coding perspective bootstrap routines are extremely simple. In its simplest form bootstrapping is essentially two nested for loops along with a function to compute the statistic of interest. This is quite an interesting case study for numpy i.e. how do we eliminate nested loops efficiently?

Imports#

from random import randint
import numpy as np

Dataset#

We will test our code using a simple sequence of 10 integer values. To avoid any overhead due to type conversion we will have one stored as a np.ndarray and the other as a list.

sample_np = np.arange(1, 11)
sample_lst = sample_np.tolist()
print(sample_np)

[ 1  2  3  4  5  6  7  8  9 10]

Benchmark: Standard python implementation#

Before we get into the numpy code here is some standard python code. We will benchmark its performance at generating 1000 boostrap datasets of the mean using the %timeit.

def bootstrap_python(data, boots):
    """
    Create bootstrap datasets that represent the distribution of the mean.
    Returns a list containing the bootstrap datasets 
    
    Parameters:
    -----------
    data: array-like
        vector of original data
    boots:int
        number of bootstrap samples to generate
        
    Returns:
    -------
    list
    """
    to_return = []
    
    # outer loop - create bootstrap datasets
    for b in range(boots):
        # resample with replacement and sum sample values
        
        # inner loop sample with replacement and calculate mean
        total = 0
        for _ in range(len(data)):
            
            # resample with replacement and increment running total. 
            total += data[randint(0, len(data)-1)]
        
        # mean of the bootstrap dataset
        to_return.append(total / len(data))

    return to_return

%timeit bootstrap_python(sample_lst, 1000)

6.24 ms ± 26.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

The exact performance might vary a bit on your machine, but you should get performance in milliseconds (ms). On my machine this was around 6.2 to 6.5ms on average.

Converting to `numpy`#

In theory, execution of loops in python is slow relative to a low level language like C. To speed up our code we need to eliminate them. This can sometimes be simple to do and sometimes require you to be much more creative!

Our first attempt will eliminate the inner loop (a list comprehension) in python

Attempt 1: the inner loop#

The function bootstrap_numpy below eliminates one of the python for loops - the inner list comprehension. Instead of iterating in python we make use numpy’s random number generator function integers. This has a size parameter that allows us to specify the number of loops that will occur ‘closer to the metal’.

Note that in bootstrap numpy we create and specify the size of our array to_return upfront.

def bootstrap_numpy(data, boots):
    """
    Create bootstrap datasets that represent the distribution of the mean.
    Makes better use of built in numpy methods for more efficient sampling
    Returns a numpy array containing the bootstrap datasets 
    
    Keyword arguments:
    data -- numpy array of systems to boostrap
    boots -- number of bootstrap 
    """
    
    # numpy random number generator
    rng = np.random.default_rng()
        
    # empty vector to return to caller. Its len=boots
    to_return = np.empty(boots)
              
    for b in range(boots):
        
        # we get numpy's randint function to do the iteration work for us.
        total = data[rng.integers(0, data.shape[0], size=data.shape[0])].sum() 

        to_return[b] = total / data.shape[0]

    return to_return

%timeit bootstrap_numpy(sample_np, 1000)

7.98 ms ± 89.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

An unexpected result! Again the exact results will vary, but on my machine the numpy implementation is actually marginally slower than the pure python implementation. Coming in a around 8ms per loop¬ There must be some overhead (hidden runtime cost) that we are missing. A hypothesis is that there may be some overhead in calling the rng.integers method from python. We do that boot times during the run.

We know that there is still one python loop in our code. Can we do better than the pure python code if we eliminate it and only call a numpy method once?

Attempt 2: eliminating the nested loops#

Bootstrapping is simply a sample with replacement. In most implementations we would use a nested loop structure as we want to create \(B\) bootstrap datasets (loop one) of size \(n\) (loop 2). But in general we are just taken \(B \times n\) samples with replacement.

One option to reimplement our algorithm in numpy is

create vector of length \(B \times n\).
Fill it will data sampled with replacement
Split it into \(B\) vectors of length \(n\)
Calculate a statistic of interest, e.g. the mean, for each of the \(B\) datasets.

We will implement steps 3 and 4 here by reshaping a 1D array into a \(B \times n\) matrix and take the mean of the rows (axis=1)

def bootstrap_better_numpy(data, boots):
    """
    Create bootstrap datasets that represent the distribution of the mean.
    Makes better use of built in numpy methods for more efficient sampling
    Returns a numpy array containing the bootstrap datasets 
    
    Parameters:
    ----------
    data: array-like
        original sample data
    boots: int
        number of bootstrap datasets to create.
        
    Returns:
    -------
    np.ndarray
    """
    # create random number generator
    rng = np.random.default_rng()
    
    # get all the bootstrap data in one go boots * len(data)
    boot_data = data[rng.integers(0, data.shape[0], size=data.shape[0]*boots)]
    
    # reshape array into (boots, len(data)) and calculate mean
    return boot_data.reshape(-1, len(data)).sum(axis=1) / len(data)

%timeit bootstrap_better_numpy(sample_np, 1000)

108 µs ± 538 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

That’s improvement! We are now seeing an order of magnitude speed up in the code. On my machine the average runtime of the function is less than the variation between runs of the pure python code.

Step by step breakdown#

Our new numpy implementation executes a lot quicker. But it is harder to understand how it works. Let’s break it down step by step.

Random number generator#

To do our sampling we create a random number generator. This has a function called integers.

# create random number generator
rng = np.random.default_rng()
rng.integers(0, 10, size=5)

array([1, 3, 9, 9, 8])

Sampling all data in one go#

The next line of code does all of the sampling with replacement that we need.

boot_data = data[rng.integers(0, data.shape[0], size=data.shape[0]*boots)]

We have set the size parameter to the length of our data vector multipled by the number of bootstraps needed. These figures are 10 and 1000 respectively. So we end up with a vector of length 10,000

data = sample_np
boots = 1000

boot_data = data[rng.integers(0, data.shape[0], size=data.shape[0]*boots)]
print(boot_data.shape)

(10000,)

Reshaping into a \(B \times n\) matrix#

The next line of code combines steps 3 and 4.

return boot_data.reshape(-1, len(data)).sum(axis=1) / len(data)

Its easier to understand if we look at the result of the reshaping first.

matrix = boot_data.reshape(-1, len(data))
print(matrix)

[[ 7  8  3 ...  3  5  3]
 [ 6  9  9 ...  9  8  2]
 [10  3  5 ...  3  9  2]
 ...
 [ 5  5  6 ...  5  7  7]
 [ 8  4  7 ...  4  6  7]
 [ 8  1  8 ...  3  8  4]]

In the reshaping we have specified that there needs to be len(data) columns while the number of rows is automatically inferred by number (-1). So the result is a matrix with the same dimensions as we would get using nested for loops (but with the performance bonus there there are no native python for loops!).

From here it just a case of summing over the columns .sum(axis=1) and dividing by the number of samples in each bootstrap dataset i.e. len(data).

Summing up.#

We saw a dramatic performance boost from numpy here. However, to get there is not always intuitive (and to be frank not as easy as this example). The result is also not as readable as the original code. The trade-off is that this may be the only way to make an algorithm feasible to execute in python.