Фиксация граничных значений на сплайне?

Одна из идей для подбора этой функции — сказать, что f(0) = 0 и f(1) = 1 — это просто наблюдения за вашими функциями. Конечно, это наблюдения, в которых вы гораздо более уверены, чем в остальных, но это все равно всего лишь наблюдения.

UnivariateSpline принимает аргумент w, чтобы придать конкретным наблюдениям больший вес. Одним из способов «обеспечения соблюдения» этих граничных условий было бы добавление этих точек как обычных, но с высоким весом.

import numpy as np
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt
## Generate noisy evaluations of the square root function.
## (In my exmaple, I don't have a closed form expression for my function)

x = np.linspace(0,1,1000)
y = np.sqrt(x) + (np.random.random(1000)-0.5)*0.05


# Add fixpoints to function
fixpoint_weight = 1000
fixpoints = np.array([
    #X, Y
    (0, 0),
    (1, 1),
]).T
x_with_fixpoints = np.concatenate([x, fixpoints[0]])
y_with_fixpoints = np.concatenate([y, fixpoints[1]])
weights = np.ones(len(x_with_fixpoints))
weights[-fixpoints.shape[1]:] = fixpoint_weight
# x must be increasing - sort x, y, and weights by x
sort_order = x_with_fixpoints.argsort()
x_with_fixpoints = x_with_fixpoints[sort_order]
y_with_fixpoints = y_with_fixpoints[sort_order]
weights = weights[sort_order]

## Fit spline
spline = UnivariateSpline(x_with_fixpoints, y_with_fixpoints, w=weights, s=0.5)

plt.figure(figsize=(7,7))
plt.scatter(x,y,s=1,label = "Monte Carlo samples")
plt.plot(x,spline(x),color='red',label = "spline")
plt.legend()
plt.title("Noisy evaluations of sqrt")
plt.grid()
plt.show()
print(spline(0) - 0, spline(1) - 1)

Вам нужна функция, которая создает фиксированный B-сплайн, который интерполирует заданные конечные точки. К сожалению, насколько мне известно, ни один scipy-сплайновый интерфейс не реализует это в точности.

Тем не менее, в разделе 9 этого онлайн-курса описан алгоритм решения именно этой задачи, обозначенный как « Глобальная аппроксимация кривой». Ниже я реализую версию этого алгоритма на Python, которая соответствует описанию на веб-сайте.

Функция bspline_least_square_with_clamped_endpoints(x, y, n_ctrl, p) возвращает объект scipy.interpolate.BSpline с коэффициентами n_ctrl, который интерполирует первую и последнюю входные точки по вашему запросу. Входные данные x,y такие же, как в вашем примере, а параметр p обозначает степень сплайна, которая по умолчанию равна 3. Параметр n_ctrl — это количество коэффициентов в результирующем B-сплайне, которое вы можете настроить под свои нужды. Обратите внимание, что по мере увеличения результата результат может переобуться.

Вызов функции, как показано ниже, при входе, с дополнительными точками (0,0) и (1,1) в начале и конце, а также с n_ctrl=6, приводит к следующему рисунку.

spline = bspline_least_square_with_clamped_endpoints(x, y, n_ctrl=6, p=3)

Реализация соответствует объяснению на веб-сайте.

from scipy.interpolate import BSpline
from scipy.linalg import solve


def bspline_least_square_with_clamped_endpoints(x, y, n_ctrl, p=3):
    """
    Build a clamped BSpline least square approximation of the points (x, y) with
    number of coefficients=n_ctrl and spline_degree=p.
    The resulting BSpline will satisfy BSpline(x[0]) == y[0] and BSpline(x[-1]) == y[-1].
    Based on the algorithm presented in: https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-APP-global.html
    """

    # Build clamped knot vector between x[0] and x[-1], (simple) implementation here builds uniform inner knots
    num_inner_knots_inc_end_pts = n_ctrl - p + 1  # including knots[p] and knots[-p-1]
    inner_knots = np.linspace(x[0], x[-1], num_inner_knots_inc_end_pts)
    knots = np.concatenate([np.array([x[0]] * p), inner_knots, np.array([x[-1]] * p)])

    N = build_basis_functions_eval_matrix(knots, x, p)
    Q = y - N[:, 0] * y[0] - N[:, -1] * y[-1]

    # now remove first and last rows and columns since we only need rows 1 to n-1, columns 1 to h-1
    N = N[1:-1, 1:-1]
    Q = Q[1:-1]

    c = solve((N.T).dot(N), (N.T).dot(Q))
    c = np.concatenate([[y[0]], c, [y[-1]]])  # append the first and last values (as the interpolating end coefficients)

    bs = BSpline(knots, c, p)
    return bs

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

Основная алгоритмическая работа выполняется в функции build_basis_functions_eval_matrix(knots, x, p) ниже, которая строит матрицу N оценок базисной функции B-сплайна на входах. Форма матрицы равна (len(x), n_ctrl), а запись N[k, i] содержит оценку базисной функции N_i(x_k). В своей реализации я использую тот факт, что i-я базисная функция B-сплайна может быть восстановлена ​​из scipy.interpolate.BSpline, установив коэффициенты в [0,..,0,1,0,...,0] - см. мой ответ здесь, где обсуждается, как это делается.

def build_basis_functions_eval_matrix(knots, x, p=3):
    """ Build N matrix, where N[k, i] contains the evaluation of basis function N_i(x_k)."""

    n = len(knots) - p - 1  # number of ctrls == number of columns
    m = len(x)  # number of samples == number of rows

    # Using the hack that scipy.interpolate.BSpline with coefficients [0,...,0, 1, 0,...,0]
    # reconstructs the i'th basis function (see https://stackoverflow.com/a/77962609/9702190).
    cols = []
    for i in range(n):
        c = np.zeros((n,))
        c[i] = 1
        N_i = BSpline(knots, c, p)
        col = N_i(x)
        cols.append(col.reshape((m, 1)))
    N = np.hstack(cols)
    return N