Квинтический Эрмит на NumPy vs SciPy BPoly.from_derivatives: разбор и выбор API

Сравниваем квинтическую эрмитову интерполяцию на NumPy и SciPy BPoly.from_derivatives: бенчмарк, причины разрыва и рекомендация — CubicHermiteSpline.

Когда чистая реализация квинтической эрмитовой интерполяции на NumPy обгоняет SciPy BPoly.from_derivatives с большим отрывом, естественно задаться вопросом, откуда берётся разница и как выбрать подходящее API. Ниже — компактное руководство, которое воспроизводит сравнение, объясняет первопричину со слов мейнтейнеров и указывает прагматичный путь вперёд.

Воспроизводим сравнение

Следующий скрипт строит синусоиду из 10 000 точек с аналитическими первой и второй производными, создаёт интерполянт с помощью SciPy BPoly.from_derivatives, а затем делает то же самое вручную на NumPy — квинтический Эрмит через поинтервальные решения 6×6. Скрипт проверяет совпадение значений интерполянтов в узлах и замеряет время построения коэффициентов и вычисления. Имена переменных и функций выбраны для ясности, а вычислительный поток соответствует исходному бенчмарку.

"""
Синусоидальная интерполяция на 10 000 точках: SciPy BPoly.from_derivatives против чистого NumPy (квинтический Эрмит).
Измеряет время построения коэффициентов и время вычисления, а также проверяет численное совпадение.
"""
import sys
import time
import numpy as np
import scipy
from scipy.interpolate import BPoly

# Отчёт об окружении
print(f"Python version: {sys.version.split()[0]}")
print(f"NumPy version: {np.__version__}")
print(f"SciPy version: {scipy.__version__}")

# Сетка и точные значения/производные на [0, 2π]
n_pts = 10_000
grid_x = np.linspace(0.0, 2 * np.pi, n_pts)
f_val = np.sin(grid_x)           # f(x)
d1_val = np.cos(grid_x)          # f'(x)
d2_val = -np.sin(grid_x)         # f''(x)

# === SciPy: BPoly.from_derivatives ===
stacked_y = np.column_stack((f_val, d1_val, d2_val))

t_start = time.perf_counter()
poly_scipy = BPoly.from_derivatives(grid_x, stacked_y)
t_end = time.perf_counter()
scipy_build = t_end - t_start

# Вычисление в исходных узлах
t_start = time.perf_counter()
y_bp = poly_scipy(grid_x)
t_end = time.perf_counter()
scipy_eval = t_end - t_start

# === Pure NumPy: quintic Hermite (per-interval linear solves) ===
def assemble_quintic(xk, f, fp, fpp):
    """
    Вычислить коэффициенты квинтического Эрмита на интервалах [xk[i], xk[i+1]].
    Возвращает массив формы (6, len(xk)-1) с коэффициентами полинома по интервалам.
    """
    m = len(xk) - 1
    params = np.zeros((6, m))
    for i in range(m):
        h = xk[i + 1] - xk[i]
        M = 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],
        ])
        rhs = np.array([f[i], fp[i], fpp[i], f[i + 1], fp[i + 1], fpp[i + 1]])
        params[:, i] = np.linalg.solve(M, rhs)
    return params


def eval_quintic(xk, coeff, xq):
    """
    Вычислить кусочно-квинтический полином в точках xq, используя заранее посчитанные coeff.
    xk: точки разбиения (n,), coeff: (6, n-1), xq: набор точек.
    """
    idx = np.searchsorted(xk, xq) - 1
    idx = np.clip(idx, 0, len(xk) - 2)
    dx = xq - xk[idx]
    out = np.zeros_like(xq)
    for j in range(6):
        out += coeff[j, idx] * dx**j
    return out

# Построение пользовательских коэффициентов
t_start = time.perf_counter()
quintic_coef = assemble_quintic(grid_x, f_val, d1_val, d2_val)
t_end = time.perf_counter()
np_build = t_end - t_start

# Вычисление в узлах
t_start = time.perf_counter()
y_np = eval_quintic(grid_x, quintic_coef, grid_x)
t_end = time.perf_counter()
np_eval = t_end - t_start

# Проверка совпадения на сетке
assert np.allclose(y_bp, y_np, atol=1e-12), "NumPy and SciPy interpolants differ"
print("Custom NumPy interpolation matches SciPy BPoly for 10 000 points.")

# Отчёт по времени
print(f"SciPy coeff time: {scipy_build:.6f} s")
print(f"SciPy eval time : {scipy_eval:.6f} s")
print(f"Custom coeff time: {np_build:.6f} s")
print(f"Custom eval time : {np_eval:.6f} s")

# Необязательно: визуализация (не требуется для замеров)
import matplotlib.pyplot as plt
plt.plot(grid_x, y_bp, label='SciPy BPoly')
plt.plot(grid_x, y_np, '--', label='Custom NumPy')
plt.legend()
plt.title('Sinusoidal interpolation: SciPy vs NumPy Quintic')
plt.xlabel('x')
plt.ylabel('sin(x)')
plt.show()

В одном из зафиксированных запусков вывод показал, что самописный вариант на NumPy даёт те же значения в узлах и быстрее строит коэффициенты. Тест выполнялся на Python 3.11.7, NumPy 2.1.3 и SciPy 1.15.2.

Что на самом деле даёт разрыв в производительности

BPoly.from_derivatives не получал заметных доработок после начальной реализации (которая, в своё время, заменила интерполятор, более медленный на порядок). На то были разные причины, одна из них — его во многом заменяет CubicHermiteSpline. В реальных проектах мы почти не видели его использования, поэтому и мотивации уделять ему внимание было мало. Если вы им пользуетесь или стали бы пользоваться, будь он быстрее, не стесняйтесь прислать pull request. Мы с радостью его рассмотрим.

Ключевой момент в том, что BPoly.from_derivatives — это общий путь, который не оптимизировали целенаправленно, в том числе потому, что его в основном заменил CubicHermiteSpline и он редко используется. Наблюдаемый разрыв — следствие инженерных приоритетов, а не ограничение всей системы интерполяции SciPy.

Дополнительное наблюдение автора бенчмарка: дальнейшая векторизация построителя поинтервальных коэффициентов позволяет убрать явный цикл, однако прирост на протестированных системах был ничтожным. Это не меняет ключевого вывода о сравнительной производительности для данной задачи.

Практический курс: используйте оптимизированный путь

Когда постановка задачи позволяет, более специализированный и поддерживаемый CubicHermiteSpline — современная альтернатива BPoly.from_derivatives. Его прямо называют преемником, и он получает больше внимания. Ниже показано, как построить и вычислить его на той же синусоиде, используя значения функции и первые производные.

from scipy.interpolate import CubicHermiteSpline

# Построить кубический эрмитов сплайн по значениям и первым производным
start = time.perf_counter()
chs = CubicHermiteSpline(grid_x, f_val, d1_val)
stop = time.perf_counter()
chs_build = stop - start

# Вычисление в исходных узлах
start = time.perf_counter()
chs_y = chs(grid_x)
stop = time.perf_counter()
chs_eval = stop - start

print(f"CubicHermiteSpline coeff time: {chs_build:.6f} s")
print(f"CubicHermiteSpline eval time : {chs_eval:.6f} s")

Опыт, приведённый вместе с сравнением, согласуется с этой рекомендацией.

Это правильный ответ. Я ещё больше оптимизировал реализацию numpy-numba, сверх того, что сделал @JeromeRichard в своём ответе. При значительных усилиях я вышел на сопоставимый уровень с CubicHermiteSpline, но на моей машине CubicHermiteSpline всё же оказался быстрее. Я бы не рекомендовал обновлять более общий from_derivatives. Было бы полезно поправить документацию SciPy, в частности добавить указатель на более оптимизированные варианты при рассмотрении общих опций интерполяции. Команда SciPy, отличная работа!

Почему это важно инженерам

Выводы о производительности зависят от того, какое именно API вы используете. Универсальная рутина, рассчитанная на широкие возможности, может быть настроена хуже, чем узкоспециализированный конструктор сплайна, а компактная реализация на NumPy легко обгонит неотрефакторенный код, когда математика проста и предсказуема. Понимание того, какие функции в библиотеке активно поддерживаются, а какие — наследие или заменены, напрямую влияет на результат ваших бенчмарков и, в конечном итоге, на производственные решения.

Выводы и рекомендации

Если вы видите большой разрыв между самописным интерполятором на NumPy и SciPy BPoly.from_derivatives, скорее всего, это отражает отсутствие недавних оптимизаций в этой конкретной функции. Предпочитайте CubicHermiteSpline, когда он подходит под ваши ограничения. Если вашему процессу действительно нужен BPoly.from_derivatives и требуется больше скорости, подумайте о вкладах в апстрим — мейнтейнеры готовы рассматривать pull request’ы. Для локальной доводки простой NumPy может быть весьма эффективен, а более глубокие микрооптимизации поверх ясной реализации могут мало что изменить — в зависимости от системы.