Efficient repeated evaluation of splines at fixed points.
Project description
cached_interpolate
Efficient evaluation of interpolants at fixed points.
When performing a Bayesian analysis with a stochastic sampler however, one sometimes has to evaluate many interpolants (representing the parameterized model) at the same set of points (the data) with the same knot points.
Evaluating interpolants typically requires two stages:
- finding the closest knot of the interpolant to the new point and the distance from that knot.
- evaluating the interpolant at that point.
When we have identical knots and evaluation points but different functions being approximated the first of these stages is done many times unnecessarily. This can be made more efficient by caching the locations of the evaluation points leaving just the evaluation of the interpolation coefficients to be done at each iteration.
A further advantage of this, is that it allows trivially parallising the interpolation using cupy
.
This package implements this caching for nearest neighbour, linear, and cubic interpolation.
Installation
Currently this is only installable by downloading the source and installing locally.
$ git clone git@github.com:ColmTalbot/cached_interpolation.git
$ cd cached_interpolation
$ pip install .
Demonstration
We can compare the interpolation to the scipy.interpolate.CubicSpline
implementation.
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from cached_interpolate import CachingInterpolant
x_nodes = np.linspace(0, 1, 10)
y_nodes = np.random.uniform(-1, 1, 10)
evaluation_points = np.sort(np.random.uniform(0, 1, 10000))
interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, kind="cubic")
sp_interpolant = CubicSpline(x=x_nodes, y=y_nodes, bc_type="natural")
figure, axes = plt.subplots(nrows=2, gridspec_kw={'height_ratios': [3, 1]}, sharex=True)
axes[0].plot(evaluation_points, interpolant(evaluation_points))
axes[0].plot(evaluation_points, sp_interpolant(evaluation_points))
axes[0].scatter(x_nodes, y_nodes, color="k")
axes[1].plot(evaluation_points, interpolant(evaluation_points) - sp_interpolant(evaluation_points))
axes[0].set_ylabel("$y$")
axes[1].set_xlabel("$x$")
axes[1].set_ylabel("$\\Delta y$")
axes[1].set_xlim(0, 1)
plt.tight_layout()
plt.show()
plt.close(figure)
I note here that we use the "natural" boundary condition.
This means that first and second derivatives of the spline vanish at the endpoints.
This is different from the default "not-a-knot" boundary condition used in scipy
.
We can now evaluate this interpolant in a loop to demonstrate the performance
of the caching.
At every iteration we change the function being interpolated by randomly setting the
y
argument.
Before the loop, we add an initial call to set up the cache as this will be much
slower than the following iterations.
import time
import numpy as np
from scipy.interpolate import CubicSpline
from cached_interpolate import CachingInterpolant
x_nodes = np.linspace(0, 1, 10)
y_nodes = np.random.uniform(-1, 1, 10)
evaluation_points = np.random.uniform(0, 1, 10000)
interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, kind="cubic")
sp_interpolant = CubicSpline(x=x_nodes, y=y_nodes, bc_type="natural")
interpolant(x=evaluation_points, y=y_nodes)
n_iterations = 1000
start = time.time()
for _ in range(n_iterations):
y_nodes = np.random.uniform(-1, 1, 10)
interpolant(x=evaluation_points, y=y_nodes)
stop = time.time()
print(f"Cached time = {(stop - start):.3f}s for {n_iterations} iterations.")
start = time.time()
for _ in range(n_iterations):
y_nodes = np.random.uniform(-1, 1, 10)
CubicSpline(x=x_nodes, y=y_nodes, bc_type="natural")(evaluation_points)
stop = time.time()
print(f"Scipy time = {(stop - start):.3f}s for {n_iterations} iterations.")
Cached time = 0.187s for 1000 iterations.
Scipy time = 0.450s for 1000 iterations.
We've gained a factor of ~2.5 over the scipy
version without caching.
If this is the dominant cost in a simulation that takes a week to run, this is a big improvement.
If we need to evaluate for a new set of points, we have to tell the interpolant to reset the cache. There are two ways to do this:
- create a new interpolant, this will require reevaluating the interplation coefficients.
- disable the evaluation point caching.
import numpy as np
from scipy.interpolant import CubicSpline
from cached_interpolate import CachingInterpolant
x_nodes = np.linspace(0, 1, 10)
y_nodes = np.random.uniform(-1, 1, 10)
evaluation_points = np.random.uniform(0, 1, 10000)
interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, kind="cubic")
interpolated_values = interpolant(evaluation_points)
new_evaluation_points = np.random.uniform(0, 1, 1000)
interpolant(x=new_evaluation_points, use_cache=False)
Using the code in this way is much slower than scipy
and so not practically very useful.
If you have access to an Nvidia
GPU and are evaluating the spline at ~ O(10^5) or more points you may want to switch
to the cupy
backend.
This uses cupy
just for the evaluation stage, not for computing the interpolation coefficients.
import cupy as cp
import numpy as np
from cached_interpolate import CachingInterpolant
x_nodes = np.linspace(0, 1, 10)
y_nodes = np.random.uniform(-1, 1, 10)
evaluation_points = np.random.uniform(0, 1, 10000)
evaluation_points = cp.asarray(evaluation_points)
interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, backend=cp)
interpolated_values = interpolant(evaluation_points)
We can now repeat our timing test. To make the use of a GPU more realistic we'll increase the number of evaluation points.
import time
import cupy as cp
import numpy as np
from scipy.interpolate import CubicSpline
from cached_interpolate import CachingInterpolant
x_nodes = np.linspace(0, 1, 10)
y_nodes = np.random.uniform(-1, 1, 10)
evaluation_points = np.random.uniform(0, 1, 100000)
cp_evaluation_points = cp.asarray(evaluation_points)
interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, kind="cubic")
cp_interpolant = CachingInterpolant(x=x_nodes, y=y_nodes, kind="cubic", backend=cp)
sp_interpolant = CubicSpline(x=x_nodes, y=y_nodes, bc_type="natural")
interpolant(x=evaluation_points)
cp_interpolant(x=cp_evaluation_points)
n_iterations = 1000
start = time.time()
for _ in range(n_iterations):
y_nodes = np.random.uniform(-1, 1, 10)
interpolant(x=evaluation_points, y=np.random.uniform(-1, 1, 10))
stop = time.time()
print(f"CPU cached time = {(stop - start):.3f}s for {n_iterations} iterations.")
start = time.time()
for _ in range(n_iterations):
y_nodes = np.random.uniform(-1, 1, 10)
CubicSpline(x=x_nodes, y=y_nodes, bc_type="natural")(evaluation_points)
stop = time.time()
print(f"Scipy time = {(stop - start):.3f}s for {n_iterations} iterations.")
start = time.time()
for _ in range(n_iterations):
y_nodes = np.random.uniform(-1, 1, 10)
cp_interpolant(x=cp_evaluation_points, y=np.random.uniform(-1, 1, 10))
stop = time.time()
print(f"GPU cached time = {(stop - start):.3f}s for {n_iterations} iterations.")
CPU cached time = 3.367s for 1000 iterations.
Scipy time = 4.213s for 1000 iterations.
GPU cached time = 0.212s for 1000 iterations.
While there are likely more optimizations that can be made and improved flexibility in the implementation, we can see that the GPU version is well over an order of magnitude faster than either of the CPU versions.
If you have any comments/questions feel free to contact me through the issue tracker or a pull request.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distributions
Hashes for cached_interpolate-0.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | e2f70edf7512b18756c38747346dd2e0aa62bfe39f81e4f9432efd6e3c680ff0 |
|
MD5 | 5297b4062e04d5046b43a366b768a53e |
|
BLAKE2b-256 | bc02fe2ff57c7667b71afc24733522b6efe8af3846eeaf50aeb1571a832fb109 |
Hashes for cached_interpolate-0.1.0-cp310-cp310-macosx_10_15_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | f83fa6017dc0f4cd174b5bc05dbcb65e940d6961da911056ced06734abe42f78 |
|
MD5 | 41a2234a4446660abc4d612f6672aa05 |
|
BLAKE2b-256 | 030c2ad5e27fa5a1862a0ee13f6161f452414213f55be5331ca6f28934a27b22 |
Hashes for cached_interpolate-0.1.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 462f8a3a4a7c929f7e6add135258ec2c34249abaf41b36d471952dd591291e4b |
|
MD5 | 6eadd7b6f53ef6783206117510c4b532 |
|
BLAKE2b-256 | 1df620c06f0ecc7b09d2355328967cc4300630afdcb47aef309db7dbc3fabb28 |
Hashes for cached_interpolate-0.1.0-cp39-cp39-macosx_10_15_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | b806a98668acc895bd358dcf3a175749599d8bb811e921a99d915b7c18fd550a |
|
MD5 | 45add81f01fee63c2423bcb8d4abbf99 |
|
BLAKE2b-256 | 904035bd8ff19b9538cfcca0a33b1d2982d0b6695addb6ebf411269ec9b8d833 |
Hashes for cached_interpolate-0.1.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 0d83bd952e5bdf454b1bda5a63b7337fd41c8948c3785232b020aac5041be373 |
|
MD5 | f239b5e88993e5b36c5e936659cee53c |
|
BLAKE2b-256 | 56cabe40e92c88874e963b41e9d840d28a8ad226bc1e1c62fb89574449ed363b |
Hashes for cached_interpolate-0.1.0-cp38-cp38-macosx_10_14_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | ca91bf39b981e659d0c80c0513bbeb2e5265e3acdfc2010b573f551202b856f0 |
|
MD5 | 0b547f3db2906067f0326a79bb38e035 |
|
BLAKE2b-256 | fcaf171706d5616858e751319372f68abca4453b851d21c492340645a701a78a |