GCV в make_smoothing_spline: почему lambda скачет и что делать

Разбираем непоследовательное сглаживание в SciPy: GCV в make_smoothing_spline, выбор lambda, влияние масштаба x и прием со шкалированием к единичному шагу

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

Воспроизводим непоследовательное сглаживание

Ниже приведённый фрагмент загружает две нормализованные серии, которые визуально почти не различаются. Обе сглаживаются с помощью make_smoothing_spline с выбором lam по умолчанию на основе GCV. Дополнительно выводятся значения lambda, которые возвращает scipy.interpolate._bsplines._compute_optimal_gcv_parameter для двух заранее подготовленных матриц дизайна и штрафных терминов.

import pickle
import matplotlib.pyplot as plt
from scipy.interpolate import make_smoothing_spline
from scipy.interpolate._bsplines import _compute_optimal_gcv_parameter

# Загружаем подготовленные массивы
with open('good_lam.pkl', 'rb') as fh_a:
    t_ok, matX_ok, wE_ok, series_ok, w_ok = pickle.load(fh_a)
with open('bad_lam.pkl', 'rb') as fh_b:
    t_ko, matX_ko, wE_ko, series_ko, w_ko = pickle.load(fh_b)

# Сглаживание с выбором по GCV
spline_ok = make_smoothing_spline(t_ok, series_ok, w=None, lam=None)
yhat_ok = spline_ok(t_ok)
spline_ko = make_smoothing_spline(t_ko, series_ko, w=None, lam=None)
yhat_ko = spline_ko(t_ko)

# Сравниваем lambda, которые использует внутренняя процедура GCV в SciPy
print('lam_ok:', _compute_optimal_gcv_parameter(matX_ok, wE_ok, series_ok, w_ok))
print('lam_ko:', _compute_optimal_gcv_parameter(matX_ko, wE_ko, series_ko, w_ko))

# Визуализация
plt.plot(series_ok, label='ok', color='green')
plt.plot(yhat_ok, label='ok_smooth', color='palegreen')
plt.plot(series_ko, label='ko', color='red')
plt.plot(yhat_ko, label='ko_smooth', color='lightpink')
plt.legend()
plt.show()

Что на самом деле происходит

Причина непоследовательности — чувствительность к масштабу внутри целевой функции GCV в реализации _compute_optimal_gcv_parameter. Для проблемного набора данных фробениусовы нормы матриц, входящих в систему, сильно различаются по величине: норма XtWX порядка 7.5, а XtE — около 1.9e+08. Алгоритм формирует левую часть

lhs = XtWX + lam * XtE

Когда XtE на порядки превосходит XtWX, последняя фактически игнорируется при широком диапазоне lam, и минимум определяется численно плохо обусловленным компромиссом. Итог сильно зависит от шага по x. На практике изменение масштаба x способно изменить форму кривой GCV: в одном наборе данных наблюдалось несколько локальных минимумов, которые исчезли после умножения x на 100 и пересчёта промежуточных массивов. Это указывает на разрежение/шаг по x как на значимый фактор. Кроме того, есть признаки, что в некоторых сериях шум неслучаен, и тогда GCV может обосновать крайне малое lam — по сути, отказ от сглаживания.

Практический приём: подбирать GCV на x с единичным шагом и масштабировать lambda обратно

Простой способ повысить численную устойчивость — перед вызовом процедуры GCV перейти к параметризации с единичным шагом по x, а затем вернуть полученное lam к исходному масштабу. Равномерная сетка по x уменьшает разрыв между нормами XtWX и XtE до куда более разумного диапазона (например, 16 и 40 в одном из тестов), благодаря чему оптимизация ведёт себя стабильнее.

Преобразование несложное: разделите x на средний шаг, выполните расчёт GCV на пересчитанной сетке и умножьте итоговую lambda на x_spacing_avg ** 3, чтобы вернуть её к исходному x. В таком подходе для проблемного набора данных и уточнённый поиск GCV, и метод с пересчётом сходятся примерно к 2.63e-08.

import numpy as np
from scipy.interpolate._bsplines import _compute_optimal_gcv_parameter
from scipy._lib._array_api import array_namespace
from scipy.interpolate._bsplines import _coeff_of_divided_diff
from scipy.interpolate import BSpline


def pick_gcv_lambda_stable(t_vals, y_vals):
    d = np.diff(t_vals)
    assert (d >= 0).all(), 'x must be sorted'
    mean_step = d.mean()
    assert mean_step != 0, 'div by zero'

    t_unit = t_vals / mean_step
    Phi, WEmat, y_arr, w_vec = assemble_spline_arrays(t_unit, y_vals)
    lam_unit = _compute_optimal_gcv_parameter(Phi, WEmat, y_arr, w_vec)
    return lam_unit * mean_step ** 3


def assemble_spline_arrays(x_in, y_in, w_in=None):
    ax = 0
    _ = array_namespace(x_in, y_in)

    x = np.ascontiguousarray(x_in, dtype=float)
    y = np.ascontiguousarray(y_in, dtype=float)

    if any(x[1:] - x[:-1] <= 0):
        raise ValueError('``x`` should be an ascending array')

    if x.ndim != 1 or x.shape[0] != y.shape[ax]:
        raise ValueError(f'``x`` should be 1D and {x.shape = } == {y.shape = }')

    if w_in is None:
        w = np.ones(len(x))
    else:
        w = np.ascontiguousarray(w_in)
        if any(w <= 0):
            raise ValueError('Invalid vector of weights')

    t_knots = np.r_[[x[0]] * 3, x, [x[-1]] * 3]
    n = x.shape[0]

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

    y = np.moveaxis(y, ax, 0)
    y_tail = y.shape[1:]
    if y_tail != ():
        y = y.reshape((n, -1))

    # Матрица признаков в базисе B-сплайнов
    bx = BSpline.design_matrix(x, t_knots, 3)
    Phi = np.zeros((5, n))
    for k in range(1, 4):
        Phi[k, 2:-2] = bx[k: k - 4, 3:-3][np.diag_indices(n - 4)]

    Phi[1, 1] = bx[0, 0]
    Phi[2, :2] = ((x[2] + x[1] - 2 * x[0]) * bx[0, 0], bx[1, 1] + bx[1, 2])
    Phi[3, :2] = ((x[2] - x[0]) * bx[1, 1], bx[2, 2])
    Phi[1, -2:] = (bx[-3, -3], (x[-1] - x[-3]) * bx[-2, -2])
    Phi[2, -2:] = (bx[-2, -3] + bx[-2, -2], (2 * x[-1] - x[-2] - x[-3]) * bx[-1, -1])
    Phi[3, -2] = bx[-1, -1]

    WEmat = np.zeros((5, n))
    WEmat[2:, 0] = _coeff_of_divided_diff(x[:3]) / w[:3]
    WEmat[1:, 1] = _coeff_of_divided_diff(x[:4]) / w[:4]
    for j in range(2, n - 2):
        WEmat[:, j] = (x[j+2] - x[j-2]) * _coeff_of_divided_diff(x[j-2:j+3]) / w[j-2:j+3]
    WEmat[:-1, -2] = -_coeff_of_divided_diff(x[-4:]) / w[-4:]
    WEmat[:-2, -1] = _coeff_of_divided_diff(x[-3:]) / w[-3:]
    WEmat *= 6
    return Phi, WEmat, y, w

Получив пересчитанное lam, используйте его напрямую при построении сглаживающего сплайна на исходных x. Больше в конвейере менять ничего не требуется.

from scipy.interpolate import make_smoothing_spline

lam_fixed = pick_gcv_lambda_stable(t_ko, series_ko)
spline_fixed = make_smoothing_spline(t_ko, series_ko, w=None, lam=lam_fixed)
yhat_fixed = spline_fixed(t_ko)

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

Цель GCV — найти баланс между точностью и сглаживанием, но численные свойства задействованных матриц способны этот баланс нарушить. Когда шаг по x делает XtE на многие порядки больше XtWX, оптимизатор может попасть в плохие локальные минимумы или фактически игнорировать часть целевой функции. Поэтому внешне похожие, нормализованные ряды дают резко разные lambda и последующие аппроксимации. Перевод x к единичному шагу регуляризует вычисление, не меняя модель и её допущения, а обратное масштабирование lam кубом среднего шага сохраняет согласованность решения с исходными данными.

Бывают случаи, когда GCV корректно возвращает очень малое lambda и почти не сглаживает. Это возможно, если вариативность ряда — не случайный шум; в таких условиях GCV предпочитает сохранить структуру, а не подавлять её. Отдельно отмечалось, что изменение шага по x может устранить лишние локальные минимумы на кривой GCV. Это поведение обсуждалось и дополнительно исследовалось; см. публичные обсуждения и отчёты об ошибках на github.com/scipy/scipy/issues, в частности https://github.com/scipy/scipy/issues/23472.

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

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

С этими корректировками вы сможете по‑прежнему автоматически подбирать lam через make_smoothing_spline и при этом улучшить численное поведение на наборах данных с разной частотой выборки и масштабами.