Longitudinal Benchmarks with asv

[mrocklin]

In conversation here it was suggested by @minrk and @ogrisel that we set up something like https://github.com/spacetelescope/asv to monitor performance of a few important computation patterns over time. This would help us to identify performance regressions as we change the scheduler.

[mrocklin]

I already have a bunch of benchmarks that I use when profiling performance. The goal of this issue should be first to set up sufficient infrastructure so that it is easy for others to add benchmarks in the future.

[dhirschfeld]

In case there's desire to have a real-world test case I'm posting below a very simple Monte-Carlo GBM simulation. I've got similar code running on a distributed cluster so have a vested-interest in making sure it runs efficiently! I'm hopeful that a real-world use-case such as this can inform future optimisation efforts (though I wonder if it's too simple to be of any use?)

multivariate_normal(covariance, ndates, nsims, nchunks)

def multivariate_normal(covariance, ndates, nsims, nchunks):
    """Returns correlated random normal variables according to the supplied ``covariance``
    
    Parameters
    ----------
    covariance : ndarray[float](ndims, ndims)
        The size of the covariance is determined by ``ndims`` - the number
        of processes being simulated
    
    ndates : int
        The size of the time axis for the simulation
    
    nsims : int
        The number of simulations for each time period
    
    nchunks : Tuple[int, int, int]
        The size of each chunk
    
    Returns
    -------
    rvs : ndarray[float](ndims, ndates, nsims)
        The correlated rvs st ``np.cov(rvs.reshape(ndims, -1)) ~= covariance``
    
    """
    ndims = covariance.shape[0]
    U, s, Vh = np.linalg.svd(covariance)
    rvs = da.random.normal(size=(ndims, ndates, nsims), chunks=(ndims, 1+ndates//nchunks, nsims))
    rvs -= rvs.mean(axis=-1)[..., None]
    rvs /= rvs.std(axis=-1)[..., None]
    rvs = da.dot(U @ np.diag(np.sqrt(s)), rvs.reshape(2, -1)).reshape(rvs.shape)
    return rvs

simulate_gbm(F, sigma, Te, rvs)

def simulate_gbm(F, sigma, Te, rvs):
    """Simulates a correlated, n-dimensional Geometric Brownian Motion process
    
    Parameters
    ----------
    F : float
        The current (forward) price
    
    sigma : float
        The annualised volatility
    
    Te : ndarray[float](ndates,)
        The annualised time-to-expiry for each date (column of ``rvs``)
    
    rvs : ndarray[float](ndims, ndates, nsims)
        The correlated rvs
    
    Returns
    -------
    sims : ndarray[float](ndims, ndates, nsims)
        The transformed ``rvs`` st each path ``sims[i, :, k]`` represents
        a realisation of a GBM process with variance ``covariance[i, i]``

    """
    variance = Te * sigma**2
    mean = -0.5*variance
    dV = np.append(variance[0], np.diff(variance))
    sims = rvs
    sims *= np.sqrt(dV)[None, :, None]
    sims.cumsum(axis=1, out=sims)
    sims += mean[None, :, None]
    np.exp(sims, out=sims)
    sims *= F
    return sims

>>> F = 50.
>>> sigma = 0.2
>>> covariance = np.array([
...     [1,  0.7],
...     [0.7,  1],
... ])
>>> Te = np.linspace(0, 3, 1 + 365*3)
>>> rvs =  multivariate_normal(covariance, Te.size, nsims=1e5, nchunks=60)
>>> sims = simulate_gbm(F, sigma, Te, rvs)
>>> sims
dask.array<mul, shape=(2, 1096, 100000), dtype=float64, chunksize=(2, 19, 100000)>
>>> %%time
... fut = cluster.compute(sims, optimize_graph=True, pure=False)
... wait(fut)
Wall time: 11.8 s
>>> %%time
... sims_ = sims.compute(scheduler='synchronous')
Wall time: 46.8 s