[Python] SCIPY BPoly.from_derivatives compared with numpy

I implemented this comparison between numpy and scipy for doing the same function interpolation. The results show how numpy crushes scipy. 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 numpy as np
import time
from scipy.interpolate import BPoly

# 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, 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
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, v, a, 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()

Continue reading...