Matrix calculus — gradient, Jacobian আর Hessian

AI-এর ভিত্তি কোর্সের পর্ব ৩৪-এ derivative, partial derivative, gradient, chain rule আর gradient descent এক বা দুই চলকে এসেছিল; পর্ব ৫৮-এ normal equations আর পর্ব ৫৯-এ logistic regression-এর gradient সূত্র হিসেবে দেখেছেন। এই পর্বের প্রশ্ন: হাজার হাজার parameter থাকলে সেই সূত্রগুলো এক লাইনে কীভাবে লেখা আর derive করা যায়, আর লেখা সূত্রটি ঠিক কি না, তা কীভাবে নিশ্চিত হবেন? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় দশ সেকেন্ড।

১ · কেন matrix calculus — একটি একটি করে নয়, পুরো vector একসঙ্গে

ML-এ প্রায় সব শেখা একটি সংখ্যা কমানো: loss L(w), যেখানে w হলো সব parameter মিলে একটি লম্বা vector. Linear regression-এ w-এর দৈর্ঘ্য feature-এর সংখ্যার সমান, ছোট neural network-এ কয়েক হাজার, বড় model-এ কোটি কোটি। প্রতিটি ঘরের partial derivative আলাদা করে লেখা অসম্ভব, আর কোডেও আমরা loop নয়, matrix গুণ লিখি। তাই derivative-ও লিখতে হয় matrix ভাষায়। এই পর্বে তিনটি বস্তু:

এই নিয়মের সবচেয়ে বড় উপকার: আকার মিলিয়ে ভুল ধরা যায়। Gradient descent-এর w ← w − η∇L লিখতে হলে ∇L-এর আকার w-এর সমান হতেই হবে; সূত্রের ডান দিকে যদি 1 × n বা n × n আসে, derivation-এ কোথাও ভুল আছে। এই পাতায় প্রতিটি সূত্রের পাশে আকারও লেখা থাকবে। চিহ্নের মানে: aᵀ মানে a-এর transpose, aᵀx মানে dot product, ‖v‖² = vᵀv মানে ঘরগুলোর বর্গের যোগফল, আর I হলো identity matrix.

২ · প্রথম তিনটি সূত্র — aᵀx, xᵀAx আর একটি নিরাপদ কৌশল

সূত্র ১: f(x) = aᵀx = a₁x₁ + … + a_nx_n. ∂f/∂x_i-তে শুধু a_ix_i পদটি x_i-এর ওপর নির্ভর করে, তাই উত্তর a_i. সব ঘর মিলে ∇(aᵀx) = a, আকার n × 1. এক চলকের “ax-এর derivative a”-র হুবহু রূপ।

সূত্র ২: f(x) = xᵀAx, যেখানে A একটি n × n matrix. একে বলে quadratic form: সব পদ দ্বিঘাত, f = ∑_i∑_j a_ijx_ix_j. এর আকার (বাটি, জিন, উপত্যকা) পর্ব ১০ (Positive definite matrix — quadratic form আর covariance)-এর বিষয়; এখানে শুধু তার ঢাল। নিচের stepper একটি 2 × 2 উদাহরণে প্রতিটি ঘর আলাদা করে derivative নেয়, তারপর ধাঁচটি matrix আকারে চিনে নেয়।

এক চলকের ax²-এর derivative 2ax; xᵀAx-এর ক্ষেত্রে তার সঠিক রূপ (A + Aᵀ)x. ঘরে ঘরে হিসাব সবসময় কাজ করে, কিন্তু বড় সূত্রে ক্লান্তিকর। এর চেয়ে নিরাপদ আর দ্রুত একটি কৌশল আছে, যাকে বলে differential পদ্ধতি।

এই কৌশলের দুটি সুবিধা: যোগফলের সূচক (i, j) লিখতে হয় না, আর শেষে gᵀdx আকার পেলে layout নিয়ে সন্দেহ থাকে না। আসলে এটি সেই পুরোনো সংজ্ঞাই: f(x + dx) ≈ f(x) + ∇f(x)ᵀdx. অর্থাৎ gradient হলো সেই vector, যার সঙ্গে dot product নিলে যেকোনো ছোট সরণের জন্য f-এর বদল পাওয়া যায়। একক দিক u-তে এই বদলের হার ∇fᵀu-কে বলে directional derivative; u যখন ∇f-এর দিকে, তখন সবচেয়ে বেশি। এজন্যই gradient সবচেয়ে খাড়া চড়াইয়ের দিক দেখায়, আর অংশ ৬-এর ছবিতে দেখবেন সেটি contour-এর সঙ্গে লম্ব।

চলক যখন একটি matrix, আর log det

Matrix চলক W (আকার p × q)-এর scalar function-এও একই কৌশল খাটে, শুধু dot product-এর জায়গায় বসে tr(Gᵀ dW) = ∑_ij G_ij dW_ij (tr মানে trace, অর্থাৎ কর্ণের যোগফল)। উদাহরণ: অনেক output-এর linear regression-এ loss f(W) = ‖XW − Y‖_F², যেখানে Frobenius norm ‖M‖_F² মানে matrix-এর সব ঘরের বর্গের যোগফল। R = XW − Y লিখলে df = 2 tr(Rᵀ dR) = 2 tr(Rᵀ X dW) = tr((2XᵀR)ᵀ dW), তাই ∇_W‖XW − Y‖_F² = 2Xᵀ(XW − Y), আকার p × q, W-এর মতোই।

আরেকটি সূত্র এই কোর্সে বারবার আসবে, বিশেষ করে multivariate Gaussian-এর likelihood-এ (পর্ব ১৯): বর্গাকার, inverse আছে এমন S-এর জন্য ∇_S log det S = S^−ᵀ (S প্রতিসম হলে সেটি S⁻¹-এর সমান)। এখানে det মানে determinant, আর log det বড় matrix-এ সরাসরি det না নিয়ে np.linalg.slogdet দিয়ে হিসাব করা হয়। এর প্রমাণ Jacobi-র একটি সূত্র থেকে আসে, d(det S) = det S · tr(S⁻¹dS); এখানে আমরা শুধু সংখ্যা দিয়ে যাচাই করব।

৩ · ‖Ax − b‖² আর chain rule — ভেতর থেকে বাইরে, বাইরে থেকে ভেতরে

Least squares-এর loss f(x) = ‖Ax − b‖², যেখানে A একটি m × n matrix (m-টি সমীকরণ, n-টি অজানা), b দৈর্ঘ্য m-এর vector. এটিকে দুটি ধাপে ভাঙা যায়: প্রথমে r = Ax − b (residual, আকার m × 1), তারপর f = rᵀr. প্রতিটি ধাপের derivative সহজ: r-এর x-এর সাপেক্ষে Jacobian হুবহু A (m × n), আর f-এর r-এর সাপেক্ষে gradient 2r. Matrix আকারের chain rule এ দুটিকে জোড়ে:

∇_xf = J_r(x)ᵀ · ∇_rf = Aᵀ · 2(Ax − b)

আকার মিলিয়ে দেখুন: Aᵀ হলো n × m, 2r হলো m × 1, গুণফল n × 1, ঠিক x-এর আকার। Transpose কেন? Scalar f-এর Jacobian হলো সারি: J_f∘r = J_f(r) · J_r(x) = (2r)ᵀA, আকার 1 × n; gradient তার transpose, Aᵀ(2r). নিচের চিত্র আর stepper এই দুই দিকের হিসাব দেখায়।

শেষ ধাপের “residual কলামগুলোর সঙ্গে লম্ব” কথাটির জ্যামিতিক ছবি, অর্থাৎ least squares মানে b-কে A-এর column space-এ projection, পর্ব ০৭ (Matrix-এর জ্যামিতি — projection, rotation আর least squares)-এর বিষয়। এখানে একই উত্তর এসেছে calculus দিয়ে: ঢাল শূন্য যেখানে, সেখানেই সর্বনিম্ন। সর্বনিম্নই কেন, সর্বোচ্চ বা জিন নয়, তার উত্তর দেবে Hessian (অংশ ৬)।

৪ · সূত্রকে বিশ্বাস নয়, যাচাই — central difference দিয়ে gradient checking

হাতে derive করা gradient-এ চিহ্ন, transpose বা ২-এর গুণক ভুল হওয়া খুব স্বাভাবিক, আর ভুল gradient দিয়েও training চলে, শুধু খারাপ ফল দেয়। তাই নিয়ম: প্রতিটি নতুন gradient সূত্র সংখ্যা দিয়ে মিলিয়ে নিন। একে বলে gradient checking. AI-এর ভিত্তি কোর্সের পর্ব ৩৪-এ এক চলকে দেখেছেন, h খুব বড় হলে সূত্রের ভুল আর খুব ছোট হলে float-এর rounding ভুল। বহু চলকে প্রতিটি ঘর i-এর জন্য আলাদা করে:

∂f/∂x_i ≈ [ f(x + h·e_i) − f(x − h·e_i) ] / (2h)

এখানে e_i হলো i-তম ঘরে 1 আর বাকি সব 0 এমন vector, আর h একটি ছোট সংখ্যা (এখানে 10⁻⁵ নেওয়া হয়েছে)। দুই দিকে সমান সরে ভাগ করাকে বলে central difference; এর সূত্রগত ভুল h²-এর অনুপাতে কমে, এক দিকের (forward) finite difference-এর ভুল কমে কেবল h-এর অনুপাতে। নিচের ব্লক একটি ছোট সাহায্যকারী function লিখে এই পর্বের ছয়টি সূত্র seed-সহ এলোমেলো ডেটায় মেলায়। ইচ্ছা করে একটি ভুল সূত্রও (অপ্রতিসম A-তে 2Ax) রাখা হয়েছে, যাতে দেখা যায় ভুল কেমন দেখায়।

import numpy as np

rng = np.random.default_rng(42)

def num_grad(f, x, h=1e-5):         # central difference, প্রতিটি ঘরে আলাদা করে
    g = np.zeros_like(x)
    for i in range(x.size):
        e = np.zeros_like(x); e.flat[i] = h
        g.flat[i] = (f(x + e) - f(x - e)) / (2 * h)
    return g

n = 4
x, a = rng.normal(size=n), rng.normal(size=n)
A = rng.normal(size=(n, n))                     # ইচ্ছা করে অপ্রতিসম
B, b = rng.normal(size=(6, n)), rng.normal(size=6)
X, Y = rng.normal(size=(5, 3)), rng.normal(size=(5, 2))
W = rng.normal(size=(3, 2))                     # এবার চলকটি নিজেই একটি matrix
S = B.T @ B + np.eye(n)                         # প্রতিসম, positive definite

checks = [
    ("a.T x",         lambda v: a @ v,               x, a),
    ("x.T A x",       lambda v: v @ A @ v,           x, (A + A.T) @ x),
    ("x.T A x vs 2Ax", lambda v: v @ A @ v,          x, 2 * A @ x),
    ("|Bx - b|^2",    lambda v: np.sum((B @ v - b)**2), x, 2 * B.T @ (B @ x - b)),
    ("|XW - Y|_F^2",  lambda V: np.sum((X @ V - Y)**2), W, 2 * X.T @ (X @ W - Y)),
    ("log det S",     lambda M: np.linalg.slogdet(M)[1], S, np.linalg.inv(S).T),
]
for name, f, at, formula in checks:
    num = num_grad(f, at)
    err = np.max(np.abs(num - formula))
    rel = err / np.max(np.abs(num))
    print(f"{name:<15} shape {str(at.shape):<7} max err {err:.1e}  rel {rel:.1e}")

a.T x           shape (4,)    max err 5.6e-12  rel 2.9e-12
x.T A x         shape (4,)    max err 4.5e-11  rel 1.0e-11
x.T A x vs 2Ax  shape (4,)    max err 2.2e+00  rel 5.1e-01
|Bx - b|^2      shape (4,)    max err 7.9e-11  rel 1.2e-11
|XW - Y|_F^2    shape (3, 2)  max err 6.7e-11  rel 1.3e-11
log det S       shape (4, 4)  max err 4.0e-11  rel 6.5e-11

৫ · Gradient শূন্য — least squares আর ridge-এর closed form

Linear regression-এ design matrix X (n × d, প্রতিটি সারি একজন গ্রাহক), label y (n × 1), weight w (d × 1). অংশ ৩-এর সূত্রে A-এর জায়গায় X আর b-এর জায়গায় y বসালেই:

L(w) = ‖Xw − y‖² ⇒ ∇L = 2Xᵀ(Xw − y) = 0 ⇒ XᵀXw = Xᵀy

এগুলোই normal equations; X-এর কলামগুলো স্বাধীন হলে XᵀX-এর inverse আছে, আর সমাধান একটিই। Ridge-এ loss-এ ওজনের আকারের শাস্তি যোগ হয়: L(w) = ‖Xw − y‖² + α‖w‖². α‖w‖² = α wᵀIw, সূত্র ২-এ A = αI বসালে gradient 2αw. তাই:

∇L = 2Xᵀ(Xw − y) + 2αw = 0 ⇒ (XᵀX + αI) w = Xᵀy

এক লাইনের এই হিসাবে ridge-এর দুটি গুণ দেখা যায়। প্রথমত, XᵀX + αI-এর প্রতিটি eigenvalue XᵀX-এর চেয়ে α বেশি, তাই প্রায় একই দুটি কলাম থাকলেও matrix-টির inverse থাকে আর condition number ছোট হয়। দ্বিতীয়ত, এটি একটি closed form: কোনো পুনরাবৃত্তি লাগে না। কিন্তু scikit-learn-এর সঙ্গে মেলাতে গেলে দুটি খুঁটিনাটি জানা দরকার। scikit-learn-এর ডকুমেন্টেশন অনুযায়ী Ridge কমায় ‖y − Xw‖² + alpha·‖w‖², অর্থাৎ n দিয়ে ভাগ নেই আর ½ নেই; তাই alpha হুবহু আমাদের α. আর fit_intercept=True (default) হলে intercept-এ শাস্তি দেওয়া হয় না: ডেটা গড় থেকে সরিয়ে (centering) ridge সমাধান করা হয়, তারপর b = ȳ − x̄ᵀw.

import numpy as np
from sklearn.linear_model import Ridge

rng = np.random.default_rng(42)        # কাল্পনিক বাজারঘর: ৪০ জন গ্রাহক
n = 40
orders = rng.poisson(6, n)                         # গত ৬ মাসে অর্ডার
basket = rng.normal(900, 250, n)                   # গড় ঝুড়ি, টাকা
days = rng.integers(1, 60, n)                      # শেষ অর্ডারের পর দিন
X = np.column_stack([orders, basket / 100, days]).astype(float)
y = 300 + 120 * orders + 45 * basket / 100 - 8 * days + rng.normal(0, 150, n)

# ১) least squares: ∇ = 2Xᵀ(Xw − y) = 0  →  XᵀX w = Xᵀy
X1 = np.column_stack([np.ones(n), X])              # intercept-এর জন্য সব-1 কলাম
w_ne = np.linalg.solve(X1.T @ X1, X1.T @ y)
w_ls = np.linalg.lstsq(X1, y, rcond=None)[0]
print("OLS w (b, orders, basket, days):", np.round(w_ls, 3))
print("solve vs lstsq       :", f"{np.max(np.abs(w_ne - w_ls)):.1e}")
print("|gradient| at w_ls   :", f"{np.max(np.abs(2 * X1.T @ (X1 @ w_ls - y))):.1e}")

# ২) ridge: ‖y − Xw‖² + α‖w‖²,  ∇ = 0  →  (XᵀX + αI) w = Xᵀy
alpha = 50.0
w_r = np.linalg.solve(X.T @ X + alpha * np.eye(3), X.T @ y)
sk = Ridge(alpha=alpha, fit_intercept=False).fit(X, y)
print("ridge, no intercept  :", f"{np.max(np.abs(w_r - sk.coef_)):.1e}")

# ৩) intercept থাকলে: কেন্দ্রে সরিয়ে solve, intercept-এ শাস্তি নেই
mx, my = X.mean(axis=0), y.mean()
Xc, yc = X - mx, y - my
w_c = np.linalg.solve(Xc.T @ Xc + alpha * np.eye(3), Xc.T @ yc)
b_c = my - mx @ w_c
sk = Ridge(alpha=alpha).fit(X, y)
print("ridge + intercept    :", f"{max(np.max(np.abs(w_c - sk.coef_)), abs(b_c - sk.intercept_)):.1e}")
print("ridge w:", np.round(w_c, 3), " b:", round(b_c, 2))
wrong = np.linalg.solve(X1.T @ X1 + alpha * np.eye(4), X1.T @ y)   # intercept-কেও শাস্তি
print("penalised intercept  :", np.round(wrong, 2))

OLS w (b, orders, basket, days): [283.455 120.362  58.024 -10.898]
solve vs lstsq       : 2.5e-12
|gradient| at w_ls   : 7.4e-09
ridge, no intercept  : 7.1e-14
ridge + intercept    : 1.4e-14
ridge w: [ 97.764  50.328 -11.384]  b: 501.29
penalised intercept  : [ 18.51 119.37  75.6   -8.26]

৬ · Hessian আর দ্বিতীয় ক্রমের Taylor — ঢালের ঢাল, আর বাঁকের মাপ

Gradient বলে কোন দিকে ঢাল, কিন্তু সেই ঢাল সামনে কতটা বদলাবে তা বলে না। বাটির তলা, জিনের মাঝখান আর সমতল মাঠ, তিন জায়গাতেই gradient শূন্য। পার্থক্য ধরে দ্বিতীয় derivative-এর matrix, Hessian:

[∇²f(x)]_ij = ∂²f / ∂x_i∂x_j, আকার n × n

Hessian আসলে gradient-এর Jacobian: ∇f: ℝⁿ → ℝⁿ-এর Jacobian নিলেই পাওয়া যায়। দ্বিতীয় derivative যথেষ্ট মসৃণ হলে ক্রম বদলালে ফল বদলায় না (Schwarz-এর উপপাদ্য), তাই Hessian প্রতিসম। কয়েকটি তৈরি উদাহরণ: aᵀx-এর Hessian 0; xᵀAx-এর A + Aᵀ (সব x-এ একই); ‖Ax − b‖²-এর 2AᵀA; ridge-এর 2(XᵀX + αI). এগুলোর সঙ্গে আসে Taylor-এর দ্বিতীয় ক্রমের approximation:

f(x + d) ≈ f(x) + ∇f(x)ᵀd + ½ dᵀ∇²f(x) d

প্রথম পদ এখনকার মান, দ্বিতীয়টি ঢাল ধরে সরলরেখায় অনুমান, তৃতীয়টি বাঁকের সংশোধন। Taylor approximation-এর ভুল প্রথম ক্রমে সাধারণত d-এর দৈর্ঘ্যের বর্গের অনুপাতে, দ্বিতীয় ক্রমে ঘনের অনুপাতে কমে। Quadratic function-এ (least squares, ridge) দ্বিতীয় ক্রমের approximation হুবহু সত্য, কোনো ভুল নেই। আর gradient শূন্য এমন বিন্দুতে তৃতীয় পদই সব বলে: Hessian positive definite হলে (সব eigenvalue ধনাত্মক) প্রতিটি দিকে f বাড়ে, তাই সেটি local minimum; কিছু eigenvalue ধনাত্মক আর কিছু ঋণাত্মক হলে saddle point. 2AᵀA সবসময় positive semidefinite, কারণ dᵀ(2AᵀA)d = 2‖Ad‖² ≥ 0; এজন্যই least squares-এর loss বাটির মতো, অর্থাৎ convex. Positive definite হওয়ার পরীক্ষাগুলো পর্ব ১০-এ, আর convexity-র পুরো তত্ত্ব পর্ব ১২ (Convexity — কেন কিছু loss সহজ আর কিছু কঠিন)-এ।

Jacobian আর Taylor, সংখ্যায়

এবার একটি অ-quadratic উদাহরণ, বাজারঘরের ডেলিভারি থেকে। শহরের তিনটি এলাকা d₁, d₂, d₃ (km-এ স্থানাঙ্ক), দিনে সেখানে ৩০, ৫০ আর ২০টি অর্ডার। নতুন হাব x-এ বসালে তিনটি দূরত্বের vector f(x) = (‖x − d₁‖, ‖x − d₂‖, ‖x − d₃‖), অর্থাৎ f: ℝ² → ℝ³, তাই Jacobian 3 × 2. ‖x − d_k‖-এর gradient হলো d_k থেকে x-এর দিকের একক vector u_k = (x − d_k)/‖x − d_k‖, তাই Jacobian-এর সারি k হলো u_kᵀ. মোট রাইড-km L(x) = ∑ w_k‖x − d_k‖-এর gradient chain rule-এ Jᵀw, আর Hessian ∑ w_k(I − u_ku_kᵀ)/‖x − d_k‖. নিচের ব্লক Jacobian সংখ্যা দিয়ে মেলায়, তারপর একটি দিকে ধাপ t ছোট করতে করতে দুই ক্রমের Taylor-এর ভুল দেখে।

import numpy as np
import matplotlib.pyplot as plt

D = np.array([[2.0, 1.0], [7.0, 2.0], [5.0, 6.0]])   # ৩টি ডেলিভারি এলাকা (km)
wt = np.array([30.0, 50.0, 20.0])                    # দিনে অর্ডার
dist = lambda x: np.linalg.norm(x - D, axis=1)        # f: R² → R³
def jac(x):                                           # সারি k = (x − d_k)ᵀ / ‖x − d_k‖
    return (x - D) / dist(x)[:, None]
L = lambda x: wt @ dist(x)                            # মোট রাইড-km
grad = lambda x: jac(x).T @ wt                        # chain rule: Jᵀ · ∂L/∂f
def hess(x):
    U = jac(x)
    return sum(w / r * (np.eye(2) - np.outer(u, u)) for w, r, u in zip(wt, dist(x), U))
x0, h = np.array([4.0, 3.0]), 1e-6                   # প্রস্তাবিত হাবের জায়গা
J_num = np.column_stack([(dist(x0 + e) - dist(x0 - e)) / (2 * h) for e in np.eye(2) * h])
print("J analytic:\n", np.round(jac(x0), 4))
print("max |J_num - J|:", f"{np.max(np.abs(J_num - jac(x0))):.1e}")

g, H = grad(x0), hess(x0)
print("L(x0) =", round(L(x0), 3), " grad =", np.round(g, 3))
u = np.array([0.6, 0.8])                              # একটি একক দিক
print("   t    true L  1st-order err  2nd-order err")
for t in (2.0, 1.0, 0.5, 0.25):
    d = t * u
    lin = L(x0) + g @ d
    quad = lin + 0.5 * d @ H @ d
    print(f"{t:>4}  {L(x0 + d):8.3f}   {L(x0 + d) - lin:11.2e}   {L(x0 + d) - quad:11.2e}")
ts = np.linspace(-3, 3, 121)
plt.figure(figsize=(6.2, 3.4), layout="constrained")
plt.plot(ts, [L(x0 + t * u) for t in ts], "k", lw=2, label="true loss")
plt.plot(ts, L(x0) + ts * (g @ u), "C0--", label="1st-order Taylor")
plt.plot(ts, L(x0) + ts * (g @ u) + 0.5 * ts**2 * (u @ H @ u), "C3", label="2nd-order Taylor")
plt.xlabel("step t along direction u (km)"); plt.ylabel("total ride-km")
plt.title("Taylor approximations near x0 = (4, 3)"); plt.legend(fontsize=8); plt.show()

J analytic:
 [[ 0.7071  0.7071]
 [-0.9487  0.3162]
 [-0.3162 -0.9487]]
max |J_num - J|: 1.8e-10
L(x0) = 306.212  grad = [-32.546  18.051]
   t    true L  1st-order err  2nd-order err
 2.0   330.897      3.49e+01      4.71e+00
 1.0   309.351      8.23e+00      6.88e-01
 0.5   305.645      1.98e+00      9.20e-02
0.25   305.423      4.83e-01      1.18e-02

৭ · Logistic regression — gradient আর Hessian matrix আকারে, তারপর Newton

AI-এর ভিত্তি কোর্সের পর্ব ৫৯-এ সূত্রটি ছিল: logistic regression-এর log loss-এর gradient Xᵀ(σ(Xw) − y)/n. এবার সেটি chain rule দিয়ে derive করি, আর এক ধাপ এগিয়ে Hessian-ও। এখানে loss গড় নয়, যোগফল (n দিয়ে ভাগ করলে gradient আর Hessian দুটিই 1/n গুণ হয়, Newton-এর ধাপ বদলায় না)। প্রতিটি উদাহরণের জন্য z_i = x_iᵀw, p_i = σ(z_i), আর

L(w) = ∑_i [ log(1 + e^z_i) − y_iz_i ]

এটি −[y log p + (1 − y) log(1 − p)]-এরই অন্য রূপ, আর কোডে np.logaddexp(0, z) দিয়ে লিখলে বড় z-এ overflow হয় না। এখন chain rule: log(1 + e^z)-এর z-এর সাপেক্ষে derivative e^z/(1 + e^z) = σ(z), তাই ∂L/∂z_i = p_i − y_i. সব উদাহরণ একসঙ্গে z = Xw, যার Jacobian হলো n × d আকারের X নিজেই। অংশ ৩-এর নিয়মে:

∇L = Xᵀ(p − y), আকার d × 1

Hessian-এর জন্য এই gradient-এর আবার Jacobian নিই। p − y-এর z-এর সাপেক্ষে Jacobian একটি কর্ণ-matrix, কারণ p_i শুধু z_i-এর ওপর নির্ভর করে, আর sigmoid-এর derivative σ′(z) = σ(z)(1 − σ(z)). তাই S = diag(p₁(1 − p₁), …, p_n(1 − p_n)) লিখলে:

∇²L = Xᵀ S X, আকার d × d

S-এর প্রতিটি ঘর ০ থেকে ০.২৫-এর মধ্যে, ধনাত্মক। তাই যেকোনো d-তে dᵀXᵀSXd = ∑ s_i(x_iᵀd)² ≥ 0: Hessian positive semidefinite, log loss convex, কোনো ভুয়া local minimum নেই। এই Hessian দিয়ে Newton-এর ধাপ: Taylor-এর দ্বিতীয় ক্রমের approximation-কে সর্বনিম্ন করা। f(w) + gᵀd + ½dᵀHd-এর d-এর সাপেক্ষে gradient g + Hd; শূন্য বসালে:

H d = −g, তারপর w ← w + d

এই পুনরাবৃত্তির আরেকটি নাম আছে: IRLS, পুরো নাম iteratively reweighted least squares; অর্থাৎ প্রতি ধাপে ওজন বদলে নতুন করে least squares করা। Hd = −g লিখলে প্রতিটি ধাপ হয় একটি weighted least squares, ওজন S-এর ঘরগুলো। ১৯৭২ সালে Nelder আর Wedderburn-এর generalized linear model-এর কাঠামোয় এই পদ্ধতি প্রচলিত হয়, আর statsmodels-এর মতো পরিসংখ্যান library-তে logistic regression আজও সাধারণত এভাবে fit হয় (GLM আসবে পর্ব ২০-এ)। এবার বড় ডেটায় দুটি পদ্ধতি পাশাপাশি: ২০০টি কাল্পনিক অর্ডার, intercept আর দুটি feature. প্রথমে gradient আর Hessian সংখ্যা দিয়ে যাচাই, তারপর ছয় ধাপ Newton বনাম ছয় ধাপ gradient descent.

import numpy as np

rng = np.random.default_rng(42)          # কাল্পনিক বাজারঘর: ২০০টি অর্ডার, ফেরত এল কি না
n = 200
X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))])   # intercept + ২টি feature
y = (rng.random(n) < 1 / (1 + np.exp(-(X @ [-1.0, 2.0, -1.5])))).astype(float)

sig = lambda z: 1 / (1 + np.exp(-z))
loss = lambda w: np.sum(np.logaddexp(0, X @ w) - y * (X @ w))  # log loss-এর যোগফল
grad = lambda w: X.T @ (sig(X @ w) - y)
hess = lambda w: X.T @ (X * (sig(X @ w) * (1 - sig(X @ w)))[:, None])   # XᵀSX

w0, h = np.array([0.3, -0.2, 0.5]), 1e-5
E = np.eye(3) * h
g_num = np.array([(loss(w0 + e) - loss(w0 - e)) / (2 * h) for e in E])
H_num = np.array([(grad(w0 + e) - grad(w0 - e)) / (2 * h) for e in E])
print("gradient check:", f"{np.max(np.abs(g_num - grad(w0))):.1e}",
      "| Hessian check:", f"{np.max(np.abs(H_num - hess(w0))):.1e}")
print("Hessian eigenvalues:", np.round(np.linalg.eigvalsh(hess(w0)), 2))

w_nt, w_gd = np.zeros(3), np.zeros(3)
eta = 4 / np.linalg.eigvalsh(X.T @ X).max()        # GD-র নিরাপদ ধাপ: 1 / (‖X‖² / 4)
print("step  Newton loss     GD loss")
for t in range(1, 7):
    w_nt = w_nt - np.linalg.solve(hess(w_nt), grad(w_nt))   # Newton
    w_gd = w_gd - eta * grad(w_gd)                           # gradient descent
    print(f"{t:>4}  {loss(w_nt):.10f}  {loss(w_gd):.6f}")
best, k = loss(w_nt), 0
while loss(w_gd) - best > 1e-6:
    w_gd, k = w_gd - eta * grad(w_gd), k + 1
print("GD needs", 6 + k, "steps to get within 1e-6 of Newton")
print("w (Newton):", np.round(w_nt, 4))

gradient check: 2.5e-09 | Hessian check: 2.0e-09
Hessian eigenvalues: [33.19 45.82 47.73]
step  Newton loss     GD loss
   1  89.1090771737  92.705837
   2  82.1390890364  86.867112
   3  81.3070395977  84.485820
   4  81.2887610069  83.265540
   5  81.2887496375  82.567522
   6  81.2887496375  82.141213
GD needs 50 steps to get within 1e-6 of Newton
w (Newton): [-1.0208  1.9266 -1.5297]

৮ · Automatic differentiation — chain rule-কে যন্ত্রের হাতে দেওয়া

এখন পর্যন্ত প্রতিটি gradient আমরা হাতে derive করে সংখ্যা দিয়ে মিলিয়েছি। কোটি parameter-এর network-এ হাতে derive করা অবাস্তব, আর finite difference-এ প্রতিটি parameter-এর জন্য দুবার পুরো model চালাতে হয়। তৃতীয় পথ automatic differentiation (autodiff): প্রোগ্রাম চলার সময় প্রতিটি ছোট হিসাব (গুণ, যোগ, exp, log) আর তার নিজস্ব ছোট Jacobian মনে রাখা, তারপর chain rule দিয়ে সেগুলো জোড়া। ফল সূত্রের মতোই নির্ভুল (শুধু float-এর rounding), finite difference-এর মতো আনুমানিক নয়। AI-এর ভিত্তি কোর্সের পর্ব ৬৮-এ PyTorch-এর autograd হাতে লেখা gradient-এর সঙ্গে মিলেছিল; এখানে প্রশ্ন, ভেতরে জোড়ার ক্রম কোনটি।

একটি হিসাব যদি x → h₁ → h₂ → … → L ধাপে চলে, তাহলে J_L = J_k ⋯ J₂J₁. গুণের ক্রম দুই রকম হতে পারে:

ML-এ loss একটি সংখ্যা আর parameter অসংখ্য, তাই reverse mode-ই স্বাভাবিক পছন্দ। অংশ ৩-এর ছবির লাল তীরগুলো ঠিক এটিই: Aᵀ(2r) মানে vector 2r-কে Jacobian A-এর বাঁ দিক থেকে গুণ। এজন্য পুরো Jacobian (যা বিশাল হতে পারে) কখনো বানাতে হয় না। ইতিহাসে forward mode-এর একটি প্রথম দিকের বর্ণনা R. E. Wengert-এর ১৯৬৪ সালের ছোট প্রবন্ধে (Communications of the ACM) পাওয়া যায়; reverse mode-এর আধুনিক রূপের কৃতিত্ব সাধারণত Seppo Linnainmaa-র ১৯৭০ সালের master's thesis-কে দেওয়া হয়, আর neural network-এ backpropagation নামে তা জনপ্রিয় হয় ১৯৮৬ সালে Rumelhart, Hinton আর Williams-এর Nature-এর প্রবন্ধে। কে প্রথম, তা নিয়ে গবেষকদের মধ্যে মতভেদ আছে, কারণ একই ধারণা নিয়ন্ত্রণ-তত্ত্ব আর সংখ্যাতত্ত্বে আলাদাভাবে এসেছিল। Neural network-এ autodiff-এর বিস্তারিত প্রয়োগ পরে Deep Learning কোর্সে আসবে।

৯ · AI-এর যুগে — AI-এর derive করা gradient আর optimizer কোড যাচাই

AI assistant gradient-এর সূত্র আত্মবিশ্বাসী আর পরিষ্কার ভাষায় লেখে। কিন্তু derivation দেখতে নির্ভুল হলেও একটি transpose, একটি ২-এর গুণক বা একটি “A প্রতিসম ধরে নিলাম” প্রায়ই চুপচাপ ঢুকে যায়। ভালো খবর: এই পর্বের প্রতিটি দাবি পাঁচ লাইনের কোডে যাচাই করা যায়।

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

• প্রতিটি পদের পাশে আকার লিখেছি, আর gradient-এর আকার চলকের আকারের সমান কি না দেখেছি
• Float64-এ central difference দিয়ে সূত্রটি অন্তত একটি এলোমেলো, অপ্রতিসম আর অ-বর্গাকার উদাহরণে মিলিয়েছি
• সবচেয়ে ছোট উদাহরণে (2 × 2) হাতে একবার derive করে AI-এর সূত্রের সঙ্গে মিলিয়েছি
• Library-র objective (গড় না যোগফল, ½ আছে কি না, intercept-এ শাস্তি আছে কি না) নিজের version-এর docs থেকে পড়েছি
• Newton বা closed form চালানোর আগে Hessian-এর eigenvalue আর condition number দেখেছি
• ডেটা আগে ভাগ করেছি, test set তালা দেওয়া, আর ফল কয়েকটি seed-এ ও একটি সরল baseline-এর সঙ্গে তুলনা করেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি