Regularization-এর তত্ত্ব — L1/L2-এর জ্যামিতি, MAP আর early stopping

AI-এর ভিত্তি কোর্সের পর্ব ৬১-এ ridge নিজের হাতে লিখে scikit-learn-এর সঙ্গে মিলিয়েছেন, lasso-কে coefficient ঠিক শূন্য করতে দেখেছেন, আর L1 হীরা বনাম L2 বৃত্তের ছবি আর “regularization মানে prior” কথাটি এক অনুচ্ছেদে পেয়েছেন। এই পর্ব সেখান থেকে এক ধাপ ভেতরে যায়: ঠিক কোন দিকে কতটা ছোট হয়, শূন্য কোন সমীকরণ থেকে আসে, আর শাস্তি ছাড়াই (শুধু আগে থামিয়ে) একই ফল কেন আসে? হাতিয়ার আগের মডিউলের: পর্ব ০৯-এর SVD আর পর্ব ১১-এর gradient. চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় দশ সেকেন্ড।

১ · কেন regularize — ill-posed সমস্যা আর variance

গণিতবিদ Jacques Hadamard বিশ শতকের শুরুতে একটি সমস্যাকে well-posed বলেছিলেন যদি তিনটি শর্ত মানে: সমাধান আছে, সমাধান একটিই, আর ডেটা সামান্য বদলালে সমাধানও সামান্য বদলায়। ML-এর least squares প্রায়ই দ্বিতীয় বা তৃতীয় শর্তে আটকায়। Feature-এর সংখ্যা p যদি সারির সংখ্যা n-এর চেয়ে বেশি হয়, অসংখ্য w একই ভুল দেয় (একক সমাধান নেই)। আর দুটি feature প্রায় একই হলে (multicollinearity) সমাধান একটি থাকলেও তা ডেটার noise-এ লাফায়। পর্ব ০৯-এর ভাষায় এর উৎস ছোট singular value.

সংখ্যায় দেখা যাক। মডেল y = Xw + ε, যেখানে ε-এর প্রতিটি ঘর স্বাধীন, গড় ০ আর variance σ_ε². Least squares-এর উত্তর ŵ = (XᵀX)⁻¹Xᵀy, আর তার covariance:

Cov(ŵ) = σ_ε² (XᵀX)⁻¹ = σ_ε² ∑_j v_jv_jᵀ / σ_j²

এখানে σ_j হলো X-এর j-তম singular value আর v_j তার ডান singular vector (X = UΣVᵀ). অর্থ: v_j দিকে ŵ-এর variance σ_ε²/σ_j². একটি σ_j যদি ০.০১ হয়, ওই দিকে noise ১০,০০০ গুণ ফুলে ওঠে। এটিই variance, আর overfitting-এর গাণিতিক চেহারা। Regularization এই ফোলা আটকায়, বিনিময়ে একটু bias নেয়। পর্ব ০৩-এর ভাষায় এটি hypothesis-এর ক্লাসকে ছোট রাখা, অর্থাৎ structural risk minimization-এর একটি রূপ।

এই ধারণার শিকড় দুই জায়গায়। পদার্থবিদ্যা আর প্রকৌশলের inverse problem-এ Andrey Tikhonov ১৯৪০ থেকে ১৯৬০-এর দশকে অস্থির সমীকরণকে স্থির করার এই কৌশল গড়ে তোলেন, তাই একে Tikhonov regularization-ও বলে। পরিসংখ্যানে Arthur Hoerl আর Robert Kennard ১৯৭০ সালে Technometrics-এ ridge regression প্রস্তাব করেন, প্রায়-collinear feature-এর অস্থির coefficient সারাতে।

২ · Ridge-কে SVD-র চোখে — প্রতিটি দিকে আলাদা ছাঁটাই

Ridge কমায় ‖y − Xw‖² + λ‖w‖², আর পর্ব ১১-এ দেখেছেন gradient শূন্য বসালে (XᵀX + λI)w = Xᵀy. ধরুন X আগে থেকে center করা (তাই intercept আলাদা), আর তার thin SVD X = UΣVᵀ: U-এর কলাম u_j (দৈর্ঘ্য n), V-এর কলাম v_j (দৈর্ঘ্য p), আর σ₁ ≥ σ₂ ≥ … ≥ ০. তাহলে XᵀX = VΣ²Vᵀ, আর XᵀX + λI = V(Σ² + λI)Vᵀ: λ যোগ করলে eigenvector বদলায় না, শুধু প্রতিটি eigenvalue σ_j² থেকে σ_j² + λ হয়। Inverse নিয়ে গুণ করলে:

w_λ = ∑_j v_j · [σ_j/(σ_j² + λ)] · u_jᵀy    ŷ = ∑_j u_j · [σ_j²/(σ_j² + λ)] · u_jᵀy

λ = ০ হলে বন্ধনীর গুণক ১, আর পাই সাধারণ least squares: প্রতিটি দিকে θ̂_j = u_jᵀy/σ_j. Ridge সেই একই সংখ্যাকে গুণ করে shrink factor f_j = σ_j²/(σ_j² + λ) দিয়ে, যা সবসময় ০ আর ১-এর মধ্যে। σ_j² ≫ λ হলে f_j ≈ ১ (প্রায় অক্ষত); σ_j² ≪ λ হলে f_j ≈ ০ (প্রায় মুছে যায়)। এই কারণেই ridge-কে বলা হয় shrinkage পদ্ধতি, আর পর্ব ০৯-এর truncated pseudo-inverse-এর সঙ্গে তুলনাটিও স্পষ্ট: সেখানে গুণক হয় ১ নয় ০ (কড়া ছাঁটাই), এখানে মসৃণ।

কেন ছোট σ-র দিকটিই কাটা উচিত? কারণ u_j একক দৈর্ঘ্যের, তাই u_jᵀε-এর noise প্রতিটি দিকে সমান, σ_ε. ভাগ করার সময় ছোট σ_j সেই noise-কে ফোলায়। নিচের stepper দুটি দিক নিয়ে হিসাবটি দেখায়।

৩ · Effective degrees of freedom — λ দিয়ে model-এর আকার মাপা

Ridge-এর ŷ = H_λy, যেখানে H_λ = X(XᵀX + λI)⁻¹Xᵀ হলো ridge-এর hat matrix. Least squares-এ H একটি projection, আর তার trace হুবহু p, অর্থাৎ parameter-এর সংখ্যা। Ridge-এর জন্য একই সংজ্ঞা নিলে:

df(λ) = tr(H_λ) = ∑_j σ_j²/(σ_j² + λ)

λ = ০-তে df = p (X-এর rank পূর্ণ হলে), আর λ → ∞-এ df → ০. মাঝে এটি একটি ভগ্নাংশ, যেমন ২.৪৭: “এই model প্রায় আড়াইটি parameter-এর সমান নমনীয়”. সংজ্ঞাটির একটি গভীর রূপও আছে (Bradley Efron-এর কাজে জনপ্রিয়): df = ∑_i Cov(ŷ_i, y_i)/σ_ε², অর্থাৎ প্রতিটি y_i নিজের অনুমান ŷ_i-কে কতটা টানে। Linear smoother-এ এটি trace-এর সমান। Lasso-র ক্ষেত্রে Zou, Hastie আর Tibshirani ২০০৭ সালে দেখান যে কিছু শর্তে অশূন্য coefficient-এর সংখ্যাই df-এর একটি unbiased অনুমান। এই সংখ্যা λ-এর তুলনা সহজ করে (বিভিন্ন ডেটায় λ-এর মান তুলনীয় নয়, df তুলনীয়), আর generalized cross-validation বা C_p-এর মতো সূত্রে বসে।

নিচের ব্লক বাজারঘরের ৪০ জন কাল্পনিক গ্রাহকের ডেটায় তিনটি কাজ করে। “দেখা পাতা” কলামটি ইচ্ছা করে “app-এ ঢোকা”-র প্রায় নকল, যাতে একটি singular value খুব ছোট হয়। (১) SVD সূত্রে ridge বানিয়ে scikit-learn-এর Ridge-এর সঙ্গে মেলায়, (২) df আর hat matrix-এর trace মেলায়, (৩) λ বদলে df-এর বক্ররেখা আঁকে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge
np.set_printoptions(precision=3, suppress=True)

rng = np.random.default_rng(42)          # কাল্পনিক বাজারঘর: ৪০ জন গ্রাহক
n = 40
visits = rng.normal(10, 3, n)                      # মাসে app-এ ঢোকা
pages = 2.5 * visits + rng.normal(0, 0.5, n)       # দেখা পাতা: visits-এর প্রায় নকল
days = rng.integers(1, 60, n).astype(float)        # শেষ অর্ডারের পর দিন
basket = rng.normal(9, 2.5, n)                     # গড় ঝুড়ি, শত টাকা
X = np.column_stack([visits, pages, days, basket])
y = 200 + 30 * visits + 4 * pages - 6 * days + 40 * basket + rng.normal(0, 60, n)
Xs = (X - X.mean(0)) / X.std(0)                    # standardize: সব কলাম একই মাপে
yc = y - y.mean()

U, s, Vt = np.linalg.svd(Xs, full_matrices=False)
lam = 10.0; f = s**2 / (s**2 + lam)               # প্রতিটি SVD দিকের shrink factor
w_svd = Vt.T @ ((s / (s**2 + lam)) * (U.T @ yc))   # = V diag(σ/(σ²+λ)) Uᵀy
w_sk = Ridge(alpha=lam).fit(Xs, y).coef_
print("singular values σ :", s)
print("shrink σ²/(σ²+λ)  :", f)
print("ridge via SVD     :", w_svd)
print("sklearn Ridge     :", w_sk, " max diff", f"{np.abs(w_svd - w_sk).max():.1e}")

lams = np.logspace(-2, 4, 60)
df = [np.sum(s**2 / (s**2 + l)) for l in lams]
H = Xs @ np.linalg.solve(Xs.T @ Xs + lam * np.eye(4), Xs.T)   # hat matrix, λ = 10
print(f"df(10) = {np.sum(f):.3f}, trace(H) = {np.trace(H):.3f}")
for l in (0.1, 10, 100, 1000):
    print(f"λ = {l:<6} df = {np.sum(s**2 / (s**2 + l)):.2f}")
plt.figure(figsize=(6, 3.4), layout="constrained")
plt.semilogx(lams, df, lw=2)
plt.xlabel("lambda (log scale)"); plt.ylabel("effective df")
plt.title("ridge: effective degrees of freedom"); plt.grid(alpha=0.3); plt.show()

singular values σ : [9.187 6.733 5.496 0.254]
shrink σ²/(σ²+λ)  : [0.894 0.819 0.751 0.006]
ridge via SVD     : [ 39.862  42.37  -72.267  84.321]
sklearn Ridge     : [ 39.862  42.37  -72.267  84.321]  max diff 8.0e-13
df(10) = 2.471, trace(H) = 2.471
λ = 0.1    df = 3.39
λ = 10     df = 2.47
λ = 100    df = 1.00
λ = 1000   df = 0.15

৪ · L1, subgradient আর soft-thresholding — lasso-র coordinate update

Lasso-র শাস্তি λ‖w‖₁ = λ∑|w_j|, অর্থাৎ L1 norm. সমস্যা: |w|-এর w = ০-তে কোনো derivative নেই, কোণা আছে। তাই gradient শূন্য বসানোর পুরোনো কৌশল খাটে না। এর সমাধান subgradient: একটি convex function-এর কোনো বিন্দুতে যে সব ঢালের রেখা পুরো function-এর নিচে থাকে, তাদের ঢালের সেট। |w|-এর জন্য:

∂|w| = {+1} যদি w > 0,   {−1} যদি w < 0,   [−1, +1] যদি w = 0

আর convex function-এর সর্বনিম্ন বিন্দুর শর্ত: ০ সেই বিন্দুর subgradient সেটের মধ্যে আছে। এবার সবচেয়ে ছোট lasso: এক চলক, min_w ½(w − z)² + λ|w|. এখানে z হলো শাস্তি ছাড়া যা উত্তর হতো। শর্ত লিখলে: ০ ∈ w − z + λ·∂│w│. তিনটি ক্ষেত্র আলাদা করে দেখলে উত্তর আসে soft-thresholding:

w* = S(z, λ) = sign(z) · max(|z| − λ, 0)

অনেক feature-এ: coordinate descent

Feature অনেক আর পরস্পর সম্পর্কিত হলে সব coefficient একসঙ্গে closed form-এ মেলে না। কিন্তু একটি বাদে বাকি সব স্থির রাখলে সমস্যাটি ঠিক উপরের এক-চলকের রূপ নেয়। scikit-learn-এর Lasso যা কমায় (ডকুমেন্টেশন অনুযায়ী), তা দিয়েই লিখি:

L(w) = (1/2n)‖y − Xw‖² + α‖w‖₁

w_j বাদে বাকিদের residual r_(j) = y − ∑_k≠j x_kw_k. তাহলে w_j-এর অংশ (1/2n)‖r_(j) − x_jw_j‖² + α|w_j|. ρ_j = x_jᵀr_(j)/n আর z_j = x_jᵀx_j/n লিখে বর্গ পূর্ণ করলে, subgradient শর্ত থেকে আসে:

w_j ← S(ρ_j, α) / z_j

Standardize করা কলামে z_j = ১. প্রতিটি j-তে এই update পালাক্রমে চালানো, যতক্ষণ না কোনো coefficient আর নড়ে: এটিই coordinate descent. Lasso-র loss convex, আর শাস্তিটি প্রতিটি coordinate-এ আলাদা হয়ে যায় (separable), তাই এই পদ্ধতি সর্বনিম্নে পৌঁছায় (Paul Tseng-এর ২০০১ সালের ফলাফল এই ধরনের সমস্যার জন্য প্রচলিত তত্ত্ব)। আরেকটি সহজ ফল: w = ০-তে subgradient শর্ত হলো প্রতিটি j-এ |x_jᵀy|/n ≤ α. তাই α_max = max_j |x_jᵀy|/n-এর ওপরে সব coefficient শূন্য; regularization path এখান থেকেই শুরু হয়।

নিচের ব্লক বাজারঘরের ৫০টি কাল্পনিক ডেলিভারিতে (৮টি feature, শুধু প্রথম তিনটি কাজের) এই update নিজে লেখে, scikit-learn-এর সঙ্গে মেলায়, তারপর α_max থেকে নিচে ৪০টি α-তে পুরো path আঁকে। প্রতিটি α-র সমাধান আগের α-র উত্তর থেকে শুরু হয়; একে warm start বলে, আর path দ্রুত হওয়ার মূল কৌশল এটিই।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso
np.set_printoptions(precision=3, suppress=True)

rng = np.random.default_rng(7)           # বাজারঘরের ৫০টি কাল্পনিক ডেলিভারি, ৮টি feature
n, p = 50, 8                              # শুধু প্রথম ৩টি (দূরত্ব, ওজন, যানজট) কাজের
X = rng.normal(size=(n, p)); X = (X - X.mean(0)) / X.std(0)
y = X @ np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0]) + rng.normal(0, 1.5, n); y -= y.mean()

def soft(z, a): return np.sign(z) * np.maximum(np.abs(z) - a, 0.0)   # S(z, a)

def lasso_cd(X, y, a, w=None, tol=1e-10):    # (1/2n)‖y − Xw‖² + a‖w‖₁
    w = np.zeros(p) if w is None else w.copy()
    z, r = (X**2).sum(0) / n, y - X @ w
    for sweep in range(1000):
        w_old = w.copy()
        for j in range(p):
            r += X[:, j] * w[j]                    # j বাদে বাকিদের residual
            w[j] = soft(X[:, j] @ r / n, a) / z[j]
            r -= X[:, j] * w[j]
        if np.abs(w - w_old).max() < tol: return w
    return w

w_cd = lasso_cd(X, y, 0.5)
w_sk = Lasso(alpha=0.5, tol=1e-12, max_iter=100000).fit(X, y).coef_
print("CD     :", w_cd + 0.0); print("sklearn:", w_sk + 0.0, " max diff", f"{np.abs(w_cd - w_sk).max():.1e}")
a_max = np.abs(X.T @ y).max() / n
alphas = a_max * np.logspace(0.01, -3, 40); w, path = None, []
for a in alphas:
    w = lasso_cd(X, y, a, w); path.append(w)       # warm start: আগের উত্তর থেকে শুরু
path = np.array(path)
print(f"alpha_max = {a_max:.3f}; non-zeros on path:", [int((r != 0).sum()) for r in path[::5]])
plt.figure(figsize=(6, 3.4), layout="constrained"); plt.semilogx(alphas, path); plt.gca().invert_xaxis()
plt.xlabel("alpha (log scale, decreasing)"); plt.ylabel("coefficient"); plt.title("lasso path by coordinate descent"); plt.show()

CD     : [ 2.132 -1.532  1.043  0.     0.     0.     0.     0.   ]
sklearn: [ 2.132 -1.532  1.043  0.     0.     0.     0.     0.   ]  max diff 1.9e-12
alpha_max = 2.515; non-zeros on path: [0, 3, 3, 4, 6, 6, 7, 7]

৫ · জ্যামিতি — শাস্তির রূপ বনাম সীমার রূপ, আর L1 বলের কোণ

একই lasso দুইভাবে লেখা যায়। শাস্তির রূপ: min ‖y − Xw‖² + λ‖w‖₁. সীমার রূপ: min ‖y − Xw‖², শর্ত ‖w‖₁ ≤ t. প্রতিটি λ > ০-এর জন্য একটি t আছে (সেই λ-এর সমাধানের ‖w‖₁ নিজেই), যেখানে দুই সমস্যার উত্তর এক; ridge-এর জন্যও একই কথা। সংযোগটি হলো Lagrange multiplier: সীমার রূপের KKT শর্ত লিখলে যে multiplier আসে, সেটিই λ. এই শর্তগুলোর পুরো তত্ত্ব পর্ব ১৬ (Constraint-সহ optimization — Lagrange multiplier আর KKT)-এর বিষয়; এখানে শুধু ফল লাগবে। দুটি সতর্কতা: t আর λ-এর সম্পর্ক ডেটার ওপর নির্ভর করে, আর lasso-য় একই সমাধান অনেক সময় এক সারি λ-এর জন্য থাকে না, কারণ path টুকরো টুকরো সরলরেখা।

সীমার রূপে জ্যামিতি পরিষ্কার। AI-এর ভিত্তি কোর্সের পর্ব ৬১-এর হীরা আর বৃত্তের ছবিটি মনে করুন। এবার সেই ছবির পেছনের হিসাবটি সংখ্যায়, সবচেয়ে সরল ক্ষেত্রে: X-এর কলাম orthonormal, তাই loss-এর সমান-মানের রেখা বৃত্তাকার আর কেন্দ্র ŵ.

তাহলে L0 (অশূন্য coefficient গোনা) সরাসরি শাস্তি দিলে কী হয়? সেটি best subset selection, যা সবচেয়ে sparse, কিন্তু non-convex আর সাধারণভাবে combinatorial: p = ৪০ হলে প্রায় ১.১ লাখ কোটি উপসেট। |w|^q (q < ১)-ও বেশি sparse, কিন্তু non-convex. L1 হলো সবচেয়ে “sparse-প্রেমী” convex শাস্তি: q = ১-এর নিচে নামলেই convexity হারায়। কেন convexity এত দামি, সেটি পর্ব ১২ (Convexity — কেন কিছু loss সহজ আর কিছু কঠিন)-এর বিষয়।

৬ · Regularization মানে MAP — Gaussian আর Laplace prior

AI-এর ভিত্তি কোর্সের পর্ব ৫০-এ Bayes-এর নিয়ম আর পর্ব ৬১-এর “আরও গভীরে” বাক্সে এই সংযোগের এক অনুচ্ছেদের সারাংশ পেয়েছেন। এবার পুরো হিসাব। মডেল: y = Xw + ε, ε-এর ঘরগুলো স্বাধীন N(০, σ²). তাহলে likelihood p(y | w) ∝ exp(−‖y − Xw‖²/(2σ²)). এখন w-এর প্রতিটি ঘরে একটি prior বসাই। Bayes-এর নিয়মে posterior ∝ likelihood × prior, আর posterior যেখানে সর্বোচ্চ, সেটিই MAP. লগ নিয়ে চিহ্ন উল্টালে, “সর্বোচ্চ posterior” মানে “সর্বনিম্ন −log posterior”:

−log p(w | y) = ‖y − Xw‖²/(2σ²) − log p(w) + ধ্রুবক

এই দৃষ্টিভঙ্গি কয়েকটি ব্যবহারিক নিয়মকে যুক্তি দেয়। Standardization: সব w_j-এ একই prior বসানো যুক্তিসঙ্গত কেবল যখন feature-গুলো একই মাপে; না হলে কিলোমিটার আর মিটারের coefficient-এর একই prior অর্থহীন। Intercept-এ শাস্তি নেই: intercept-এর prior সমতল (প্রায় অসীম τ), কারণ গড় স্তর নিয়ে আমাদের কোনো পূর্ব ধারণা নেই। λ-এর অর্থ: noise আর prior-এর অনুপাত; ডেটা যত বাড়ে, likelihood-এর পদ (n-টি বর্গের যোগফল) তত বড় হয়, আর একই λ-এর প্রভাব আপেক্ষিকভাবে কমে।

৭ · Early stopping — শাস্তি ছাড়াই regularization

AI-এর ভিত্তি কোর্সের পর্ব ৩৪-এর gradient descent চালান ½‖y − Xw‖²-এ, শুরু w₀ = ০ থেকে, ধাপের মাপ η: w_t+1 = w_t + ηXᵀ(y − Xw_t). কোনো শাস্তি নেই, তবু আগে থামালে w ছোট থাকে। কতটা ছোট? SVD-র ভাষায় প্রতিটি দিক আলাদা চলে। v_j দিকের অংশ θ_j(t) আর least squares-এর মান θ̂_j-এর পার্থক্য প্রতি ধাপে (1 − ησ_j²) গুণ হয়, আর শুরুতে পার্থক্য −θ̂_j. তাই:

θ_j(t) = [1 − (1 − ησ_j²)^t] · θ̂_j    বনাম ridge: [σ_j²/(σ_j² + λ)] · θ̂_j

দুটিই একটি filter: বড় σ-র দিক আগে শেখা হয়, ছোট σ-র দিক দেরিতে। মিলটি সংখ্যায় দেখতে λ = 1/(ηt) বসান। ησ_j²t ছোট হলে দুটি গুণকই প্রায় ησ_j²t; বড় হলে দুটিই প্রায় ১. মাঝখানে পার্থক্য আছে: x = ησ²t লিখলে GD-র গুণক প্রায় 1 − e^−x, ridge-এর x/(1 + x); x = ১-এ ০.৬৩ বনাম ০.৫. তাই early stopping মোটামুটি ridge, হুবহু নয়; ধাপের সংখ্যা t উল্টো λ-এর মতো কাজ করে। Early stopping এক ধরনের implicit regularization: loss-এ কিছু যোগ না করেও optimizer-এর আচরণ থেকে regularization আসে। শর্ত মনে রাখুন: η < 2/σ₁² না হলে |1 − ησ²| > ১ হয়ে GD ছিটকে যায়; ধাপের মাপের পুরো আলোচনা পর্ব ১৩ (Gradient descent-এর পরিবার)-এ।

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(3)            # কাল্পনিক ৪০টি সারি, ৬টি standardized feature
n, p = 40, 6
X = rng.normal(size=(n, p)) @ np.diag([3, 2, 1.5, 1, 0.5, 0.2])   # দিকভেদে আলাদা σ
X -= X.mean(0)
y = X @ np.array([1.0, -1.0, 0.8, 0.5, 0.0, 0.0]) + rng.normal(0, 1, n); y -= y.mean()
U, s, Vt = np.linalg.svd(X, full_matrices=False)
eta = 0.9 / s[0]**2                        # ধাপের মাপ: η < 2/σ_max² হলে স্থির
w_ols = np.linalg.lstsq(X, y, rcond=None)[0]

def ridge(lam):                            # ‖y − Xw‖² + λ‖w‖²
    return np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

w, gd = np.zeros(p), []
for t in range(1, 2001):                   # (1/2)‖y − Xw‖²-এ gradient descent, শুরু 0 থেকে
    w = w + eta * X.T @ (y - X @ w); gd.append(w.copy())
print(" t   λ=1/(ηt)   ‖w_t‖   ‖w_ridge‖   ‖w_t − w_ridge‖/‖w_ols‖")
for t in (5, 20, 100, 500, 2000):
    lam = 1 / (eta * t); wr = ridge(lam)
    print(f"{t:4d} {lam:9.2f} {np.linalg.norm(gd[t-1]):7.3f} {np.linalg.norm(wr):9.3f}"
          f" {np.linalg.norm(gd[t-1] - wr) / np.linalg.norm(w_ols):12.3f}")
print("σ:", np.round(s, 2)); t = 20; print("filter GD   :", np.round(1 - (1 - eta * s**2)**t, 3))
print("filter ridge:", np.round(s**2 / (s**2 + 1 / (eta * t)), 3))
gd = np.array(gd); rp = np.array([ridge(1 / (eta * t)) for t in range(1, 2001)])
fig, ax = plt.subplots(1, 2, figsize=(7, 3.2), sharey=True, layout="constrained")
ax[0].plot(np.linalg.norm(gd, axis=1), gd); ax[0].set_title("gradient descent, t = 1..2000")
ax[1].plot(np.linalg.norm(rp, axis=1), rp); ax[1].set_title("ridge, lambda = 1/(eta t)")
for a in ax: a.set_xlabel("norm of w")
ax[0].set_ylabel("coefficient"); plt.show()

 t   λ=1/(ηt)   ‖w_t‖   ‖w_ridge‖   ‖w_t − w_ridge‖/‖w_ols‖
   5     76.20   1.553     1.237        0.195
  20     19.05   1.648     1.510        0.098
 100      3.81   1.653     1.620        0.025
 500      0.76   1.655     1.647        0.012
2000      0.19   1.657     1.654        0.011
σ: [18.52 12.57  9.61  7.21  3.11  0.98]
filter GD   : [1.    1.    0.996 0.947 0.403 0.049]
filter ridge: [0.947 0.892 0.829 0.732 0.337 0.048]

বাস্তবে t বাছা হয় validation দিয়ে: প্রতি epoch-এর পর validation loss দেখুন, কয়েক epoch (patience) উন্নতি না হলে থামুন, আর সেরা অবস্থার ওজন ফেরত নিন। একটি বাড়তি সুবিধা: একটিমাত্র training-এ পুরো “path” পাওয়া যায়, প্রতিটি λ-এর জন্য আলাদা fit লাগে না। Neural network-এ এটি সবচেয়ে প্রচলিত regularization-গুলোর একটি; সেখানে loss non-convex, তাই উপরের হুবহু সূত্র খাটে না, শুধু ধারণাটি। এর বিস্তার পরে Deep Learning কোর্সে।

৮ · Elastic net, একসঙ্গে চলা feature আর cross-validation দিয়ে λ বাছাই

Lasso-র দুটি পরিচিত দুর্বলতা। প্রথম, p > n হলে সে বড়জোর n-টি feature বেছে নিতে পারে। দ্বিতীয়, প্রায় একই কয়েকটি feature থাকলে সে প্রায় খেয়ালখুশিমতো একটি রাখে, বাকিগুলো শূন্য করে, আর ডেটা একটু বদলালে অন্যটি রাখে। Hui Zou আর Trevor Hastie ২০০৫ সালে এর প্রতিকার হিসেবে elastic net প্রস্তাব করেন: দুটি শাস্তি একসঙ্গে। scikit-learn-এর ElasticNet (ডকুমেন্টেশন অনুযায়ী) কমায়:

(1/2n)‖y − Xw‖² + α·ρ·‖w‖₁ + ½·α·(1 − ρ)·‖w‖₂²

এখানে ρ হলো l1_ratio: ρ = ১ মানে lasso, ρ = ০ মানে ridge. L2 অংশটি loss-কে strictly convex করে, তাই সমাধান একটিই থাকে; আর Zou ও Hastie দেখান যে দুটি feature-এর সম্পর্ক ১-এর কাছে গেলে তাদের coefficient-ও কাছাকাছি হয়, যাকে তাঁরা grouping effect বলেছেন। L1 অংশ sparsity রাখে। নিচের ব্লকে তিনটি কলাম (app-এ ঢোকা, দেখা পাতা, কাটানো মিনিট) একটি অদৃশ্য “আগ্রহ”-এর তিনটি মাপ, পরস্পর সম্পর্ক প্রায় ০.৯৯; বাকি পাঁচটি অকাজের। λ (আর elastic net-এ ρ) বাছা হয় পাঁচ-fold cross-validation-এ, তারপর একই কাজ পাঁচটি bootstrap নমুনায়, কে টিকে থাকে দেখতে। max_iter=10000 দেওয়া হয়েছে যাতে solver পুরোপুরি মিলে যায় আর কোনো ConvergenceWarning না আসে।

import numpy as np
from sklearn.linear_model import LassoCV, ElasticNetCV
from sklearn.model_selection import KFold
np.set_printoptions(precision=2, suppress=True)

rng = np.random.default_rng(11)         # কাল্পনিক বাজারঘর: ৮০ জন গ্রাহকের মাসিক কেনাকাটা
n = 80
engage = rng.normal(size=n)                        # অদৃশ্য "আগ্রহ", তিনটি মাপে ধরা পড়ে
group = np.column_stack([engage + rng.normal(0, 0.1, n) for _ in range(3)])  # visits, pages, minutes
noise = rng.normal(size=(n, 5))                    # ৫টি অকাজের feature
X = np.column_stack([group, noise]); X = (X - X.mean(0)) / X.std(0)
y = 3 * engage + rng.normal(0, 1, n)
print("corr inside group:", np.round(np.corrcoef(group.T)[np.triu_indices(3, 1)], 2))

cv = KFold(5, shuffle=True, random_state=42)
las = LassoCV(cv=cv, max_iter=10000).fit(X, y)
enet = ElasticNetCV(l1_ratio=[0.1, 0.3, 0.5, 0.7, 0.9], cv=cv, max_iter=10000).fit(X, y)
print(f"LassoCV      alpha={las.alpha_:.3f}             coef:", las.coef_ + 0.0)
print(f"ElasticNetCV alpha={enet.alpha_:.3f} l1_ratio={enet.l1_ratio_}  coef:", enet.coef_ + 0.0)
print(f"CV MSE  lasso {las.mse_path_.mean(1).min():.3f}   enet {enet.mse_path_.mean(2).min():.3f}")

for seed in range(5):                   # একই সমস্যা, ৫টি bootstrap নমুনা: কে টিকে থাকে?
    i = np.random.default_rng(seed).integers(0, n, n)
    a = LassoCV(cv=cv, max_iter=10000).fit(X[i], y[i]).coef_[:3]
    b = ElasticNetCV(l1_ratio=0.5, cv=cv, max_iter=10000).fit(X[i], y[i]).coef_[:3]
    print(f"bootstrap {seed}: lasso group {a + 0.0}   enet(0.5) group {b + 0.0}")

corr inside group: [0.99 0.99 0.99]
LassoCV      alpha=0.159             coef: [ 1.23  0.19  0.85  0.   -0.08  0.    0.    0.  ]
ElasticNetCV alpha=0.190 l1_ratio=0.9  coef: [ 0.9   0.54  0.8   0.   -0.07  0.    0.    0.  ]
CV MSE  lasso 1.229   enet 1.208
bootstrap 0: lasso group [0.   0.   2.42]   enet(0.5) group [0.65 0.64 1.1 ]
bootstrap 1: lasso group [1.19 0.22 0.96]   enet(0.5) group [0.84 0.64 0.84]
bootstrap 2: lasso group [0.66 0.   1.6 ]   enet(0.5) group [0.78 0.46 0.98]
bootstrap 3: lasso group [0.39 0.   1.9 ]   enet(0.5) group [0.72 0.59 0.98]
bootstrap 4: lasso group [0.72 0.   1.62]   enet(0.5) group [0.77 0.56 1.02]

λ বাছাইয়ের নিয়ম

• পুরো path, warm start: LassoCV আর ElasticNetCV নিজেরাই α_max থেকে α_max/১০০০ পর্যন্ত (default eps=1e-3) ১০০টি α-র grid বানায় আর প্রতিটি fold-এ warm start দিয়ে path চালায়। ডকুমেন্টেশন l1_ratio-এর জন্য ১-এর কাছের মানগুলো (যেমন ০.১, ০.৫, ০.৭, ০.৯, ০.৯৫, ০.৯৯, ১) বেশি চেষ্টা করার পরামর্শ দেয়।
• One-standard-error rule: সর্বনিম্ন CV ভুলের α না নিয়ে, তার এক standard error-এর মধ্যে থাকা সবচেয়ে বড় α নেওয়া (সরলতম model). Breiman ও সহকর্মীদের CART বই (১৯৮৪) থেকে আসা এই নিয়ম প্রায় একই নির্ভুলতায় বেশি sparse model দেয়; তবে এটি একটি প্রথা, উপপাদ্য নয়।
• Leakage এড়ানো: standardization আর feature বাছাই প্রতিটি fold-এর ভেতরে, pipeline দিয়ে (AI-এর ভিত্তি কোর্সের পর্ব ৬৪). CV দিয়ে বাছা α-র score আশাবাদী; সৎ অনুমান চাইলে nested CV বা তালাবন্ধ test set (পর্ব ৬৬).
• Path পড়া: basis function আর interaction-সহ বড় linear model-এ regularization path কীভাবে পড়বেন, সেটি পর্ব ২৩ (Linear model গভীরে — basis function, interaction আর regularization path)-এর বিষয়।

৯ · AI-এর যুগে — regularization-এর কোড আর দাবি যাচাই

Regularization-এর কোড ছোট, কিন্তু ভুলগুলো নিঃশব্দ: ভুল objective-এর সংজ্ঞা, ভুল জায়গায় standardization, বা test set দিয়ে λ বাছাই। কোড চলে, সংখ্যাও আসে, শুধু সেই সংখ্যা যা বলছে তা সত্য নয়। AI assistant এই জায়গাগুলোতেই প্রায়ই আত্মবিশ্বাসী ভুল করে।

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

• Library-র objective (1/(2n) আছে কি না, l1_ratio-র সংজ্ঞা) নিজের version-এর docs থেকে পড়েছি
• ছোট একটি ডেটায় নিজের closed form বা coordinate descent-এর সঙ্গে library-র coefficient মিলিয়েছি
• Standardization আর feature বাছাই pipeline-এর ভেতরে, CV-র প্রতিটি fold-এ আলাদা করে হয়েছে
• λ বেছেছি CV-তে; test set তালাবন্ধ ছিল, শেষে একবার মেপেছি
• একটি সরল baseline (শুধু গড়, বা OLS) কয়েকটি seed বা fold-এ তুলনা করেছি
• Lasso-র বাছাই bootstrap-এ কতটা স্থির, দেখেছি; শূন্য coefficient-কে কারণগত প্রমাণ বলিনি
• কোনো warning চাপা দিইনি, আর কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি