Variational inference-এর ধারণা — ELBO আর approximation

AI-এর ভিত্তি কোর্সের পর্ব ৫০-এ prior, likelihood, evidence আর posterior এসেছে, আর পর্ব ৩৮-এ KL divergence. এই কোর্সের পর্ব ৩৯ EM-এর প্রমাণে ELBO ব্যবহার করেছে। পর্ব ৪৫ (Bayesian linear regression — posterior আর predictive অনিশ্চয়তা) আর পর্ব ৪৬ (Gaussian process — kernel দিয়ে অনিশ্চয়তাসহ অনুমান) এমন model দেখায় যেখানে posterior সূত্রে লেখা যায়। পর্ব ৪৭ (MCMC — Metropolis-Hastings আর Gibbs sampling) দেখায় সূত্র না থাকলে নমুনা দিয়ে কীভাবে চলে। এই পর্ব তৃতীয় পথ: posterior-কে একটি সহজ বণ্টন দিয়ে আসন্ন করা, আর সেই খোঁজকে optimization বানানো। চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python: NumPy, SciPy, matplotlib আর scikit-learn সহ)। প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় ১০ সেকেন্ড।

১ · Inference-এর সমস্যা — posterior কেন সরাসরি পাওয়া যায় না

Bayesian model-এ অজানা জিনিসগুলোকে (parameter বা লুকানো চলক) একসঙ্গে z বলি। Bayes-এর নিয়মে:

p(z | x) = p(x | z) p(z) ÷ p(x), p(x) = ∫ p(x | z) p(z) dz

লবটি সহজ: likelihood আর prior গুণ করলেই হয়। সমস্যা হরে। p(x)-কে বলে evidence বা marginal likelihood, আর এটি সব সম্ভাব্য z-এর ওপর একটি integral. Conjugate prior থাকলে integral-টি সূত্রে আসে, যেমন পর্ব ৪৫-এর Gaussian prior আর Gaussian noise. কিন্তু logistic regression-এর likelihood-এ sigmoid আছে, তাই Gaussian prior-এর সঙ্গে গুণ করলে চেনা কোনো বণ্টন আসে না। Mixture model-এ লুকানো দল-পরিচয়ের সম্ভাব্য বিন্যাস Kⁿ-টি; ৩টি দল আর ৩৬০ জন ক্রেতায় সংখ্যাটি প্রায় ১০^১৭২. দশ মাত্রার parameter-এ প্রতি অক্ষে ১০০টি grid বিন্দু নিলে integral-এ ১০^২০-টি বিন্দু লাগে।

Posterior হাতে না থাকলে তিনটি পথ খোলা থাকে:

VI-তে প্রশ্নটি দাঁড়ায় এভাবে: q* = arg min_{q ∈ Q} KL(q(z) ‖ p(z | x)). Q-কে বলে variational family, আর q-এর parameter-গুলো (যেমন একটি Gaussian-এর গড় আর sd) variational parameter. “Variational” নামটি এসেছে calculus of variations থেকে, যেখানে একটি সংখ্যা নয়, একটি পুরো function বদলে কিছু সর্বনিম্ন করা হয়। এখানে সেই function হলো বণ্টন q.

২ · KL-এর দুই দিক — mode-seeking বনাম mass-covering

AI-এর ভিত্তি কোর্সের পর্ব ৩৮ বলেছে: KL divergence কখনো ঋণাত্মক নয়, শূন্য শুধু দুটি বণ্টন এক হলে, আর এটি প্রতিসম নয়। প্রতিসম না হওয়াই এখানে মূল কথা:

KL(q ‖ p) = ∑ q(z) log [q(z) / p(z)], KL(p ‖ q) = ∑ p(z) log [p(z) / q(z)]

প্রথমটিতে গড় নেওয়া হয় q-এর ওপর। যেখানে q কিছু ভর রাখে কিন্তু p প্রায় শূন্য, সেখানে log(q/p) বিশাল। তাই KL(q ‖ p) কমাতে q এমন জায়গা এড়ায় যেখানে p নেই। একে বলে zero-forcing: p শূন্য হলে q-কেও শূন্য হতে হয়। দ্বিতীয়টিতে গড় নেওয়া হয় p-এর ওপর। যেখানে p ভর রাখে কিন্তু q প্রায় শূন্য, সেখানে log(p/q) বিশাল। তাই KL(p ‖ q) কমাতে q-কে p-এর সব ভর ঢাকতে হয়, প্রয়োজনে ফাঁকা জায়গাতেও ছড়িয়ে। নিচের stepper তিনটি অবস্থার একটি ছোট উদাহরণে দুটি দিক হাতে হিসাব করে।

এবার অবিরাম রূপ। বাজারঘরের অর্ডার-মূল্যের একটি কাল্পনিক দুই-চূড়ার বণ্টন p (শত টাকায়): ৬০% সাধারণ দিনের অর্ডার, গড় ৯, sd ২; ৪০% উৎসবের, গড় ২৪, sd ৩. Q হলো সব এক-চূড়ার Gaussian. নিচের ব্লক একটি grid-এ integral-কে যোগফল দিয়ে আসন্ন করে, আর scipy.optimize.minimize দিয়ে (m, log s) বদলে দুই দিকের KL আলাদা করে সর্বনিম্ন করে। KL(q ‖ p)-কে তিনটি ভিন্ন শুরু থেকে চালানো হয়েছে।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize
x = np.linspace(-10, 45, 4001); dx = x[1] - x[0]   # অর্ডারের মূল্য, শত টাকায়
p = 0.6 * norm.pdf(x, 9, 2) + 0.4 * norm.pdf(x, 24, 3)  # দুই চূড়ার লক্ষ্য p
def kl(a, b):                              # grid-এ ∫ a log(a/b) dx
    m = a > 1e-300
    return np.sum(a[m] * (np.log(a[m]) - np.log(b[m] + 1e-300))) * dx
def q_of(t):                               # t = (m, log s)
    return norm.pdf(x, t[0], np.exp(t[1]))
fits = {}
for start in (5.0, 15.0, 28.0):            # KL(q‖p): শুরু বদলালে ফল বদলায়
    r = minimize(lambda t: kl(q_of(t), p), [start, 1.0], method="Nelder-Mead")
    fits[f"KL(q||p), start {start:g}"] = r
r = minimize(lambda t: kl(p, q_of(t)), [15.0, 1.0], method="Nelder-Mead")
fits["KL(p||q)"] = r
for name, r in fits.items():
    print(f"{name:20s} m = {r.x[0]:6.2f}  s = {np.exp(r.x[1]):5.2f}  KL = {r.fun:.3f}")
mean = np.sum(x * p) * dx; sd = np.sqrt(np.sum((x - mean) ** 2 * p) * dx)
print(f"mean and sd of p:    m = {mean:6.2f}  s = {sd:5.2f}")
fig, ax = plt.subplots(figsize=(6, 3.6), layout="constrained")
ax.fill_between(x, p, alpha=0.3, color="gray", label="target p (bimodal)")
for name, st in (("KL(q||p), start 5", "-"), ("KL(q||p), start 28", "-"), ("KL(p||q)", "--")):
    ax.plot(x, q_of(fits[name].x), st, label="q from " + name)
ax.set_xlim(0, 40); ax.set_xlabel("order value (100 taka)"); ax.set_ylabel("density")
ax.set_title("Mode-seeking KL(q||p) vs mass-covering KL(p||q)"); ax.legend(fontsize=8)
plt.show()

KL(q||p), start 5    m =   9.02  s =  2.03  KL = 0.508
KL(q||p), start 15   m =   9.02  s =  2.03  KL = 0.508
KL(q||p), start 28   m =  23.92  s =  3.12  KL = 0.910
KL(p||q)             m =  15.00  s =  7.75  KL = 0.523
mean and sd of p:    m =  15.00  s =  7.75

VI কেন KL(q ‖ p) নেয়? কারণ ব্যবহারিক। KL(q ‖ p)-তে প্রত্যাশা q-এর ওপর, আর q আমাদের বানানো, তাই তার থেকে নমুনা তোলা আর গড় নেওয়া সহজ। অংশ ৩ দেখাবে যে এই দিকে অজানা p(x) একটি ধ্রুবক হয়ে আলাদা হয়ে যায়। KL(p ‖ q)-তে প্রত্যাশা posterior-এর ওপর, যা থেকে নমুনা তোলাই আসল সমস্যা। মূল্যটা হলো: VI-এর উত্তর একটি চূড়ায় আটকাতে পারে, আর অনিশ্চয়তা কম দেখায় (অংশ ৪ আর ৫)। Expectation propagation (Thomas Minka, ২০০১) নামের আরেকটি পদ্ধতি স্থানীয়ভাবে KL(p ‖ q)-ধরনের মাপ ব্যবহার করে, তাই তার উত্তর বেশি ছড়ানো হয়।

৩ · ELBO — log evidence = ELBO + KL

KL(q ‖ p(z | x)) সরাসরি হিসাব করা যায় না, কারণ তার ভেতরে অজানা posterior. সংজ্ঞা খুলে p(z | x) = p(x, z) / p(x) বসালে:

KL(q ‖ p(z | x)) = E_q[log q(z)] − E_q[log p(x, z)] + log p(x)

log p(x) z-এর ওপর নির্ভর করে না, তাই প্রত্যাশার বাইরে চলে আসে। সাজিয়ে লিখলে পর্বের কেন্দ্রীয় সমীকরণ:

log p(x) = ELBO(q) + KL(q ‖ p(z | x)), ELBO(q) = E_q[log p(x, z)] − E_q[log q(z)]

বাঁ দিক q-এর ওপর নির্ভর করে না, এটি ডেটা আর model-এর একটি স্থির সংখ্যা। তাই q বদলে ELBO বাড়ালে KL ঠিক ততটাই কমে। ELBO-তে আছে শুধু যৌথ p(x, z) = p(x | z)p(z), যা আমরা লিখতে পারি, আর q, যা আমাদের বানানো। KL ≥ ০ বলে ELBO ≤ log p(x), তাই নাম evidence lower bound. পর্ব ৩৯ একই নিম্নসীমা Jensen's inequality দিয়ে পেয়েছিল। এখানের পথে ফাঁকটির পরিচয়ও পাওয়া যায়: ফাঁক = KL.

ELBO-কে আরেকভাবে লেখা যায়: ELBO = E_q[log p(x, z)] + H(q), যেখানে H(q) হলো q-এর entropy. প্রথম পদ q-কে টানে সেই z-গুলোর দিকে যেখানে যৌথ সম্ভাবনা বেশি, অর্থাৎ চূড়ার দিকে। একা থাকলে q সংকুচিত হয়ে MAP বিন্দুতে বসত। দ্বিতীয় পদ q-কে ছড়াতে উৎসাহ দেয়। দুইয়ের টানাটানিতে q চূড়ার চারপাশে একটি যুক্তিসংগত চওড়ায় থামে।

এবার সংখ্যায়। বাজারঘরের একটি অর্ডার ৩ দিন দেরিতে পৌঁছেছে (x). কোন কুরিয়ার পাঠিয়েছিল, তা রেকর্ডে নেই (লুকানো z). দ্রুত কুরিয়ার অর্ডারের ৭০% নেয়, আর তার ১০% অর্ডার এমন দেরি করে; ধীর কুরিয়ার ৩০% নেয়, তার ৪০% দেরি করে। এখানে posterior হাতেই বের করা যায়, তাই সমীকরণটি হুবহু যাচাই করা যায়।

EM আর VI-এর সম্পর্ক। পর্ব ৩৯-এ θ ছিল একটি সংখ্যা (point estimate), আর লুকানো z-এর posterior হুবহু হিসাব করা যেত, তাই E-step-এ KL শূন্য হতো। VI-তে parameter-ও অনিশ্চিত, তাই সেগুলোও z-এর অংশ, আর posterior হুবহু পাওয়া যায় না। Q-এর ভেতরে KL সবচেয়ে কম করাই তখন সেরা চেষ্টা। দুটির মাঝামাঝি পদ্ধতির নাম variational EM: লুকানো চলকের জন্য VI, parameter-এর জন্য point estimate. Model বাছাইয়ে ELBO-কে প্রায়ই log p(x)-এর বদলে ব্যবহার করা হয়, তবে সতর্কতা দরকার: দুই model-এর ফাঁক (KL) আলাদা হতে পারে, তাই ELBO-র তুলনা আর evidence-এর তুলনা এক জিনিস নয়।

৪ · Mean-field আর CAVI — অজানা গড় আর precision-এর Normal model

Q কতটা বড় হবে, তা-ই VI-এর প্রধান নকশা-সিদ্ধান্ত। সবচেয়ে প্রচলিত পছন্দ mean-field approximation: অজানাগুলো q-এর অধীনে পরস্পর স্বাধীন, q(z) = q₁(z₁) q₂(z₂) ⋯ q_m(z_m). নামটি পদার্থবিদ্যা থেকে এসেছে, যেখানে প্রতিটি কণা বাকিদের প্রভাবকে একটি গড় “ক্ষেত্র” হিসেবে দেখে। এই গঠনের একটি বড় সুবিধা আছে। বাকি সব অংশ স্থির রাখলে একটি অংশের সেরা রূপ সরাসরি লেখা যায়:

log q_j*(z_j) = E_−j[log p(x, z)] + ধ্রুবক

E_−j মানে z_j ছাড়া বাকি সব অংশের বর্তমান q দিয়ে প্রত্যাশা। প্রতিটি অংশ পালাক্রমে এভাবে হালনাগাদ করাকে বলে CAVI (coordinate ascent variational inference). প্রতিটি ধাপ নিজের অংশে ELBO-কে হুবহু সর্বোচ্চ করে, তাই ELBO কখনো কমে না। এটি পর্ব ১৭-এর lasso-তে দেখা coordinate descent-এর ভাই, শুধু চলকগুলো এখানে পুরো বণ্টন। Conjugate model-এ প্রতিটি q_j* চেনা পরিবারের হয়, তাই হালনাগাদ মানে কয়েকটি সংখ্যা বদলানো।

উদাহরণ: অর্ডার-মূল্যের গড় আর ছড়ানো, দুটিই অজানা। বাজারঘরের ১০টি নতুন অর্ডার (হাজার টাকায়) x_i ~ N(μ, 1/τ), যেখানে τ = 1/σ² হলো precision. Prior: τ ~ Gamma(a₀, b₀) (Gamma distribution, b₀ হলো rate), আর μ | τ ~ N(μ₀, 1/(λ₀τ)). এই prior-এর আসল posterior জানা আছে (Normal-Gamma), তাই VI-এর ভুল সরাসরি মাপা যায়। এর ভেতরে μ আর τ নির্ভরশীল: τ বড় হলে (ছড়ানো কম) μ-এর অনিশ্চয়তাও কমে। Mean-field নেয় q(μ, τ) = q(μ)q(τ), আর উপরের সূত্র থেকে আসে (Bishop-এর ২০০৬ সালের বইয়ের অধ্যায় ১০-এর পরিচিত উদাহরণ):

দুটি একে অপরের ওপর নির্ভর করে, তাই পালাক্রমে হালনাগাদ। নিচের ব্লক q(τ)-কে ইচ্ছা করে খারাপ জায়গায় শুরু করায় (E[τ] = ০.১৫, অর্থাৎ sd প্রায় ২.৬ হাজার টাকা), প্রতি পাকে ELBO ছাপে, শেষে আসল log p(x)-এর সঙ্গে ফাঁক মাপে, আর দুই posterior-এর contour আঁকে। ELBO-র পদগুলো Gaussian আর Gamma-র পরিচিত প্রত্যাশা থেকে আসে; E[log τ] = ψ(a_N) − log b_N, যেখানে ψ হলো digamma.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammaln, digamma
from scipy.stats import norm, gamma
rng = np.random.default_rng(42)
x = rng.normal(1.2, 0.4, 10); N, xb = len(x), x.mean()   # ১০টি অর্ডার, হাজার টাকায়
m0, l0, a0, b0 = 1.0, 1.0, 2.0, 0.2       # prior: μ|τ ~ N(m0, 1/(l0·τ)), τ ~ Gamma(a0, b0)
mN = (l0 * m0 + N * xb) / (l0 + N)         # q(μ)-এর গড় প্রথম ধাপেই স্থির
aN, bN = a0 + (N + 1) / 2, 50.0            # q(τ)-এর শুরু: ইচ্ছা করে খারাপ, E[τ] = 0.15
def elbo(lN, aN, bN):
    Et, Elt = aN / bN, digamma(aN) - np.log(bN)
    S = np.sum((x - mN) ** 2) + N / lN; P = (mN - m0) ** 2 + 1 / lN
    return (0.5 * (N + 1) * (Elt - np.log(2 * np.pi)) + 0.5 * np.log(l0) - 0.5 * Et * (S + l0 * P)
            + a0 * np.log(b0) - gammaln(a0) + (a0 - 1) * Elt - b0 * Et
            + 0.5 * np.log(2 * np.pi * np.e / lN) + aN - np.log(bN) + gammaln(aN) + (1 - aN) * digamma(aN))
for it in range(1, 7):                     # CAVI: পালাক্রমে q(μ), তারপর q(τ)
    lN = (l0 + N) * aN / bN
    bN = b0 + 0.5 * (np.sum((x - mN) ** 2) + N / lN + l0 * ((mN - m0) ** 2 + 1 / lN))
    print(f"iter {it}: E[tau] = {aN / bN:6.3f}  ELBO = {elbo((l0 + N) * aN / bN, aN, bN):.4f}")
lN = (l0 + N) * aN / bN
A, B = a0 + N / 2, b0 + 0.5 * np.sum((x - xb) ** 2) + l0 * N * (xb - m0) ** 2 / (2 * (l0 + N))
logpx = (gammaln(A) - gammaln(a0) + a0 * np.log(b0) - A * np.log(B)
         + 0.5 * np.log(l0 / (l0 + N)) - N / 2 * np.log(2 * np.pi))       # আসল log evidence
print(f"exact log p(x) = {logpx:.4f}   gap (= KL) = {logpx - elbo(lN, aN, bN):.4f}")
print(f"sd of mu : exact {np.sqrt(B / ((l0 + N) * (A - 1))):.4f}  VI {1 / np.sqrt(lN):.4f}")
print(f"sd of tau: exact {np.sqrt(A) / B:.3f}   VI {np.sqrt(aN) / bN:.3f}")
mu, t = np.meshgrid(np.linspace(0.7, 1.45, 200), np.linspace(1, 19, 200))
exact = norm.pdf(mu, mN, 1 / np.sqrt((l0 + N) * t)) * gamma.pdf(t, A, scale=1 / B)
vi = norm.pdf(mu, mN, 1 / np.sqrt(lN)) * gamma.pdf(t, aN, scale=1 / bN)
fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
lv = exact.max() * np.array([0.1, 0.3, 0.6, 0.9])      # দুটিতে একই উচ্চতার রেখা
ax.contour(mu, t, exact, lv, colors="k"); ax.contour(mu, t, vi, lv, colors="r", linestyles="--")
ax.set_xlabel("mean mu (1000 taka)"); ax.set_ylabel("precision tau")
ax.set_title("Exact posterior (black) vs mean-field q (red dashed)"); plt.show()

iter 1: E[tau] =  1.802  ELBO = -11.0578
iter 2: E[tau] =  6.776  ELBO = -5.9109
iter 3: E[tau] =  8.304  ELBO = -5.7547
iter 4: E[tau] =  8.431  ELBO = -5.7537
iter 5: E[tau] =  8.439  ELBO = -5.7537
iter 6: E[tau] =  8.440  ELBO = -5.7537
exact log p(x) = -5.7184   gap (= KL) = 0.0353
sd of mu : exact 0.1121  VI 0.1038
sd of tau: exact 3.190   VI 3.082

কেন অনিশ্চয়তা কম দেখায়। একটি দুই-মাত্রার Gaussian লক্ষ্য নিন, যার দুই চলকের sd σ₁, σ₂ আর correlation ρ. KL(q ‖ p)-তে mean-field Gaussian-এর গড় ঠিক আসে, কিন্তু প্রতিটি অক্ষের sd হয় 1/√Λ_ii, যেখানে Λ = Σ⁻¹ হলো precision matrix. হিসাব করলে এটি σ_i√(1 − ρ²), অর্থাৎ শর্তসাপেক্ষ sd, প্রান্তিক sd নয় (পর্ব ১৯-এর শর্তসাপেক্ষ বণ্টন)। ρ = ০.৯ হলে √(১ − ০.৮১) = ০.৪৩৬: VI-এর sd আসলের অর্ধেকেরও কম। Zero-forcing-এর যুক্তিতে q হেলানো উপবৃত্তের বাইরে যেতে পারে না, তাই সে উপবৃত্তের ভেতরের সবচেয়ে বড় অক্ষ-বরাবর উপবৃত্তটি নেয়। পরের অংশের ব্লকে এই সূত্রটি সংখ্যায় মিলে যাবে।

৫ · Black-box VI আর reparameterisation trick

CAVI-এর জন্য প্রতিটি model-এ আলাদা derivation লাগে, আর conjugacy না থাকলে E_−j[…] সূত্রে আসে না। বাজারঘরের প্রশ্নটি নিন: ডেলিভারি যত দিন দেরি, ফেরতের সম্ভাবনা কি তত বাড়ে? Model: P(ফেরত) = σ(w₀ + w₁·দেরি), prior w ~ N(0, 3²I). এখানে ELBO-র প্রত্যাশা E_q[log p(y, w)]-এ sigmoid-এর log আছে, যার Gaussian গড় সূত্রে আসে না। সমাধান: প্রত্যাশাকে Monte Carlo নমুনার গড় দিয়ে আসন্ন করা, আর তার gradient দিয়ে পর্ব ১৩-এর মতো stochastic ধাপ নেওয়া। Model থেকে শুধু log p(x, z) আর তার gradient লাগে, তাই একে বলে black-box VI.

কঠিন অংশ হলো q-এর parameter-এর সাপেক্ষে gradient, কারণ নমুনাগুলো q থেকেই আসে। দুটি পরিচিত কৌশল আছে। Score-function estimator (REINFORCE নামেও পরিচিত) ∇ log q ব্যবহার করে। এটি যেকোনো q-এ চলে, কিন্তু gradient-এর variance খুব বেশি। Reparameterisation trick নমুনাকে একটি নির্দিষ্ট রূপান্তর হিসেবে লেখে: ε ~ N(0, 1) তুলে w = m + s·ε. এখন এলোমেলোতা ε-এ, যা m বা s-এর ওপর নির্ভর করে না। তাই chain rule সরাসরি চলে: ∂/∂m = ∇_wlog p, আর ∂/∂s = ∇_wlog p · ε. Entropy-র অংশ Gaussian-এ সূত্রে আসে, H = log s + ধ্রুবক, তাই ∂H/∂s = 1/s. নিচের stepper একটি নমুনায় পুরো হিসাব দেখায়।

এবার আসল কাজ। ৪০টি কাল্পনিক অর্ডারের দেরি (০ থেকে ৮ দিন) আর ফেরত (১৯টি ফেরত), q(w) = N(m, diag(s²)), অর্থাৎ mean-field Gaussian. প্রতি ধাপে ১০টি ε, ৮০০টি ধাপ, আর পর্ব ১৪-এর Adam (নিজে লেখা, ছোট হতে থাকা ধাপসহ)। তুলনার জন্য একই posterior-এ ২০,০০০ ধাপের random-walk Metropolis, যার বিস্তারিত পর্ব ৪৭-এর বিষয়; প্রথম ২,০০০টি নমুনা বাদ।

import numpy as np, time
rng = np.random.default_rng(42)            # ৪০টি অর্ডার: ডেলিভারি কত দিন দেরি, ফেরত এল কি না
d = rng.integers(0, 9, 40); y = (rng.random(40) < 1 / (1 + np.exp(-(-2.5 + 0.5 * d)))).astype(float)
X = np.column_stack([np.ones(40), d])
def logp(w):                               # log prior N(0, 3²I) + log likelihood
    z = X @ w; return y @ z - np.logaddexp(0, z).sum() - w @ w / 18
def grad(W):                               # প্রতিটি সারি একটি w-এর জন্য gradient
    return (y - 1 / (1 + np.exp(-W @ X.T))) @ X - W / 9
t0 = time.perf_counter()
m, rho = np.zeros(2), np.full(2, -1.0)     # q(w) = N(m, diag(s²)), s = exp(rho)
gm, gr, vm, vr = (np.zeros(2) for _ in range(4))
for step in range(1, 801):                 # Adam দিয়ে Monte Carlo ELBO বাড়ানো
    eps = rng.standard_normal((10, 2)); s = np.exp(rho)
    G = grad(m + s * eps)                  # reparameterisation: w = m + s·ε
    dm, dr = G.mean(0), (G * eps * s).mean(0) + 1     # +1 আসে entropy-র log s থেকে
    for g, a, v in ((dm, gm, vm), (dr, gr, vr)):
        a[:] = 0.9 * a + 0.1 * g; v[:] = 0.999 * v + 0.001 * g ** 2
    lr = 0.05 / np.sqrt(step / 100 + 1)
    m += lr * (gm / (1 - 0.9 ** step)) / (np.sqrt(vm / (1 - 0.999 ** step)) + 1e-8)
    rho += lr * (gr / (1 - 0.9 ** step)) / (np.sqrt(vr / (1 - 0.999 ** step)) + 1e-8)
t_vi = time.perf_counter() - t0
t0 = time.perf_counter()
w, lw, chain, acc = np.zeros(2), logp(np.zeros(2)), [], 0
for i in range(20000):                     # random-walk Metropolis, একই posterior
    prop = w + 0.25 * rng.standard_normal(2); lp = logp(prop)
    if np.log(rng.random()) < lp - lw: w, lw, acc = prop, lp, acc + 1
    chain.append(w)
C = np.array(chain[2000:]); t_mh = time.perf_counter() - t0
print(f"VI  : mean {m.round(3)} sd {np.exp(rho).round(3)} corr 0 (mean-field)  {t_vi:.2f} s")
print(f"MCMC: mean {C.mean(0).round(3)} sd {C.std(0).round(3)} corr {np.corrcoef(C.T)[0, 1]:.2f}"
      f"  {t_mh:.2f} s, accept {acc / 20000:.2f}")

VI  : mean [-2.541  0.581] sd [0.379 0.077] corr 0 (mean-field)  0.13 s
MCMC: mean [-2.698  0.605] sd [0.915 0.191] corr -0.90  0.67 s, accept 0.35

৬ · VI বনাম MCMC — কখন কোনটি

দুটি পদ্ধতির লক্ষ্য একই posterior, কিন্তু ভুলের ধরন আলাদা। MCMC-র ভুল আসে সীমিত নমুনা থেকে, আর সময় বাড়ালে কমে। VI-এর ভুল আসে Q-এর আকার থেকে, আর optimization যত ভালোই হোক, Q-এর বাইরের আকার ধরা যায় না। এই পার্থক্যকে বলে variance বনাম bias.

সম্প্রদায়ে এ নিয়ে মতভেদ আছে। অনেক পরিসংখ্যানবিদ চূড়ান্ত বিশ্লেষণে MCMC (বিশেষ করে HMC/NUTS) পছন্দ করেন আর VI-কে দ্রুত খসড়া মনে করেন। বড় মাপের ML-এ, যেমন লাখ লাখ নথির topic model বা neural network-এর ভেতরের লুকানো চলক, MCMC প্রায়ই ব্যবহারিক নয়, তাই VI-ই সাধারণ পছন্দ। মাঝামাঝি একটি অভ্যাস হলো ছোট নমুনায় দুটি চালিয়ে মিলিয়ে দেখা, ঠিক অংশ ৫-এর ব্লকের মতো। VI-এর ফল কতটা বিশ্বাসযোগ্য, তা মাপার কিছু diagnostic-ও প্রস্তাব করা হয়েছে, যেমন Yao ও সহকর্মীদের (২০১৮) Pareto-smoothed importance sampling-ভিত্তিক মাপ।

৭ · VI কোথায় চলে — BayesianGaussianMixture, LDA আর VAE

আপনি হয়তো আগেই VI ব্যবহার করেছেন, না জেনে। scikit-learn-এর BayesianGaussianMixture ওজন, গড় আর covariance-এর ওপর prior বসায় আর mean-field VI দিয়ে fit হয়। গঠনটি পর্ব ৩৯-এর EM-এর মতো পালাক্রমে চলে, শুধু point estimate-এর বদলে প্রতিটি parameter-এর একটি q. গুরুত্বপূর্ণ সেটিং weight_concentration_prior: ছোট মান prior-কে বলে “অল্প কয়েকটি component-এ ভার জমুক”, বড় মান বলে “সবাইকে ভাগ দাও”। weight_concentration_prior_type ডিফল্টে "dirichlet_process" (Dirichlet process-এর stick-breaking রূপ, n_components পর্যন্ত কাটা); বিকল্প "dirichlet_distribution", যা সাধারণ Dirichlet distribution. নিচের ব্লক পর্ব ৩৯-এর তিন দলের কাল্পনিক ক্রেতা-ডেটায় ৮টি component দিয়ে শুরু করে তিনটি সেটিং তুলনা করে। আগের পর্বের মতোই init_params="random_from_data", কারণ ডিফল্ট "kmeans" ব্রাউজারের Python-এ একটি অপ্রাসঙ্গিক warning ছাপে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import BayesianGaussianMixture
rng = np.random.default_rng(42)            # পর্ব ৩৯-এর কাল্পনিক ক্রেতা-ডেটা, ৩টি আসল দল
X = np.vstack([rng.multivariate_normal([6, 3], [[4, 1], [1, 1]], 150),
               rng.multivariate_normal([15, 7], [[9, 4], [4, 4]], 150),
               rng.multivariate_normal([28, 4.5], [[30, -6], [-6, 3]], 60)])
fig, ax = plt.subplots(figsize=(6, 3.4), layout="constrained")
for j, (kind, c) in enumerate([("dirichlet_process", 0.001), ("dirichlet_distribution", 0.001),
                               ("dirichlet_distribution", 100.0)]):
    bgm = BayesianGaussianMixture(n_components=8, weight_concentration_prior_type=kind,
                                  weight_concentration_prior=c, init_params="random_from_data",
                                  max_iter=1000, random_state=0).fit(X)
    w = np.sort(bgm.weights_)[::-1]
    print(f"{kind[10:]:12s} prior {c:g}: weights {w.round(3)}")
    print(f"   components above 1%: {(w > 0.01).sum()}, iterations {bgm.n_iter_}")
    ax.bar(np.arange(8) + 0.27 * (j - 1), w, width=0.27, label=f"{kind[10:]}, {c:g}")
ax.set_xlabel("component (sorted by weight)"); ax.set_ylabel("weight")
ax.set_title("BayesianGaussianMixture with 8 components"); ax.legend(fontsize=8)
plt.show()

process      prior 0.001: weights [0.428 0.394 0.166 0.003 0.003 0.003 0.003 0.   ]
   components above 1%: 3, iterations 46
distribution prior 0.001: weights [0.426 0.409 0.165 0.    0.    0.    0.    0.   ]
   components above 1%: 3, iterations 37
distribution prior 100: weights [0.135 0.128 0.128 0.125 0.122 0.122 0.122 0.118]
   components above 1%: 8, iterations 116

LDA. পর্ব ৪৪-এর scikit-learn LatentDirichletAllocation variational Bayes দিয়ে শেখে (Blei, Ng আর Jordan, ২০০৩-এর মূল প্রবন্ধের পদ্ধতি, আর learning_method="online"-এ Hoffman, Blei আর Bach-এর ২০১০ সালের online রূপ)। প্রতিটি রিভিউয়ের topic-অনুপাত θ-এর জন্য একটি Dirichlet q, আর প্রতিটি শব্দের topic-এর জন্য একটি categorical q, mean-field-এ আলাদা। সেখানে দেখা perplexity এই ELBO থেকেই আসে, তাই সেটি আসন্ন মান। Stochastic VI: Hoffman ও সহকর্মীরা (২০১৩) দেখান যে ELBO-র gradient mini-batch থেকে অনুমান করা যায়, ফলে কোটি নথিতেও VI চলে।

Variational autoencoder (VAE). Kingma আর Welling (২০১৩ সালে প্রকাশিত, ২০১৪-এর ICLR) আর প্রায় একই সময়ে Rezende, Mohamed আর Wierstra (২০১৪) অংশ ৫-এর reparameterisation trick-কে neural network-এর সঙ্গে জোড়েন। প্রতিটি ডেটা-বিন্দুর জন্য আলাদা q-এর parameter optimize না করে একটি network (encoder) x থেকে সরাসরি q(z | x)-এর গড় আর sd দেয়। একে বলে amortized inference. আরেকটি network (decoder) z থেকে x ফিরিয়ে আনে, আর দুটিকে একসঙ্গে ELBO বাড়িয়ে শেখানো হয়। পরের deep learning কোর্স এটি বিস্তারিত দেখাবে। এখানে মনে রাখার কথা: VAE-র loss-এর “reconstruction + KL” দুই অংশ ঠিক অংশ ৩-এর ELBO.

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

VI-এর কোডে ভুল সহজে লুকায়, কারণ একটি ভুল ELBO-ও মসৃণভাবে বাড়তে পারে। AI assistant-এর উত্তরে সাধারণ ভুলগুলো:

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

• KL-এর দিক আর ELBO-র চিহ্ন একটি ছোট discrete উদাহরণে হাতে মিলিয়েছি (ELBO + KL = log p(x))
• Gradient-কে finite difference দিয়ে যাচাই করেছি, আর ELBO-র বক্ররেখা দেখেছি (CAVI-তে কখনো কমবে না; stochastic VI-তে গড়ে বাড়বে)
• কয়েকটি seed আর শুরুতে চালিয়ে দেখেছি ফল একটি চূড়ায় আটকাচ্ছে কি না
• ছোট নমুনায় MCMC-র সঙ্গে গড়, sd আর correlation তুলনা করেছি; VI-এর sd সরাসরি রিপোর্ট করার আগে সতর্ক থেকেছি
• প্রতিটি scikit-learn parameter-এর নাম আর ডিফল্ট নিজের version-এর docs-এ মিলিয়েছি, আর prior-এর কয়েকটি মান চালিয়েছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি