GCV के साथ smoothing spline में असंगत lambda: कारण और स्थिर उपाय

जानें क्यों GCV SciPy make_smoothing_spline में अस्थिर lambda चुनता है, x को यूनिट स्पेसिंग पर स्केल करें, lambda को मूल स्केल पर लौटाएँ और स्थिर स्मूदिंग पाएं.

जब आप make_smoothing_spline पर भरोसा करते हैं कि वह GCV के जरिये खुद ही स्मूदिंग की ताकत चुने, तो आप समान प्रकार की टाइम सीरीज़ पर एक जैसा व्यवहार अपेक्षित करते हैं। फिर भी, अलग‑अलग सीरीज़ को एक ही पैमाने पर नॉर्मलाइज़ करने के बाद भी, स्मूदिंग में तीखा अंतर दिख सकता है: एक कर्व वाजिब लगता है, जबकि दूसरा लगभग सपाट रेखा में सिमट जाता है। नीचे दिया गया उदाहरण उसी विसंगति को दोहराता है और फिर समझाता है कि यह क्यों होता है तथा मॉडल की धारणाएँ बदले बिना GCV चरण को संख्यात्मक रूप से कैसे अधिक स्थिर बनाया जा सकता है।

असंगत स्मूदिंग को दोहराकर दिखाना

नीचे दिया गया स्निपेट दो नॉर्मलाइज़्ड सीरीज़ लोड करता है जो देखने में बहुत मिलती-जुलती हैं। दोनों को make_smoothing_spline से स्मूद किया गया है, जहाँ lam का चयन डिफ़ॉल्ट GCV-आधारित तरीके से होता है। कोड scipy.interpolate._bsplines._compute_optimal_gcv_parameter द्वारा दो पहले से बने डिज़ाइन मैट्रिसेज़ और पेनल्टी टर्म्स के लिए लौटाए गए lambda मान भी प्रिंट करता है।

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)
# SciPy की GCV रूटीन द्वारा अंदरूनी तौर पर उपयोग किए गए lam मानों की तुलना करें
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 की एक बड़ी रेंज में XtWX व्यावहारिक रूप से नज़रअंदाज़ हो जाता है, और मिनिमाइज़र संख्यात्मक रूप से कम कंडीशन्ड ट्रेड-ऑफ से संचालित होता है। यह x के स्पेसिंग के प्रति बेहद संवेदनशील है। दरअसल, x को रिस्केल करने से GCV कर्व का आकार बदल सकता है: एक डेटासेट में कई स्थानीय मिनिमा दिखे, जो x को 100 से गुणा करके और मध्यवर्ती ऐरे दोबारा गणना करने के बाद गायब हो गए। यह अवलोकन बताता है कि x का स्पेसिंग एक बड़ा कारक है। साथ ही कुछ सीरीज़ में शोर रैंडम नहीं होता, और तब GCV अत्यंत छोटा lam उचित ठहरा सकता है — यानी व्यावहारिक रूप से स्मूदिंग न करना।

व्यावहारिक उपाय: x को इकाई अंतराल में बदलकर GCV चलाएँ और lam को मूल पैमाने पर लौटाएँ

गणना को अधिक स्थिर बनाने का सीधा तरीका है कि GCV चलाने से पहले समस्या को x में यूनिट स्पेसिंग में री-पैरामीटराइज़ कर लें, फिर प्राप्त lam को मूल स्केल पर वापस कन्वर्ट करें। समान रूप से स्पेस किए गए x के साथ XtWX और XtE के नॉर्म के बीच का अंतर काफी समझदारी की सीमा में आ जाता है (उदा., एक केस में 16 और 40), जिससे ऑप्टिमाइज़ेशन ज्यादा पूर्वानुमेय ढंग से चलता है।

कन्वर्ज़न सरल है: x को उसके औसत स्पेसिंग से भाग दें, उस रिस्केल्ड ग्रिड पर GCV गणना करें, और प्राप्त lambda को मूल x पर मैप करने के लिए x_spacing_avg ** 3 से गुणा करें। इसी तरीके से, चुनौतीपूर्ण डेटासेट पर परिष्कृत 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-spline आधार में डिज़ाइन मैट्रिक्स
    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 पर कई क्रमों से भारी बना देता है, तो ऑप्टिमाइज़र खराब स्थानीय मिनिमा में फँस सकता है या ऑब्जेक्टिव के किसी हिस्से को परोक्ष रूप से नज़रअंदाज़ कर सकता है। इसलिए देखने में समान और नॉर्मलाइज़्ड सीरीज़ भी बहुत अलग lam और नतीजतन फिट निकाल सकती हैं। x को यूनिट स्पेसिंग में रिस्केल करना, बिना आधारभूत मॉडल/मान्यताओं को बदले, गणना को नियमित करता है; और औसत स्टेप के घन से lam को वापस स्केल करना समाधान को आपके मूल डेटा के अनुरूप रखता है।

कभी‑कभी GCV सचमुच बहुत छोटा lambda लौटाता है और लगभग कोई स्मूदिंग नहीं करता। यह तब होता है जब सीरीज़ में भिन्नता यादृच्छिक शोर नहीं होती; ऐसी स्थिति में GCV संरचना को बनाए रखना पसंद करेगा, न कि उसे कम करना। अलग से यह भी देखा गया है कि x का स्पेसिंग बदलने से GCV कर्व में अतिरिक्त स्थानीय मिनिमा गायब हो सकते हैं। इस व्यवहार पर चर्चा और आगे की जाँच हुई है; संदर्भ के लिए github.com/scipy/scipy/issues और विशेष रूप से https://github.com/scipy/scipy/issues/23472 देखें।

निष्कर्ष और मार्गदर्शन

यदि GCV के जरिए स्वचालित lam चयन समान किस्म की टाइम सीरीज़ पर असंगत नतीजे दे रहा है, तो सबसे पहले x के स्पेसिंग पर नज़र डालें। x की यूनिट-स्पेस्ड री-पैरामीटराइज़ेशन पर lam निकालें, फिर उसे औसत स्पेसिंग के घन से गुणा करके मूल पैमाने पर लौटाएँ। यह छोटा सा बदलाव सिस्टम मैट्रिक्स के नॉर्म संतुलन को स्थिर कर देता है और विकृत मिनिमा की संभावना घटाता है। यदि कभी आउटपुट इनपुट से लगभग एक जैसा लगे, तो ध्यान रखें कि आपकी सीरीज़ की संरचना के आधार पर अत्यंत छोटा lam भी GCV का उचित निष्कर्ष हो सकता है। अंत में, यदि फिर भी तेज असंगतियाँ या कई स्थानीय मिनिमा के संकेत मिलें, तो SciPy के इश्यू ट्रैकर में चल रही चर्चाओं और सुधारों पर नज़र रखना उपयोगी है।

इन समायोजनों के साथ, आप make_smoothing_spline से lam को स्वतः चुनते हुए भी भिन्न सैंपलिंग दरों और स्केल वाली डेटासेट्स में संख्यात्मक व्यवहार बेहतर कर सकते हैं।