Gaussian process — kernel দিয়ে অনিশ্চয়তাসহ অনুমান

পর্ব ৩০-এ kernel ridge regression-এর অনুমান k(x_*, X)(K + λI)⁻¹y দেখেছেন, আর সেখানে বলা হয়েছিল এটি একটি Gaussian process-এর গড়ও। এই পর্বে সেই কথাটি খুলে দেখা হবে। লাভ একটাই, কিন্তু বড়: প্রতিটি অনুমানের সঙ্গে একটি সৎ অনিশ্চয়তা, যা ডেটার কাছে কম আর ডেটা থেকে দূরে বেশি। ভিত্তি তিনটি: পর্ব ১৯-এর conditional Gaussian, পর্ব ২৯-এর kernel, আর পর্ব ১০-এর Cholesky. Bayesian ভাবনার মূল কথা (prior, likelihood, posterior) AI-এর ভিত্তি কোর্সের পর্ব ৫০-এ আছে। পর্ব ৪৫ (Bayesian linear regression — posterior আর predictive অনিশ্চয়তা) weight-এর ওপর prior রাখে; এখানে prior বসবে সরাসরি function-এর ওপর। চারটি চালানো যায় এমন ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy, scikit-learn আর matplotlib); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · Function-এর ওপর বণ্টন — Gaussian process কী

পর্ব ১৯-এ multivariate Gaussian দেখেছেন: কয়েকটি সংখ্যার যৌথ বণ্টন, যা একটি গড় vector আর একটি covariance matrix দিয়ে পুরো বর্ণনা করা যায়। এখন সংখ্যাগুলো হোক একটি অজানা function f-এর মান: f(x₁), f(x₂), …, f(x_n). একটি Gaussian process হলো এমন একটি নিয়ম, যা যেকোনো সসীম সংখ্যক x-এর জন্য এই মানগুলোর একটি multivariate Gaussian দেয়। লেখা হয়:

f ~ GP(m(x), k(x, x′))

এখানে m(x) হলো গড় function (প্রায়ই ০ বা ডেটার গড়), আর k(x, x′) হলো covariance function, অর্থাৎ kernel: Cov(f(x), f(x′)) = k(x, x′). n টি বিন্দু বেছে নিলে তাদের covariance matrix K-এর (i, j) ঘর k(x_i, x_j), যা পর্ব ২৯-এর Gram matrix. Kernel বৈধ (positive semidefinite) বলে যেকোনো বাছাইয়ে K একটি বৈধ covariance. তাই GP-কে ভাবা যায় “অসীম মাত্রার Gaussian” হিসেবে, যার প্রতিটি সসীম টুকরো একটি সাধারণ Gaussian.

এই পর্বের প্রধান kernel হলো RBF kernel (squared exponential নামেও পরিচিত):

k(x, x′) = σ_f² · exp(−(x − x′)² / (২ℓ²))

দুটি সংখ্যা এর আচরণ ঠিক করে। σ_f² (signal variance) বলে f গড় থেকে সাধারণত কতটা সরে: প্রতিটি বিন্দুতে f-এর SD σ_f. আর length-scale ℓ বলে কত দূর পর্যন্ত দুটি বিন্দুর মান একসঙ্গে চলে। x আর x′-এর দূরত্ব ℓ হলে correlation e^−০.৫ ≈ ০.৬১, দূরত্ব ২ℓ হলে ০.১৪, ৩ℓ হলে ০.০১১. পর্ব ৩০-এর γ আর এই ℓ একই জিনিস, শুধু ভিন্ন রূপে: γ = ১/(২ℓ²). RBF একটি stationary kernel: মান নির্ভর করে শুধু x − x′-এর ওপর, x কোথায় তার ওপর নয়।

২ · Prior থেকে function টানা — length-scale, Cholesky আর jitter

একটি prior বুঝতে সবচেয়ে ভালো উপায় তা থেকে কয়েকটি নমুনা টেনে দেখা। ধাপগুলো সরল। ০ থেকে ১০ km-এর মধ্যে ২০০টি x-এর একটি grid নিন। তাদের ২০০ × ২০০ covariance matrix K বানান। তারপর পর্ব ১০-এর কৌশল: K = LLᵀ (Cholesky), আর z যদি ২০০টি স্বাধীন standard normal সংখ্যা হয়, তাহলে f = Lz-এর covariance LLᵀ = K. প্রতিটি z থেকে একটি পুরো function পাওয়া যায়।

একটি সমস্যা আছে। ঘন grid-এ পাশাপাশি বিন্দুর মান প্রায় এক, তাই K-এর অনেক eigenvalue শূন্যের অতি কাছে, আর floating-point-এর গোলমালে কিছু সামান্য ঋণাত্মকও হয়ে যায়। তখন Cholesky ব্যর্থ হয়। প্রতিকার jitter: কর্ণে খুব ছোট একটি সংখ্যা যোগ করা, যেমন ১০^−১০. এটি function-কে প্রায় বদলায় না, কিন্তু matrix-কে সংখ্যাগতভাবে positive definite রাখে। নিচের ব্লক তিনটি length-scale-এ একই z দিয়ে চারটি করে function টানে, প্রতিটি K-এর ক্ষুদ্রতম eigenvalue ছাপায়, আর দেখে সবচেয়ে ছোট কোন jitter-এ Cholesky চলে। এখানে scipy.linalg.cholesky নেওয়া হয়েছে, কারণ ব্রাউজারের NumPy-তে np.linalg.cholesky ব্যর্থ হলে error না দিয়ে চুপচাপ NaN ফেরায়। ছবির লেখা ইংরেজিতে, কারণ পাতার matplotlib বাংলা অক্ষর ঠিকমতো আঁকে না।

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cholesky, eigvalsh, LinAlgError

def rbf(a, b, ell, sf=1.0):          # sf^2 exp(-(a-b)^2 / (2 ell^2))
    return sf**2 * np.exp(-(a[:, None] - b[None, :]) ** 2 / (2 * ell**2))

x = np.linspace(0, 10, 200)          # 200টি বিন্দুর grid (km)
rng = np.random.default_rng(46)
z = rng.standard_normal((200, 4))    # একই z, তিনটি length-scale-এ
fig, ax = plt.subplots(1, 3, figsize=(10, 3), sharey=True)
fig.subplots_adjust(bottom=0.18, wspace=0.08)
print("ell  corr(1 km)  crossings  min eig     jitter")
for i, ell in enumerate([0.3, 1.0, 3.0]):
    K = rbf(x, x, ell)
    for jit in [0.0, 1e-10, 1e-8, 1e-6]:   # সবচেয়ে ছোট কোনটিতে চলে
        try:
            L = cholesky(K + jit * np.eye(200), lower=True)
            break
        except LinAlgError:
            continue
    f = L @ z                        # f ~ N(0, K): চারটি function
    cross = np.mean(np.sum(np.diff(np.sign(f), axis=0) != 0, axis=0))
    print(f"{ell:3.1f} {np.exp(-1 / (2 * ell**2)):10.4f} {cross:9.1f}",
          f"{eigvalsh(K)[0]:10.1e} {jit:8.0e}")
    ax[i].plot(x, f, lw=1.2)
    ax[i].set_title(f"length-scale = {ell}")
    ax[i].set_xlabel("x (km)")
ax[0].set_ylabel("f(x) sampled from the prior")
plt.show()

ell  corr(1 km)  crossings  min eig     jitter
0.3     0.0039      11.2   -4.4e-15    1e-10
1.0     0.6065       3.8   -1.7e-14    1e-10
3.0     0.9460       1.0   -6.6e-14    1e-10

৩ · Posterior — conditional Gaussian থেকে GP regression

এবার ডেটা আসে। ধরুন n টি অর্ডারের দূরত্ব X আর ডেলিভারির সময় y জানা, আর প্রতিটি মাপে স্বাধীন Gaussian noise আছে, যার variance σ_n²: y_i = f(x_i) + ε_i. নতুন একটি দূরত্ব x_*-এ f_* = f(x_*) জানতে চাই। GP-র সংজ্ঞা বলে y আর f_* একসঙ্গে Gaussian (গড় ০ ধরে):

[y, f_*] ~ N(০, [[K + σ_n²I, k_*], [k_*ᵀ, k_**]])

এখানে K হলো training বিন্দুগুলোর n × n kernel matrix, k_* হলো x_*-এর সঙ্গে প্রতিটি training বিন্দুর kernel মানের vector, আর k_** = k(x_*, x_*). পর্ব ১৯-এর conditional distribution-এর সূত্র (μ_a + Σ_abΣ_bb⁻¹(x_b − μ_b), আর Σ_aa − Σ_abΣ_bb⁻¹Σ_ba) সরাসরি বসালে কোনো নতুন optimization ছাড়াই posterior পাওয়া যায়:

গড় μ_* = k_*ᵀ(K + σ_n²I)⁻¹y
variance σ_*² = k_** − k_*ᵀ(K + σ_n²I)⁻¹k_*

দুটি সূত্র পড়তে শিখুন। গড়টি training-এর y-গুলোর একটি ওজনদার যোগফল, আর ওজন আসে kernel থেকে: কাছের অর্ডারের কথা বেশি শোনা হয়। Variance-এ k_** হলো prior-এর অনিশ্চয়তা, আর বিয়োগের অংশটি হলো “ডেটা দেখে যতটা অনিশ্চয়তা কমল”। x_* সব training বিন্দু থেকে দূরে হলে k_* ≈ ০, তাই গড় ≈ ০ (prior-এর গড়) আর variance ≈ k_** (prior-এর variance-এর সমান)। আর একটি খুঁটিনাটি: σ_*² হলো f_*-এর variance; একটি নতুন পর্যবেক্ষণ y_*-এর variance তার সঙ্গে σ_n² যোগ করে। প্রথমটিকে বলে epistemic uncertainty (model-এর না-জানা), দ্বিতীয়টিকে aleatoric uncertainty (ডেটার নিজস্ব এলোমেলোতা)। নিচের stepper দুটি training বিন্দুতে পুরো হিসাব হাতে করে।

পুরো ছবিটা এমন দাঁড়ায়: ডেটার কাছে পরিসর সরু, দুটি বিন্দুর মাঝের বড় ফাঁকে চওড়া, আর ডেটার বাইরে prior-এর পরিসরে মিলিয়ে যায়। নিচের ছবিটি পাঁচটি বিন্দুর একটি সত্যিকারের GP posterior থেকে আঁকা।

শেষ বাক্যটি একটি গুরুত্বপূর্ণ সীমা। σ_*²-এর সূত্রে y নেই। অর্থাৎ hyperparameter স্থির থাকলে GP-র অনিশ্চয়তা শুধু “এখানে ডেটা কত ঘন” বলে, “এখানে ডেটা কত এলোমেলো” নয়। দামি এলাকায় ডেলিভারির সময় বেশি ওঠানামা করলেও সাধারণ GP সব জায়গায় এক σ_n ধরে। এখন length-scale-এর প্রভাব দেখুন: দুটি বিন্দু x = ০ আর ২, দুটিরই y = ১, σ_f = ১, σ_n² = ০.০১; প্রশ্ন মাঝের x = ১-এ আর দূরের x = ৫-এ অনুমান কী।

৪ · নিজে বানানো GP — Cholesky, scikit-learn আর kernel ridge

বাস্তবে (K + σ_n²I)⁻¹ কখনো সরাসরি বানানো হয় না। Rasmussen আর Williams-এর বই Gaussian Processes for Machine Learning (২০০৬)-এর Algorithm 2.1 যে পথ দেখায়, তা এই: একবার Cholesky, K + σ_n²I = LLᵀ (প্রায় n³/৬ গুণ, আর সমান সংখ্যক যোগ); তারপর দুটি ত্রিভুজাকার সমাধানে α = (K + σ_n²I)⁻¹y (প্রতিটি প্রায় n² কাজ)। Variance-এর জন্য প্রতিটি test বিন্দুতে আরও একটি ত্রিভুজাকার সমাধান। SciPy-তে cho_factor আর cho_solve ঠিক এই কাজ করে। scikit-learn-এর GaussianProcessRegressor-ও ভেতরে একই algorithm চালায় (তার source code-এ এই Algorithm 2.1-এর উল্লেখ আছে)।

নিচের module দুই ধরনের কাল্পনিক বাজারঘর ডেটা seed দিয়ে বানায়: hub থেকে দূরত্ব বনাম ডেলিভারির মিনিট, আর দিনে অর্ডারের সংখ্যা। আসল ধারা model জানে না।

# বাজারঘর (কাল্পনিক দোকান): দুই ধরনের seeded ডেটা
import numpy as np

def eta_curve(km):                       # আসল ধারা (model জানে না)
    return 25 + 18 * np.log1p(km) + 4 * np.sin(km / 1.3)

def deliveries(n=7, seed=46):            # hub থেকে km -> ডেলিভারির মিনিট
    rng = np.random.default_rng(seed)
    km = np.sort(rng.uniform(0.3, 7, n)).round(2)
    mins = eta_curve(km) + rng.normal(0, 2, n)
    return km.reshape(-1, 1), mins.round(1)

def daily_orders(days=70, seed=46):      # দিনে অর্ডারের সংখ্যা, সাপ্তাহিক ছাঁদ
    rng = np.random.default_rng(seed)
    t = np.arange(days)
    week = np.array([0, 4, 6, 3, 5, 10, -22])   # দিন 0 = শনিবার; বৃহস্পতি বেশি, শুক্র কম
    slow = 10 * np.sin(2 * np.pi * t / 45)        # ধীর ওঠানামা
    y = 140 + slow + week[t % 7] + rng.normal(0, 4, days)
    return t.reshape(-1, 1).astype(float), y.round()

দ্বিতীয় ব্লক মাত্র ৭টি অর্ডারে GP চালায়, hyperparameter হাতে বেঁধে: σ_f = ১০ মিনিট, ℓ = ২ km, σ_n = ২ মিনিট। Prior-এর গড় হিসেবে training-এর y-এর গড় নেওয়া হয়েছে (y থেকে বাদ দিয়ে, শেষে ফেরত যোগ)। তারপর একই হিসাব GaussianProcessRegressor(optimizer=None)-এ, যাতে scikit-learn hyperparameter না বদলায়; এখানে alpha মানে কর্ণে যোগ করা noise variance σ_n². আর KernelRidge-এ: তার RBF-এ σ_f² নেই, তাই λ = σ_n²/σ_f² আর γ = ১/(২ℓ²) দিলে একই গড় আসার কথা।

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cho_factor, cho_solve
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel as C, RBF
from sklearn.kernel_ridge import KernelRidge
from bazar_gp import deliveries

X, y = deliveries()                  # 7টি অর্ডার: km -> মিনিট
ym = y.mean(); yc = y - ym           # prior-এর গড় = training গড়
sf, ell, sn = 10.0, 2.0, 2.0         # amplitude, length-scale, noise sd
k = lambda A, B: sf**2 * np.exp(-(A - B.T) ** 2 / (2 * ell**2))
xs = np.linspace(0, 12, 241).reshape(-1, 1)

c = cho_factor(k(X, X) + sn**2 * np.eye(len(y)))   # K + sn^2 I = L L^T
alpha = cho_solve(c, yc)
Ks = k(xs, X)
mean = ym + Ks @ alpha               # posterior গড়
var = sf**2 - np.sum(Ks * cho_solve(c, Ks.T).T, axis=1)
sd = np.sqrt(np.maximum(var, 0))
gp = GaussianProcessRegressor(C(sf**2) * RBF(ell), alpha=sn**2,
                              optimizer=None).fit(X, yc)
m2, s2 = gp.predict(xs, return_std=True)
krr = KernelRidge(alpha=sn**2 / sf**2, kernel="rbf",
                  gamma=1 / (2 * ell**2)).fit(X, yc)
print("max |mean - GP|  =", f"{np.max(np.abs(mean - ym - m2)):.1e}")
print("max |sd - GP|    =", f"{np.max(np.abs(sd - s2)):.1e}")
print("max |mean - KRR| =", f"{np.max(np.abs(mean - ym - krr.predict(xs))):.1e}")
for x0 in [0.85, 3.3, 9.0, 12.0]:
    i = int(round(x0 / 0.05))
    print(f"x = {x0:5.2f} km: mean {mean[i]:5.1f}, 95% +/- {1.96 * sd[i]:4.1f}")
plt.fill_between(xs[:, 0], mean - 1.96*sd, mean + 1.96*sd, alpha=0.25, label="95% band")
plt.plot(xs, mean, label="posterior mean"); plt.plot(X, y, "ko", label="orders")
plt.axhline(ym, ls=":", c="gray", label="prior mean"); plt.legend()
plt.xlabel("distance from hub (km)"); plt.ylabel("delivery time (min)"); plt.show()

max |mean - GP|  = 2.3e-14
max |sd - GP|    = 1.4e-14
max |mean - KRR| = 2.1e-14
x =  0.85 km: mean  38.8, 95% +/-  2.7
x =  3.30 km: mean  54.6, 95% +/-  4.2
x =  9.00 km: mean  53.2, 95% +/- 16.6
x = 12.00 km: mean  49.7, 95% +/- 19.6

৫ · Hyperparameter শেখা — log marginal likelihood

σ_f, ℓ আর σ_n কীভাবে বাছবেন? Cross-validation চলে, কিন্তু GP একটি সম্ভাবনার model বলে আরেকটি স্বাভাবিক পথ আছে: hyperparameter θ দিলে ডেটা y কতটা সম্ভাব্য, সেটি মাপা। f-কে “integrate out” করলে, অর্থাৎ সব সম্ভাব্য function-এর ওপর গড় নিলে, y নিজেই একটি Gaussian: y ~ N(০, K_y), যেখানে K_y = K + σ_n²I. এর log density হলো log marginal likelihood:

log p(y | X, θ) = −½ yᵀK_y⁻¹y − ½ log|K_y| − (n/২) log ২π

প্রথম পদ data-fit: ডেটা kernel-এর বলা “আকারের” সঙ্গে কতটা মেলে। দ্বিতীয় পদ complexity-র শাস্তি: |K_y| (determinant) বড় মানে prior অনেক ধরনের ডেটা সম্ভব বলে, আর তখন যেকোনো একটি নির্দিষ্ট ডেটার সম্ভাবনা কম। ছোট ℓ-এর model প্রায় সবকিছু ব্যাখ্যা করতে পারে, তাই প্রতিটির ভাগে সম্ভাবনা কম পড়ে; বড় ℓ-এর model কম ধরনের ডেটা সম্ভব বলে, কিন্তু ডেটা না মিললে data-fit-এ বড় শাস্তি খায়। এই টানাপোড়েন একটি স্বয়ংক্রিয় Occam's razor, যা David MacKay ১৯৯২ সালের “Bayesian interpolation” প্রবন্ধে স্পষ্ট করে দেখান। Cholesky থাকলে দুটিই সস্তা: yᵀK_y⁻¹y = yᵀα, আর log|K_y| = ২ Σ log L_ii (log-determinant, overflow ছাড়া)। নিচের stepper তিনটি বিন্দুতে দুটি পদ হাতে মেলায়।

θ-এর ওপর এই পরিমাণ সর্বোচ্চ করাকে বলে empirical Bayes বা type-II maximum likelihood: f-এর ওপর পুরো Bayesian গড়, কিন্তু hyperparameter-এর জন্য একটি বিন্দু-অনুমান। scikit-learn এর gradient বন্ধ রূপে গুনে পর্ব ১৫-এর L-BFGS চালায়, θ রাখে লগ-স্কেলে (kernel_.theta), আর প্রতিটি hyperparameter-এর একটি সীমা (bounds) দিতে হয়। Likelihood-এর তলে একাধিক চূড়া থাকতে পারে, তাই n_restarts_optimizer দিয়ে কয়েকটি এলোমেলো শুরু থেকে চালানো ভালো। তৃতীয় ব্লক ৩০টি অর্ডারে আগে ℓ-এর grid-এ দুটি পদ নিজে গোনে, তারপর ConstantKernel * RBF + WhiteKernel শেখায়। WhiteKernel কর্ণে একটি শেখার যোগ্য noise variance যোগ করে।

import numpy as np
from scipy.linalg import cho_factor, cho_solve
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel as C
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
from bazar_gp import deliveries

X, y = deliveries(n=30, seed=47)     # 30টি অর্ডার
yc = y - y.mean(); n = len(y)

def lml_parts(sf, ell, sn):          # log p(y | X, theta) দুই ভাগে
    K = sf**2 * np.exp(-(X - X.T) ** 2 / (2 * ell**2)) + sn**2 * np.eye(n)
    c = cho_factor(K, lower=True)
    fit = -0.5 * yc @ cho_solve(c, yc)       # data-fit: -1/2 y^T K^-1 y
    comp = -np.sum(np.log(np.diag(c[0])))    # complexity: -1/2 log|K|
    return fit, comp, fit + comp - n / 2 * np.log(2 * np.pi)

print("  ell  data-fit  complexity    total")
for ell in [0.2, 0.5, 1, 2, 4, 8]:
    f, cpx, tot = lml_parts(10.0, ell, 2.0)
    print(f"{ell:5}  {f:8.1f}  {cpx:10.1f}  {tot:7.1f}")

kernel = (C(100.0, (1.0, 1e4)) * RBF(2.0, (0.1, 50.0))
          + WhiteKernel(1.0, (1e-2, 1e2)))
gp = GaussianProcessRegressor(kernel, n_restarts_optimizer=3,
                              random_state=46).fit(X, yc)
p = np.exp(gp.kernel_.theta)         # theta লগ-স্কেলে রাখা
sf, ell, sn = np.sqrt(p[0]), p[1], np.sqrt(p[2])
print("\nlearned:", gp.kernel_)
print(f"sigma_f = {sf:.2f}, ell = {ell:.2f} km, noise sd = {sn:.2f} min")
print(f"sklearn LML = {gp.log_marginal_likelihood_value_:.3f}",
      f"| scratch = {lml_parts(sf, ell, sn)[2]:.3f}")

  ell  data-fit  complexity    total
  0.2      -5.9       -53.3    -86.8
  0.5      -7.9       -42.5    -78.0
    1      -8.6       -35.4    -71.6
    2     -10.0       -30.6    -68.1
    4     -15.4       -27.7    -70.7
    8     -25.9       -26.2    -79.6

learned: 11**2 * RBF(length_scale=2.45) + WhiteKernel(noise_level=2.57)
sigma_f = 10.96, ell = 2.45 km, noise sd = 1.60 min
sklearn LML = -66.876 | scratch = -66.876

৬ · Kernel জোড়া — সাপ্তাহিক অর্ডার আর ডেটার বাইরের অনুমান

Kernel-এর একটি সুন্দর গুণ: দুটি বৈধ kernel-এর যোগফল আর গুণফলও বৈধ kernel. যোগ মানে “f দুটি স্বাধীন অংশের যোগ”, গুণ মানে “দুটি শর্তই একসঙ্গে”। এই kernel composition দিয়ে ডেটার জানা গঠন model-এ বলা যায়। scikit-learn-এ কয়েকটি মূল উপাদান:

বাজারঘরের দৈনিক অর্ডারে তিনটি অংশ দেখা যায়: সপ্তাহের ছাঁদ (বৃহস্পতিবার বেশি, শুক্রবার অনেক কম), কয়েক সপ্তাহের ধীর ওঠানামা, আর প্রতিদিনের এলোমেলোতা। তাই kernel হবে periodic (p = ৭ দিন, স্থির) + RBF (ধীর অংশ) + WhiteKernel (noise), প্রতিটির সামনে একটি মাপ। normalize_y=True দিলে scikit-learn y-কে গড় ০ আর SD ১-এ আনে, তাই মাপের সীমা ১-এর আশপাশে দেওয়া হয়েছে। Bounds-এর মাপ ভুল হলে শেখা মান সীমায় আটকে যায় আর scikit-learn একটি ConvergenceWarning দেয়, যা চাপা না দিয়ে সীমা ঠিক করাই উচিত (এই ব্লকের সীমা এভাবেই ঠিক করা হয়েছে)। ১০ সপ্তাহের ডেটা দিয়ে শিখে আরও ৮ সপ্তাহ সামনে অনুমান করা হবে। শেষে একটি আলাদা অংশে n বাড়িয়ে শেখা kernel-এর n × n matrix-এর Cholesky-র সময় মাপা হবে।

import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cho_factor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ExpSineSquared, RBF, WhiteKernel
from sklearn.gaussian_process.kernels import ConstantKernel as C
from bazar_gp import daily_orders

t, y = daily_orders(days=70)         # 10 সপ্তাহের দৈনিক অর্ডার
kernel = (C(1.0, (1e-2, 1e2)) * ExpSineSquared(1.0, 7.0, (0.01, 10), "fixed")
          + C(1.0, (1e-2, 1e2)) * RBF(20.0, (5, 200))
          + WhiteKernel(0.1, (1e-3, 1)))      # y scale করা, তাই মাপ ~1
gp = GaussianProcessRegressor(kernel, normalize_y=True, n_restarts_optimizer=2,
                              random_state=46).fit(t, y)
print("learned:", gp.kernel_)
tf = np.arange(0, 126).reshape(-1, 1).astype(float)   # আরও 8 সপ্তাহ
m, s = gp.predict(tf, return_std=True)
wd = "Sat Sun Mon Tue Wed Thu Fri".split()      # দিন 0 = শনিবার
for d in [68, 69, 75, 76, 124, 125]:
    print(f"day {d:3d} ({wd[d % 7]}): {m[d]:6.1f} +/- {1.96 * s[d]:4.1f}")
plt.figure(figsize=(9, 3.4)); plt.subplots_adjust(bottom=0.16)
plt.fill_between(tf[:, 0], m - 1.96*s, m + 1.96*s, alpha=0.25, label="95% band")
plt.plot(tf, m, lw=1, label="GP mean"); plt.plot(t, y, "k.", label="orders")
plt.axvline(69.5, ls="--", c="gray"); plt.axhline(y.mean(), ls=":", c="gray")
plt.xlabel("day"); plt.ylabel("orders per day"); plt.legend(fontsize=8); plt.show()

print("\n    n  Cholesky ms")        # O(n^3): n দ্বিগুণ হলে প্রায় 8 গুণ
rng = np.random.default_rng(0)
for n in [250, 500, 1000, 2000]:
    xx = np.sort(rng.uniform(0, 100, n)).reshape(-1, 1)
    K = gp.kernel_(xx)               # শেখা kernel-এর n x n matrix
    t0 = time.perf_counter(); cho_factor(K)
    print(f"{n:5d} {(time.perf_counter() - t0) * 1e3:11.1f}")

learned: 0.826**2 * ExpSineSquared(length_scale=0.0512, periodicity=7) + 0.49**2 * RBF(length_scale=7.06) + WhiteKernel(noise_level=0.113)
day  68 (Thu):  150.2 +/-  8.8
day  69 (Fri):  118.3 +/-  9.0
day  75 (Thu):  150.4 +/- 12.1
day  76 (Fri):  119.4 +/- 12.6
day 124 (Thu):  152.4 +/- 14.6
day 125 (Fri):  121.0 +/- 14.6

    n  Cholesky ms
  250         0.9
  500         5.0
 1000        36.2
 2000       265.6

৭ · দাম, বিস্তার আর সীমা — O(n³), sparse GP আর classification

দাম। Exact GP-তে n × n matrix রাখতে মেমরি O(n²), Cholesky-তে সময় O(n³) (time complexity), আর hyperparameter শেখার প্রতিটি ধাপে আবার Cholesky. Fit হয়ে গেলে প্রতিটি নতুন বিন্দুর গড় O(n), variance O(n²). তাই কয়েক হাজার বিন্দু পর্যন্ত exact GP আরামে চলে, কয়েক লাখে নয়।

Sparse GP. প্রধান প্রতিকার inducing points: m টি (m ≪ n) প্রতিনিধি বিন্দু Z, আর পুরো function-এর তথ্য তাদের মধ্য দিয়ে চালানো। তখন খরচ নামে প্রায় O(nm²)-এ। এটি পর্ব ৩১-এর Nyström ধারণার আত্মীয়: K ≈ K_nmK_mm⁻¹K_mn. কয়েকটি মাইলফলক: Snelson আর Ghahramani-র FITC (NIPS ২০০৫), Quiñonero-Candela আর Rasmussen-এর একীভূত দৃষ্টি (JMLR ২০০৫), Titsias-এর variational inducing points (AISTATS ২০০৯), আর Hensman, Fusi ও Lawrence-এর stochastic রূপ (UAI ২০১৩), যা mini-batch-এ লাখো বিন্দুতে চলে। এসব sparse GP scikit-learn-এ নেই; GPyTorch বা GPflow-এর মতো library-তে পাওয়া যায়। Variational অংশটি পর্ব ৪৮ (Variational inference-এর ধারণা — ELBO আর approximation)-এর ভাষায় লেখা।

Classification. Label ০/১ হলে likelihood আর Gaussian নয় (যেমন logistic), তাই posterior-এর বন্ধ রূপ থাকে না। scikit-learn-এর GaussianProcessClassifier একটি latent function-এর posterior-কে তার শীর্ষবিন্দুর চারপাশে একটি Gaussian দিয়ে আনুমানিক করে, যাকে বলে Laplace approximation (GPML বইয়ের Algorithm 3.1, 3.2 আর 5.1; ধারণাটি Williams আর Barber-এর ১৯৯৮ সালের প্রবন্ধের)। বহু শ্রেণিতে এটি one-vs-rest বা one-vs-one চালায়। বাজারঘরে “এই অর্ডার ফেরত আসবে কি না”-র মতো প্রশ্নে এটি অনিশ্চয়তাসহ সম্ভাবনা দেয়, তবে পর্ব ২৬-এর মতো calibration যাচাই করতেই হবে।

সীমা। (১) Stationary kernel সব জায়গায় একই মসৃণতা আর একই noise ধরে। (২) RBF খুব মসৃণ; Michael Stein তাঁর ১৯৯৯ সালের বই Interpolation of Spatial Data-তে বাস্তব ডেটার জন্য Matérn kernel পছন্দ করার যুক্তি দেন। (৩) Extrapolation শুধু kernel যা বলে তা-ই: RBF বলে “গড়ে ফেরো”, periodic বলে “চিরকাল একই ছাঁদ”; কোনোটিই বাস্তব জগতের নিশ্চয়তা নয়। (৪) পরিসর তখনই বিশ্বাসযোগ্য, যখন model ঠিক; আলাদা ডেটায় coverage মেপে দেখুন, আর model-নির্ভরতা ছাড়া পরিসর চাইলে পর্ব ৫৮ (Conformal prediction — নিশ্চয়তাসহ অনুমানের পরিসর) দেখবেন।

৮ · AI-এর যুগে — GP কোড যাচাই

“ডেলিভারির সময়ের জন্য Gaussian process বানাও, অনিশ্চয়তাসহ” বললে AI assistant দ্রুত কোড দেয়, আর ছবিতে সুন্দর নীল ছায়াও আসে। কিন্তু এই বিষয়ের ভুল প্রায়ই ছবিতে ধরা পড়ে না।

AI-এর দেওয়া GP কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• দুটি বিন্দুর একটি খেলনা উদাহরণে গড় আর variance হাতে বা নিজের Cholesky কোডে বের করে scikit-learn-এর সঙ্গে মিলিয়েছি
• gp.kernel_ ছাপিয়ে দেখেছি কোনো hyperparameter সীমায় আটকে নেই, আর warning চাপা দিইনি
• noise একবারই গোনা হয়েছে, আর পরিসর f-এর না নতুন পর্যবেক্ষণের, তা জানি
• সময়ের ডেটায় ভবিষ্যতের অংশ তালাবদ্ধ রেখে backtest করেছি, আর পরিসরের coverage মেপেছি
• একটি সরল baseline (গড়, seasonal naive, linear model) আর কয়েকটি seed-এর সঙ্গে তুলনা করেছি
• কোনো আসল গ্রাহকের ঠিকানা, ফোন নম্বর বা বিক্রির গোপন সংখ্যা AI-তে পেস্ট করিনি