Кубические сплайны

Урок про кубические сплайны — стандартный способ гладко интерполировать данные без бед высокой степени.

Кубический сплайн — кусочная функция, на каждом отрезке между узлами равная своему полиному 3-й степени, причём в узлах склейки совпадают значения, первые и вторые производные — кривая получается гладкой.

Идея: много маленьких кубиков

Феномен Рунге наказывает за высокую степень. Решение радикальное: вообще не повышать степень. Вместо одного полинома степени n через все узлы построим на каждом отрезке между соседними узлами свой полином 3-й степени, а затем гладко склеим их. Степень всегда 3, сколько бы ни было узлов — Рунге исключён. Кубик выбран потому, что это минимальная степень, позволяющая обеспечить гладкость до второй производной (непрерывную кривизну) — глаз воспринимает такую кривую как «идеально плавную».

Условия склейки

На n отрезках — n кубиков, у каждого 4 коэффициента, итого 4n неизвестных. Их задают условия:

Эти условия сводятся к трёхдиагональной системе на вторые производные в узлах — и решаются прогонкой за O(n) (вот где пригодился метод из раздела про СЛАУ!).

def кубический_сплайн(xs, ys):
    # натуральный сплайн: S''=0 на концах. Возвращает коэффициенты b,c,d
    n = len(xs) - 1
    h = [xs[i+1] - xs[i] for i in range(n)]
    alpha = [0.0] * (n + 1)
    for i in range(1, n):
        alpha[i] = 3*((ys[i+1]-ys[i])/h[i] - (ys[i]-ys[i-1])/h[i-1])
    # прогонка (трёхдиагональная система на c = S''/2)
    l = [1.0] + [0.0]*n; mu = [0.0]*(n+1); z = [0.0]*(n+1)
    for i in range(1, n):
        l[i] = 2*(xs[i+1]-xs[i-1]) - h[i-1]*mu[i-1]
        mu[i] = h[i]/l[i]
        z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i]
    c = [0.0]*(n+1); b = [0.0]*n; d = [0.0]*n
    for j in range(n-1, -1, -1):
        c[j] = z[j] - mu[j]*c[j+1]
        b[j] = (ys[j+1]-ys[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3
        d[j] = (c[j+1]-c[j])/(3*h[j])
    return b, c, d

def значение_сплайна(xs, ys, b, c, d, x):
    i = 0
    while i < len(xs)-2 and x > xs[i+1]:
        i += 1
    dx = x - xs[i]
    return ys[i] + b[i]*dx + c[i]*dx*dx + d[i]*dx*dx*dx

xs = [0.0, 1.0, 2.0, 3.0]
ys = [0.0, 1.0, 0.0, 1.0]     # «зигзаг», полином 3-й степени дал бы выброс
b, c, d = кубический_сплайн(xs, ys)
print("проверка узлов:", [round(значение_сплайна(xs, ys, b, c, d, t), 4) for t in xs])
for x in [0.5, 1.5, 2.5]:
    print(f"  S({x}) = {значение_сплайна(xs, ys, b, c, d, x):.6f}")

Вывод:

проверка узлов: [0.0, 1.0, -0.0, 1.0]
  S(0.5) = 0.750000
  S(1.5) = 0.500000
  S(2.5) = 0.250000

Сплайн точно проходит через все узлы и плавно интерполирует между ними, без выбросов. На функции Рунге он дал бы маленькую ошибку там, где полином высокой степени «взрывался».

Как работает под капотом

Натуральный сплайн (с S''=0 на концах) — лишь один из вариантов краевых условий. Бывают: зажатый (clamped) — задаём наклон на концах; «not-a-knot» — требуем, чтобы у первых двух (и последних двух) кубиков совпадала и третья производная; именно его по умолчанию использует scipy.interpolate.CubicSpline, так как он точнее у краёв. Замечательное свойство сплайна: среди всех дважды дифференцируемых интерполянтов натуральный кубический сплайн минимизирует интеграл от квадрата второй производной — то есть он «самый негибкий», самый плавный из возможных. Именно поэтому он так хорош визуально и используется в графике, CAD и шрифтах.

Частые ошибки

• Путать сплайн с регрессией. Сплайн проходит через узлы (интерполяция); для зашумлённых данных это плохо — там нужно сглаживание (МНК/сглаживающий сплайн).
• Игнорировать краевые условия. У краёв разные условия дают заметно разный результат; для гладких данных «not-a-knot» обычно лучше натурального.
• Решать систему сплайна плотным Гауссом. Она трёхдиагональная — используйте прогонку O(n).

Итоги

• Кубический сплайн — гладкая склейка кубиков по отрезкам; степень всегда 3, поэтому феномена Рунге нет.
• В узлах совпадают значения, первые и вторые производные — кривая плавная (непрерывная кривизна).
• Коэффициенты находят из трёхдиагональной системы прогонкой за O(n).
• Натуральный сплайн минимизирует «изогнутость» — отсюда его непревзойдённая гладкость; на практике часто берут «not-a-knot».