Регулировка сглаживающего сплайна в SciPy: смещение lam по GCV

Почему make_smoothing_spline в SciPy не возвращает lam и как масштабировать GCV-параметр: smooth_bias помогает управлять сглаживанием, сохранив автонастройку.

Контроль степени того, насколько «агрессивно» сглаживающий сплайн подгоняется под данные, часто критичен. Функция SciPy make_smoothing_spline() автоматически подбирает параметр регуляризации lam (по критерию GCV), но не предоставляет это значение для последующей ручной настройки. Если вы хотите сместить результат в сторону более сильного сглаживания или, наоборот, большей точности, нужен способ слегка подкорректировать lam, не отказываясь от автоматического выбора.

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

В make_smoothing_spline() значение lam по GCV вычисляется внутри и не возвращается. Значит, вы не можете умножить его на коэффициент постфактум, чтобы изменить баланс между гладкостью и следованием данным.

from scipy.interpolate import make_smoothing_spline

# GCV выбирает lam внутри функции
spline_obj = make_smoothing_spline(x_values, y_values)

# Здесь нет доступа к автоматически найденному `lam`, чтобы масштабировать его
# (оно не сохраняется и не возвращается make_smoothing_spline)

Просмотр исходников подтверждает: lam нигде не сохраняется так, чтобы его можно было извлечь после вызова, то есть публичного способа масштабировать его постфактум нет. Если вы хотите оставить критерий GCV, но при этом управлять результатом, масштабирование нужно внести до решения системы.

Почему так устроено

Функция внутри вычисляет lam, оптимальный по GCV, и сразу решает задачу для коэффициентов сплайна. Значение не сохраняется в возвращаемом объекте и не отдаётся отдельно, поэтому перехватить его нельзя. Иными словами, это намеренное решение: lam считается внутренней деталью реализации.

Решение: переопределить с коэффициентом смещения

Практичный обходной путь — воспроизвести make_smoothing_spline и добавить множитель к автоматически оценённому lam. Коэффициент больше 1 смещает результат к более сильному сглаживанию; значение между 0 и 1, напротив, прижимает аппроксимацию к данным. Единица не меняет поведение. Реализация ниже сохраняет исходную логику, включая оценку по GCV, и умножает полученный lam на выбранный коэффициент.

import numpy as np
from scipy.interpolate._bsplines import _coeff_of_divided_diff as _dd_coef, _compute_optimal_gcv_parameter as _gcv_param
from scipy._lib._array_api import array_namespace as _xp_ns
from scipy.interpolate import BSpline as _BS
from scipy.linalg import solve_banded as _solve_banded
from numpy.core.multiarray import normalize_axis_index as _norm_axis

def fit_spline_gcv_tilt(x_coords, y_data, smooth_bias=1.0, weights=None, alpha=None, *, axis=0):
    """
    Build a cubic smoothing spline using the same logic as SciPy's implementation,
    but multiply the automatically estimated GCV parameter by `smooth_bias`.

    Inputs:
    - smooth_bias >= 0. Values > 1 favor smoother output; values < 1 favor fidelity.
    - alpha: if provided, acts as the regularization parameter directly (non-negative).
    """

    xp_mod = _xp_ns(x_coords, y_data)

    x_arr = np.ascontiguousarray(x_coords, dtype=float)
    y_arr = np.ascontiguousarray(y_data, dtype=float)

    if smooth_bias < 0:
        raise ValueError("`smooth_bias` must be non-negative")

    if any(x_arr[1:] - x_arr[:-1] <= 0):
        raise ValueError("`x` must be strictly increasing")

    if x_arr.ndim != 1 or x_arr.shape[0] != y_arr.shape[axis]:
        raise ValueError(f"`x` must be 1D and {x_arr.shape = } == {y_arr.shape = }")

    if weights is None:
        w_vec = np.ones(len(x_arr))
    else:
        w_vec = np.ascontiguousarray(weights)
        if any(w_vec <= 0):
            raise ValueError("Invalid vector of weights")

    knts = np.r_[[x_arr[0]] * 3, x_arr, [x_arr[-1]] * 3]
    n_pts = x_arr.shape[0]

    if n_pts <= 4:
        raise ValueError("`x` and `y` length must be at least 5")

    axis_idx = _norm_axis(axis, y_arr.ndim)
    y_arr = np.moveaxis(y_arr, axis_idx, 0)
    y_tail_shape = y_arr.shape[1:]
    if y_tail_shape != ():
        y_arr = y_arr.reshape((n_pts, -1))

    # Проектирующая матрица в базисе B-сплайнов (ленточный вид)
    dm = _BS.design_matrix(x_arr, knts, 3)
    band_A = np.zeros((5, n_pts))
    for ii in range(1, 4):
        band_A[ii, 2:-2] = dm[ii: ii - 4, 3: -3][np.diag_indices(n_pts - 4)]

    band_A[1, 1] = dm[0, 0]
    band_A[2, :2] = ((x_arr[2] + x_arr[1] - 2 * x_arr[0]) * dm[0, 0],
                     dm[1, 1] + dm[1, 2])
    band_A[3, :2] = ((x_arr[2] - x_arr[0]) * dm[1, 1], dm[2, 2])
    band_A[1, -2:] = (dm[-3, -3], (x_arr[-1] - x_arr[-3]) * dm[-2, -2])
    band_A[2, -2:] = (dm[-2, -3] + dm[-2, -2],
                      (2 * x_arr[-1] - x_arr[-2] - x_arr[-3]) * dm[-1, -1])
    band_A[3, -2] = dm[-1, -1]

    band_B = np.zeros((5, n_pts))
    band_B[2:, 0] = _dd_coef(x_arr[:3]) / w_vec[:3]
    band_B[1:, 1] = _dd_coef(x_arr[:4]) / w_vec[:4]
    for jj in range(2, n_pts - 2):
        band_B[:, jj] = (x_arr[jj + 2] - x_arr[jj - 2]) * _dd_coef(x_arr[jj - 2: jj + 3]) / w_vec[jj - 2: jj + 3]
    band_B[:-1, -2] = -_dd_coef(x_arr[-4:]) / w_vec[-4:]
    band_B[:-2, -1] = _dd_coef(x_arr[-3:]) / w_vec[-3:]
    band_B *= 6

    # Определить регуляризацию и применить смещение
    if alpha is None:
        lam_val = _gcv_param(band_A, band_B, y_arr, w_vec)
        lam_val *= smooth_bias
    elif alpha < 0.0:
        raise ValueError("Regularization parameter should be non-negative")
    else:
        lam_val = alpha

    # Решение в базисе натуральных сплайнов
    if np.ndim(lam_val) == 0:
        coef = _solve_banded((2, 2), band_A + lam_val * band_B, y_arr)
    elif np.ndim(lam_val) == 1:
        coef = np.empty((n_pts, lam_val.shape[0]))
        for ii in range(lam_val.shape[0]):
            coef[:, ii] = _solve_banded((2, 2), band_A + lam_val[ii] * band_B, y_arr[:, ii])
    else:
        raise RuntimeError("Internal error, please report it to SciPy developers.")

    coef = coef.reshape((coef.shape[0], *y_tail_shape))

    p0, p1 = coef[0:1, ...], coef[1:2, ...]
    pm0, pm1 = coef[-1:-2:-1, ...], coef[-2:-3:-1, ...]

    coeff_full = np.r_[p0 * (knts[5] + knts[4] - 2 * knts[3]) + p1,
                       p0 * (knts[5] - knts[3]) + p1,
                       coef[1: -1, ...],
                       pm0 * (knts[-4] - knts[-6]) + pm1,
                       pm0 * (2 * knts[-4] - knts[-5] - knts[-6]) + pm1]

    knts_xp, coeff_full_xp = xp_mod.asarray(knts), xp_mod.asarray(coeff_full)
    return _BS.construct_fast(knts_xp, coeff_full_xp, 3, axis=axis_idx)

Эта реализация опирается на функции из приватного API SciPy и была успешно проверена на SciPy v1.16.1. Если вам нужно получить эффективное значение lam вместе со сплайном, простой приём — изменить возврат так, чтобы включить это значение, и распаковать его на месте вызова.

# Пример шаблона, если вы решите также возвращать `lam_val` из функции
# return _BS.construct_fast(knts_xp, coeff_full_xp, 3, axis=axis_idx), lam_val
# Тогда вызов будет таким: spline_obj, lam_used = fit_spline_gcv_tilt(...)

Как применять смещение на практике

Если smooth_bias равен 1, поведение полностью совпадает с дефолтным. Увеличение коэффициента ведёт к более гладким кривым, уменьшение — к большей приверженности данным. Вы сохраняете выбор по GCV, лишь масштабируя его перед решением.

# Больше сглаживания, чем по умолчанию
s_smooth = fit_spline_gcv_tilt(x_values, y_values, smooth_bias=2.0)

# Больше соответствия данным, чем по умолчанию
s_faithful = fit_spline_gcv_tilt(x_values, y_values, smooth_bias=0.5)

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

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

Итоги

make_smoothing_spline не предоставляет автоматически выбранный lam, поэтому после вызова функцию изменить его нельзя. Если такая возможность нужна, воспроизведите логику и умножьте GCV-оптимальный параметр на коэффициент. Значения больше единицы отдают приоритет гладкости, меньше единицы — точности. Показанный подход опирается на приватный API SciPy и протестирован на SciPy 1.16.1 — учитывайте это в своей среде.