Интерполяция сплайном

Определение сплайна

Сплайн - это кусочно-многочленная функция. Для сетки \(\left[x_0, x_1, \ldots, x_n\right]\) сплайн задается в виде \(n\) многочленов \(s_i(x)\):

\[ s(x)=s_i(x), \quad x \in\left[x_{i-1}, x_i\right], \quad i=1,2, \ldots, n \]

В узлах интерполяции \(x_i\) сплайн принимает заданные значения \(f\left(x_i\right)\):

\[ s\left(x_i\right)=s_i\left(x_i\right)=s_{i+1}\left(x_i\right)=f\left(x_i\right) \]

Характеристики сплайна

Сплайн характеризуется следующими параметрами:

• Степень - это степень многочленов \(s_i(x)\)

• Гладкость - это количество непрерывных производных у \(s(x)\) (в точках склейки)

• Дефект - это разность между степенью и гладкостью.

Легко показать, что условие «дефект \(=0\) » приводит к тому, что все \(s_i(x)\) совпадают, а сплайн превращается в интерполяционный многочлен.

Интерполяционный полином Эрмита

Как мы понимаем, чтобы сплайн имел нужную нам гладкость, нам надо «сшить» производные в узлах у соседних многочленов. Для такого интерполяционного полинома (Лагранжа) не хватит.

Пусть имеем узлы \(\{ x_i \}_{i=0}^s\), на которых будем строить многочлен. В этих узлах надо знать значение функции и некоторое количество её производных - обозначим как \(\{ f_i^{(j)} \}_{j=0}^{\alpha_i - 1}\) для конкретного узла \(x_i\). Здесь, например, \(f_i^{(0)}=f(x_i)\) - известное значение функции в точке \(x_i\), а \(f_i^{(1)}=f'(x_i)\) - известное значение производной в точке \(x_i\).

Алгебраический многочлен степени не выше \(m=\alpha_{0}+\alpha_{1}+\cdots+\alpha_{s}-1\), имеющий в узлах \(x_i\) нужные производные \(f_i^{(j)}\), называется интерполяционным многочленом Эрмитa \(H_{m}\). Он может иметь следующее представление:

\[ H_{m}(x)=\sum_{i=0}^{s} \sum_{j=0}^{a_{i}-1} \sum_{k=0}^{a_{i}-j-1} f_{i}^{(j)}\frac{1}{k ! \cdot j !}\left[\frac{\left(x-x_{i}\right)^{\alpha_{i}}}{\left(x-x_{0}\right)^{\alpha_{0}} \ldots\left(x-x_{s}\right)^{\alpha_{s}}}\right]_{x=x_{i}}^{(k)} \frac{\left(x-x_{0}\right)^{\alpha_{0}} \ldots\left(x-x_{s}\right)^{\alpha_{s}}}{\left(x-x_{i}\right)^{\alpha_{i}-j-k}} . \]

Остаточный многочлен для интерполяционного многочлена Эрмита имеет вид

\[ R_{m}(x)=\frac{f_{x}^{(m+1)}(\xi)}{(m+1) !}\left(x-x_{0}\right)^{\alpha_{0}} \cdots\left(x-x_{s}\right)^{\alpha_{s}}, \]

где \(\xi \in\left(x_{0}, x_{s}\right)\).

Несложно понять по смыслу многочлена Эрмита, что для случая \(s=n\) (берём все имеющиеся узлы) и \(\alpha_i=1\) (нужно только одно число - значение функции в точке, никаких производных) получается, что интерполяционный многочлен Эрмита становится интерполяционным многочленом Лагранжа.

Таким образом, чтобы задать сплайн, можно на каждом \([x_{i}, x_{i+1}]\) задать полином Эрмита нужной степени \(m\), а производные взять из условия равенства соседним. Так останется только явно задать производные на границах отрезка.

Пример построения кубического сплайна дефекта 1

Существует несколько способов определения функций \(s_i(x)\) для кубического сплайна, все они сводят задачу к решению трехдиагональной системы.

Можем полностью определить \(s_i(x)\) из условий

\[\begin{split} \begin{aligned} s_i\left(x_{i-1}\right) & =f_{i-1} \\ s_i\left(x_i\right) & =f_i \\ s_i^{\prime}\left(x_{i-1}\right) & =m_{i-1} \\ s_i^{\prime}\left(x_i\right) & =m_i \end{aligned} \end{split}\]

где \(m_i\) - заданная производная функции f(x) в узлах \(x_i\). Эта задача называется задачей Эрмитовой интерполяции.

По заданным константам и тем фактом, что \(s_i(x)\) - кубическая парабола, можем выписать саму формулу (вспоминаем разделенные разности):

\[ s_i(x)=f_i+m_i\left(x-x_i\right)+\frac{2 m_i+m_{i-1}-3 f\left(x_{i-1}, x_i\right)}{x_i-x_{i-1}}\left(x-x_i\right)^2+\frac{m_i+m_{i-1}-2 f\left(x_{i-1}, x_i\right)}{\left(x_i-x_{i-1}\right)^2}\left(x-x_i\right)^2 \]

Условие на гладкость 2 означает требование \(s_i^{\prime \prime}\left(x_i\right)=s_{i+1}^{\prime \prime}\left(x_i\right)\), что переходит в

\[ \frac{2 m_i+m_{i-1}-3 f\left(x_{i-1}, x_i\right)}{h_{i-1 / 2}}=-\frac{2 m_i+m_{i+1}-3 f\left(x_i, x_{i+1}\right)}{h_{i+1 / 2}} \]

упрощяя,

\[ \frac{m_{i-1}}{h_{i-1 / 2}}+\left(\frac{2}{h_{i-1 / 2}}+\frac{2}{h_{i+1 / 2}}\right) m_i+\frac{m_{i+1}}{h_{i+1 / 2}}=3\left(\frac{f\left(x_{i-1}, x_i\right)}{h_{i-1 / 2}}+\frac{f\left(x_i, x_{i+1}\right)}{h_{i+1 / 2}}\right) \]

что является трехдиагональной системой при \(i=1, \ldots, n-1\). Здесь \(h_{i-1 / 2}=x_i-x_{i-1}\). Вспоминаем метод прогонки - нужны ещё первая и последние строчки системы - граничные условия.

Для естественного сплайна

\[ \frac{s_1^{\prime \prime}\left(x_0\right)}{2}=-\frac{2 m_1+m_0-3 f\left(x_0, x_1\right)}{h_{1 / 2}}=0, \quad \frac{s_n^{\prime \prime}\left(x_n\right)}{2}=\frac{2 m_n+m_{n-1}-3 f\left(x_{n-1}, x_n\right)}{h_{n-1 / 2}}=0 \]

дают граничные уравнения в системе

\[\begin{split} \begin{aligned} \frac{m_0+2 m_1}{h_{1 / 2}} & =3 \frac{f\left(x_0, x_1\right)}{h_{1 / 2}} \\ \frac{2 m_n+m_{n-1}}{h_{n-1 / 2}} & =3 \frac{f\left(x_{n-1}, x_n\right)}{h_{n-1 / 2}} \end{aligned} \end{split}\]

Матрица такой системы будет симметричной и положительно определенной. А значит все \(m_i\) определяются методом прогонки. Итого, все компоненты сплайна определены.

Сплайн на практике через scipy

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 4, 30)
y = np.sin(x**3/3+4)

x_dens = np.linspace(0, 4, 200)
y_dens = np.sin(x_dens**3/3+4)

from scipy import interpolate
f1 = interpolate.interp1d(x, y, kind = 'linear') # Определяем вид интерполяции, т.е. степень многочлена
f2 = interpolate.interp1d(x, y, kind = 'cubic') # И что удобно, это сразу же функции


plt.plot(x, y, 'o', x_dens, y_dens, 'g-', x_dens, f1(x_dens), '-', x_dens, f2(x_dens), 'r--')
plt.legend(['data','true func', 'linear', 'cubic'], loc = 'best')
plt.show()

Сплайнить можно не только изначальные данные

from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt


W = np.array([0.435,0.606,0.814,1.05,1.25,1.40,1.60])
sum_all = np.array([  3.87282732e+32,   8.79993191e+32,   1.74866333e+33,   1.59946687e+33,
   9.08556547e+33,   6.70458731e+33,   9.84832359e+33])
sum_can = np.array([  2.98381061e+28,   1.26194810e+28,   3.30328780e+28,   2.90254609e+29,
   3.65117723e+29,   3.46256846e+29,   3.64483736e+29])

fall = CubicSpline(W,sum_all)
newallx=np.arange(0.435,1.6,0.001)
newally=fall(newallx)



fcan = CubicSpline(W,sum_can)
newcanx=np.arange(0.435,1.6,0.001)
newcany=fcan(newcanx)



#----plot



plt.plot(newallx,newally)
plt.plot(newcanx,newcany)
plt.plot(W,sum_all,marker='o',color='r',linestyle='')
plt.plot(W,sum_can,marker='o',color='b',linestyle='')
plt.yscale("log")
plt.ylabel("Flux S$_v$ [erg s$^-$$^1$ cm$^-$$^2$ Hz$^-$$^1$]")
plt.xlabel("Wavelength [n$\lambda$]")
plt.show()

Здесь показана довольно типичная ошибка, которая появляется при подобной обработке данных. Сплайн проходит ниже 0, что не определено в масштабе шкалы.

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

from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt


W = np.array([0.435,0.606,0.814,1.05,1.25,1.40,1.60])
sum_all = np.array([  3.87282732e+32,   8.79993191e+32,   1.74866333e+33,   1.59946687e+33,
   9.08556547e+33,   6.70458731e+33,   9.84832359e+33])
sum_can = np.array([  2.98381061e+28,   1.26194810e+28,   3.30328780e+28,   2.90254609e+29,
   3.65117723e+29,   3.46256846e+29,   3.64483736e+29])



fall = CubicSpline(W,np.log10(sum_all))
newallx=np.arange(0.435,1.6,0.001)
newally=fall(newallx)



fcan = CubicSpline(W,np.log10(sum_can))
newcanx=np.arange(0.435,1.6,0.01)
newcany=fcan(newcanx)



plt.plot(newallx,10**newally)
plt.plot(newcanx,10**newcany)
plt.plot(W,sum_all,marker='o',color='r',linestyle='')
plt.plot(W,sum_can,marker='o',color='b',linestyle='')
plt.yscale("log")



plt.ylabel("Flux S$_v$ [erg s$^-$$^1$ cm$^-$$^2$ Hz$^-$$^1$]")
plt.xlabel("Wavelength [n$\lambda$]")
plt.show()

UnivariateSpline для регрессии - лучшее от обоих миров

На практике, в данных часто есть сильный шум, так что сплайн, не смотря на свою простоту в реализации, не подходит для задачи аппроксимации, в отличие от линейной или не очень регрессии. Однако её надо тонко настраивать, подбирая и обучая модель.

Было бы идеально, если бы существовал вид сплайна, который мог бы немного отклоняться от точек, сохраняя свойства сплайна иметь нужные производные. Такое уже реализовано и называется UnivariateSpline - это кубический сплайн, который может откланяться от целевых точек.

Зачастую, это самый простой способ что-то зафитить. Главное, поиграться с параметром s - ожидаемый MSE.

О работе со сплайном

UnivariateSpline всё ещё остаётся сплайном. Если мы будем брать от него производные, рано или поздно он занулится, а непосредственно до этого станет кусочно-линейной функцией. Например, для кубического сплайна, первая производная будет сплайном степени 2, а вот вторая уже сплайном степени 1 (кусочно-линейная функция).

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

# Сгенерируем искусственные данные: небольшая функция + шум
np.random.seed(0)
x = np.linspace(0, 10, 50)
y = np.sin(x) + 0.1 * np.random.randn(len(x))

# Четыре разных значения параметра s (сильная недофит -> сильная переобученность)
s_values = [1e-6, 0.5, 2.5, 10]  # s=1e-6 - "слишком маленький", s=10 - "слишком большой", два других "где-то посередине"

fig, axes = plt.subplots(2, 2, figsize=(10, 8))
axes = axes.ravel()  # Раскладываем двумерный массив осей в один

for i, s in enumerate(s_values):
    spline = UnivariateSpline(x, y, s=s)
    
    # Рисуем исходные точки
    axes[i].scatter(x, y, label='Данные', color='red', s=20, alpha=0.6)
    
    # Рисуем аппроксимирующую кривую
    xs = np.linspace(0, 10, 200)
    ys = spline(xs)
    axes[i].plot(xs, ys, label=f'Spline, s={s}', color='blue')
    
    axes[i].set_title(f'UnivariateSpline, s={s}')
    axes[i].legend()
    axes[i].set_xlabel('x')
    axes[i].set_ylabel('y')

plt.tight_layout()
plt.show()