I implemented this comparison between numpy and scipy for doing the same function interpolation. The results show how numpy crushes scipy.
Python version: 3.11.7
NumPy version: 2.1.3
SciPy version: 1.15.2
Custom NumPy interpolation matches SciPy BPoly for 10 000 points.
SciPy coeff time: 0.218046 s
SciPy eval time : 0.000725 s
Custom coeff time: 0.061066 s
Custom eval time : 0.000550 s
edit: likely, I was too pessimistic below when giving the 4x to 10x, after latest streamlining of code, seems very significant on average.
Varying on system, I get somewhere between 4x to 10x speedup with numpy. This is not the first time I encounter this. So in general, I wonder:
why this huge difference in perf? should we do well in viewing scipy as just a reference implementation, and go to other pathways for peformance (numpy, numba)?
# BPoly.from_derivatives example with 10 000 sinusoidal points and timing comparison
"""
Creates a sinusoidal dataset of 10 000 points over [0, 2π].
Interpolates using SciPy's BPoly.from_derivatives and a custom pure NumPy quintic Hermite.
Verifies exact match, compares timing for coefficient computation and evaluation, and visualizes both.
"""
import sys
import numpy as np
import scipy
import time
from scipy.interpolate import BPoly
# Environment versions
print(f"Python version: {sys.version.split()[0]}")
print(f"NumPy version: {np.__version__}")
print(f"SciPy version: {scipy.__version__}")
# Generate 10 000 sample points over one period
n = 10_000
x = np.linspace(0.0, 2*np.pi, n)
# Analytical sinusoidal values and derivatives
y = np.sin(x) # y(x)
v = np.cos(x) # y'(x)
a = -np.sin(x) # y''(x)
# === SciPy implementation with timing ===
y_and_derivatives = np.column_stack((y, v, a))
t0 = time.perf_counter()
bp = BPoly.from_derivatives(x, y_and_derivatives)
t1 = time.perf_counter()
scipy_coeff_time = t1 - t0
# Evaluation timing
t0 = time.perf_counter()
y_scipy = bp(x)
t1 = time.perf_counter()
scipy_eval_time = t1 - t0
# === Pure NumPy implementation ===
def compute_quintic_coeffs(x, y, v, a):
"""
Compute quintic Hermite coefficients on each interval [x[i], x[i+1]].
Returns coeffs of shape (6, n-1).
"""
m = len(x) - 1
coeffs = np.zeros((6, m))
for i in range(m):
h = x[i+1] - x[i]
A = np.array([
[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 2, 0, 0, 0],
[1, h, h**2, h**3, h**4, h**5],
[0, 1, 2*h, 3*h**2, 4*h**3, 5*h**4],
[0, 0, 2, 6*h, 12*h**2, 20*h**3],
])
b = np.array([y[i], v[i], a[i], y[i+1], v[i+1], a[i+1]])
coeffs[:, i] = np.linalg.solve(A, b)
return coeffs
def interp_quintic(x, coeffs, xx):
"""
Evaluate quintic Hermite using precomputed coeffs.
x: breakpoints (n,), coeffs: (6, n-1), xx: query points.
"""
idx = np.searchsorted(x, xx) - 1
idx = np.clip(idx, 0, len(x) - 2)
dx = xx - x[idx]
yy = np.zeros_like(xx)
for j in range(6):
yy += coeffs[j, idx] * dx**j
return yy
# Coefficient estimation timing for custom
t0 = time.perf_counter()
custom_coeffs = compute_quintic_coeffs(x, y, v, a)
t1 = time.perf_counter()
cust_coeff_time = t1 - t0
# Evaluation timing for custom
t0 = time.perf_counter()
y_custom = interp_quintic(x, custom_coeffs, x)
t1 = time.perf_counter()
cust_eval_time = t1 - t0
# Verify exact match
assert np.allclose(y_scipy, y_custom, atol=1e-12), "Custom interp deviates"
print("Custom NumPy interpolation matches SciPy BPoly for 10 000 points.")
# Print timing results
print(f"SciPy coeff time: {scipy_coeff_time:.6f} s")
print(f"SciPy eval time : {scipy_eval_time:.6f} s")
print(f"Custom coeff time: {cust_coeff_time:.6f} s")
print(f"Custom eval time : {cust_eval_time:.6f} s")
# Visualization
import matplotlib.pyplot as plt
plt.plot(x, y_scipy, label='SciPy BPoly')
plt.plot(x, y_custom, '--', label='Custom NumPy')
plt.legend()
plt.title('Sinusoidal interpolation: SciPy vs NumPy Quintic')
plt.xlabel('x')
plt.ylabel('y(x) = sin(x)')
plt.show()
BPoly.from_derivatives has not seen much work after an initial implementation (which, back in the day, replaced an order of magnitude slower interpolator). This is for various reasons, one of which it is mostly superseded by CubicHermiteSpline. We barely saw it used in the wild, hence had little incentive to work on it. If you are using it or would use it if it were faster, please don't hesitate sending a pull request. It will be gladly reviewed.