Second-order পদ্ধতি — Newton, L-BFGS আর line search

AI-এর ভিত্তি কোর্সের পর্ব ৩৪-এ gradient descent আর learning rate দেখেছেন, আর পর্ব ৪৪-এ scipy.optimize.minimize-কে method না বলে দিলে সে BFGS চালায়, সেটিও দেখেছেন। এই কোর্সের পর্ব ১১-এ logistic regression-এর Hessian XᵀSX বের করে Newton-এর কয়েকটি ধাপ চালিয়েছিলাম, আর দেখেছিলাম ৫০ ধাপের gradient descent-এর কাজ Newton চার-পাঁচ ধাপে সারে। এই পর্বের প্রশ্ন তার পরের: Newton কেন এত দ্রুত, কখন সে ভুল দিকে যায়, তাকে কীভাবে নিরাপদ করা হয়, আর Hessian বানানো যখন অসম্ভব তখন L-BFGS কীভাবে প্রায় একই সুবিধা দেয়? Gradient descent-এর নিজের পরিবার (batch, SGD, schedule) পর্ব ১৩ (Gradient descent-এর পরিবার — batch, SGD, mini-batch আর learning-rate schedule)-এর বিষয়, আর momentum ও Adam পর্ব ১৪ (Momentum, RMSProp আর Adam — derive করে আর তুলনা করে)-এর। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় দশ সেকেন্ড।

১ · ঢাল দেখে হাঁটা বনাম বাঁক দেখে লাফ — কেন second-order

ML-এ প্রায় সব শেখা একটি loss L(w) কমানো। Gradient descent-এর ধাপ w ← w − η∇L, যেখানে η হলো learning rate. এর দুটি দুর্বলতা আছে। প্রথমত, η নিজে বাছতে হয়: বড় হলে loss লাফিয়ে বাড়ে, ছোট হলে হাজার হাজার ধাপ লাগে। দ্বিতীয়ত, loss-এর বাটি এক দিকে সরু আর অন্য দিকে চওড়া হলে, অর্থাৎ Hessian-এর condition number κ = λ_max/λ_min বড় হলে, gradient-এর দিক তলার দিকে নয়, খাড়া দেয়ালের দিকে হেলে থাকে; তাই পথ আঁকাবাঁকা হয়। এই আঁকাবাঁকা পথের ছবি আর হার পর্ব ১৩-এ আসবে। Second-order পদ্ধতি দুটি সমস্যাই একসঙ্গে ধরে: Hessian দিয়ে প্রতিটি দিকের ধাপ আলাদা মাপে বাঁকিয়ে নেয়, আর সঠিক বাঁক জানা থাকলে ধাপের স্বাভাবিক দৈর্ঘ্য ১, কোনো η লাগে না।

টেবিলের “গতি” কথাটি একটু পরে নিখুঁতভাবে সংজ্ঞায়িত হবে। আপাতত মনে রাখুন: মেমরি আর প্রতি ধাপের খরচ যত বেশি, ধাপের সংখ্যা তত কম। কোনটি জিতবে তা নির্ভর করে d কত বড়, ডেটা কত noisy, আর loss কতটা মসৃণ, তার ওপর।

২ · Newton-এর ধাপ — Taylor-এর বাটির তলায় লাফ

পর্ব ১১-এ দেখেছেন, x-এর কাছে মসৃণ f-এর দ্বিতীয় ক্রমের Taylor approximation:

f(x + d) ≈ m(d) = f(x) + gᵀd + ½ dᵀHd

এখানে g = ∇f(x) হলো gradient, H = ∇²f(x) হলো Hessian, আর d হলো সরণ। m(d)-কে বলা হয় quadratic model: আসল f-এর জায়গায় একটি বাটি, যা x-এর কাছে f-এর মান, ঢাল আর বাঁক তিনটিই মেলায়। H positive definite হলে বাটিটির একটিই তলা। m-এর d-এর সাপেক্ষে gradient g + Hd; শূন্য বসালে Newton-এর ধাপ:

H d = −g, তারপর x ← x + d   (এক চলকে: x ← x − f′(x)/f″(x))

কোডে কখনো H-এর inverse বানানো হয় না; সমীকরণটি সমাধান করা হয় (np.linalg.solve), আর H positive definite জানা থাকলে Cholesky দিয়ে, যা পর্ব ১০-এ দেখেছেন। Quadratic function-এ m আর f হুবহু এক, তাই Newton এক ধাপেই সর্বনিম্নে পৌঁছায়; পর্ব ১১-এর ছবি ৬.১-এ লাল তীর ঠিক এটিই দেখিয়েছিল। অ-quadratic function-এ মডেলটি শুধু স্থানীয়ভাবে ঠিক, তাই ধাপ বারবার নিতে হয়। নিচের stepper একটি অ-quadratic এক-চলকের function-এ প্রথম দুটি ধাপ হাতে করে দেখায়।

Newton-এর একটি সুন্দর গুণ হলো affine invariance: চলক x-কে একটি invertible matrix দিয়ে বদলে (যেমন একটি feature-কে টাকা থেকে হাজার টাকায়) নতুন চলকে Newton চালালে ঠিক একই বিন্দুগুলো পাওয়া যায়, শুধু নতুন এককে লেখা। কারণ gradient আর Hessian দুটিই বদলে যায়, আর H⁻¹g-তে সেই বদল কেটে যায়। Gradient descent-এ এমন হয় না; সেখানে feature scaling বদলালে পথ আর গতি দুটিই বদলায়। এজন্য Newton-এর পরিবারের পদ্ধতি ill-conditioned সমস্যায় অনেক কম ভোগে। তবে float-এর rounding আর সংখ্যাগত স্থিতিশীলতার কারণে বাস্তবে scaling তবুও সাহায্য করে।

৩ · Quadratic convergence আর বিশুদ্ধ Newton-এর দুই ব্যর্থতা

গতির ভাষা নিখুঁত করা যাক। e_k = ‖x_k − x*‖ লিখুন, অর্থাৎ k-তম ধাপে সর্বনিম্ন থেকে দূরত্ব। পুনরাবৃত্তিটি converge করছে ধরে:

Newton-এর quadratic convergence-এর তিনটি শর্ত আছে: সর্বনিম্নে Hessian positive definite, Hessian যথেষ্ট মসৃণভাবে বদলায় (Lipschitz), আর শুরুর বিন্দু সর্বনিম্নের যথেষ্ট কাছে। শেষ শর্তটিই সবচেয়ে বিপজ্জনক, কারণ “যথেষ্ট কাছে” কত, আগে থেকে সাধারণত জানা যায় না। নিচের ব্লক তিনটি function-এ বিশুদ্ধ Newton চালায়: একটিতে সব ঠিক, একটিতে বাঁক ভুল চিহ্নের, আর একটিতে শুরু অনেক দূরে। শেষে একটি সরল প্রতিকার, damped Newton: ধাপ যতক্ষণ f কমায় না, ততক্ষণ তাকে অর্ধেক করা।

import numpy as np

def newton(fp, fpp, x, steps):              # 1-D Newton: x ← x − f′(x)/f″(x)
    xs = [x]
    for _ in range(steps):
        x = x - fp(x) / fpp(x)
        xs.append(x)
    return xs

# ১) মসৃণ, convex: f(x) = eˣ − 2x, সর্বনিম্ন x* = ln 2
xs = newton(lambda x: np.exp(x) - 2, np.exp, 0.0, 5)
print(" k   x_k                 error")
for k, x in enumerate(xs):
    print(f"{k:2d}  {x:.15f}  {abs(x - np.log(2)):.1e}")

# ২) f(x) = x⁴ − 2x², f″ < 0 যখন |x| < 0.577: Newton গর্ত নয়, চূড়া খোঁজে
fp, fpp = lambda x: 4*x**3 - 4*x, lambda x: 12*x**2 - 4
print("x^4-2x^2, x0=0.3:", ", ".join(f"{v:.2g}" for v in newton(fp, fpp, 0.3, 4)))
print("  f''(0) =", fpp(0.0), "-> x = 0 is a local MAXIMUM (minima at +-1)")

# ৩) f(x) = √(1 + x²): Newton-এর ধাপ x ← −x³, |x0| > 1 হলে ছিটকে যায়
fp = lambda x: x / np.sqrt(1 + x**2)
fpp = lambda x: (1 + x**2) ** -1.5
print("sqrt(1+x^2), x0=0.5:", ", ".join(f"{v:.3g}" for v in newton(fp, fpp, 0.5, 4)))
print("sqrt(1+x^2), x0=1.5:", ", ".join(f"{v:.3g}" for v in newton(fp, fpp, 1.5, 4)))

# প্রতিকার: damped Newton, ধাপ t অর্ধেক করতে থাকুন যতক্ষণ f না কমে
f, x = lambda x: np.sqrt(1 + x**2), 1.5
for k in range(6):
    d, t = -fp(x) / fpp(x), 1.0
    while f(x + t * d) >= f(x):
        t /= 2
    x = x + t * d
    print(f"damped step {k+1}: t = {t:<9} x = {x:.3e}")

 k   x_k                 error
 0  0.000000000000000  6.9e-01
 1  1.000000000000000  3.1e-01
 2  0.735758882342885  4.3e-02
 3  0.694042299918915  9.0e-04
 4  0.693147581059771  4.0e-07
 5  0.693147180560025  8.0e-14
x^4-2x^2, x0=0.3: 0.3, -0.074, 0.00082, -1.1e-09, 0
  f''(0) = -4.0 -> x = 0 is a local MAXIMUM (minima at +-1)
sqrt(1+x^2), x0=0.5: 0.5, -0.125, 0.00195, -7.45e-09, 0
sqrt(1+x^2), x0=1.5: 1.5, -3.38, 38.4, -5.68e+04, 1.83e+14
damped step 1: t = 0.5       x = -9.375e-01
damped step 2: t = 1.0       x = 8.240e-01
damped step 3: t = 1.0       x = -5.594e-01
damped step 4: t = 1.0       x = 1.751e-01
damped step 5: t = 1.0       x = -5.366e-03
damped step 6: t = 1.0       x = 1.545e-07

খরচের হিসাব — কেন বড় model-এ বিশুদ্ধ Newton অচল

d সংখ্যক parameter-এ Hessian-এ d² ঘর, আর Cholesky দিয়ে সমাধানে প্রায় d³/3 গুণ-যোগ। float64-এ প্রতিটি ঘর ৮ বাইট। d = ১,০০০ হলে Hessian ৮ মেগাবাইট আর সমাধানে প্রায় ৩.৩ × ১০^৮ গুণ-যোগ, একটি ল্যাপটপে এক সেকেন্ডেরও কম। d = ১০ লাখ (একটি TF-IDF ফিচারের logistic regression-এ খুব সাধারণ) হলে Hessian ৮ × ১০^১২ বাইট, অর্থাৎ ৮ টেরাবাইট: শুধু রাখাই অসম্ভব। তার ওপর Hessian নিজে হিসাব করতে পুরো ডেটার ওপর d² মাপের কাজ। এই দুই বাধা, ভুল বাঁক আর বিশাল খরচ, এই পর্বের বাকি অংশের দুটি সুতো: line search প্রথমটির উত্তর, quasi-Newton দ্বিতীয়টির।

৪ · Line search — Armijo-র “যথেষ্ট কমা”র শর্ত

Damped Newton-এর ধারণা সাধারণ: দিক d ঠিক রেখে ধাপের দৈর্ঘ্য t আলাদা করে বাছা, x ← x + t·d. দিক বাছা আর দূরত্ব বাছা দুটি আলাদা কাজ। দিক যে কোনো পদ্ধতির হতে পারে (gradient-এর উল্টো, Newton, BFGS), শর্ত শুধু একটি: d যেন descent direction হয়, অর্থাৎ gᵀd < ০. তাহলে খুব ছোট t-তে f অবশ্যই কমে। দূরত্ব বাছার কাজটির নাম line search, কারণ x থেকে d-এর দিকে একটি সরলরেখার ওপর এক চলকের function φ(t) = f(x + t·d) দেখা হয়।

শুধু “f কমেছে কি না” দেখা যথেষ্ট নয়: খুব সামান্য কমা বারবার নিলে পদ্ধতি সর্বনিম্নে না পৌঁছে আটকে থাকতে পারে। Larry Armijo ১৯৬৬ সালের একটি প্রবন্ধে যে শর্ত দেন, এখন তাকে বলে Armijo condition বা sufficient decrease:

f(x + t·d) ≤ f(x) + c·t·gᵀd,   ০ < c < ১

ডান দিকটি একটি সরলরেখা, শুরুর ঢাল gᵀd-এর c গুণ। c = ১ হলে রেখাটি φ-এর tangent; ছোট c (প্রচলিত মান ১০^−৪) মানে “ঢাল যা প্রতিশ্রুতি দেয়, তার সামান্য একটি অংশ কমলেই চলবে”। Backtracking হলো এই শর্ত পরীক্ষার সবচেয়ে সরল কৌশল: t = ১ দিয়ে শুরু, শর্ত না মিললে t ← β·t (প্রায়ই β = ০.৫), যতক্ষণ না মেলে। Newton আর quasi-Newton-এ t = ১ দিয়ে শুরু করা জরুরি, কারণ সর্বনিম্নের কাছে পুরো ধাপই quadratic বা superlinear গতি দেয়; line search তখন চুপচাপ t = ১ মেনে নেয়।

Wolfe শর্ত — খুব ছোট ধাপও আটকানো

Armijo খুব বড় ধাপ আটকায়, কিন্তু খুব ছোট ধাপ আটকায় না; backtracking-এ সমস্যা কম, কারণ t বড় থেকে ছোট হয় আর প্রথম গ্রহণযোগ্য মানেই থামে। অন্য কৌশলে (যেমন t বাড়িয়ে খোঁজা) দ্বিতীয় একটি শর্ত লাগে, যাকে বলে curvature condition: φ′(t) ≥ c₂·φ′(০), যেখানে c < c₂ < ১. মানে নতুন বিন্দুতে ঢাল যেন শুরুর মতো খাড়া নিচের দিকে না থাকে; থাকলে আরও এগোনো উচিত ছিল। দুই শর্ত মিলে Wolfe conditions, Philip Wolfe-এর ১৯৬৯ সালের কাজের নামে। Nocedal আর Wright-এর পাঠ্যবইয়ে Newton ও quasi-Newton-এর জন্য c₂ = ০.৯ প্রচলিত বলা হয়েছে। উপরের stepper-এ t = ০.১২৫-এ φ′ = −১০১ + ১০০১ × ০.১২৫ = ২৪.১, যা −৯০.৯-এর চেয়ে বড়, তাই Wolfe-ও মেলে। BFGS-এর জন্য এই শর্ত বিশেষ গুরুত্বপূর্ণ, অংশ ৬-এ কারণ দেখবেন।

৫ · Rosenbrock-এর কলার মতো উপত্যকা — fixed step, Armijo আর damped Newton

Optimizer পরীক্ষার সবচেয়ে পরিচিত function, Rosenbrock function, Howard Rosenbrock-এর ১৯৬০ সালের প্রবন্ধ থেকে: f(x, y) = (1 − x)² + 100(y − x²)². সর্বনিম্ন (১, ১)-এ, মান ০. সমস্যা হলো তলাটি একটি সরু, বাঁকা, কলার মতো উপত্যকা (y ≈ x²): দেয়াল খাড়া (১০০ গুণক), কিন্তু উপত্যকার মেঝে প্রায় সমতল। সর্বনিম্নে Hessian-এর condition number প্রায় ২,৫০০. প্রচলিত শুরু (−১.২, ১). নিচের ব্লক তিনটি পদ্ধতি চালায়, আর গোনে কতবার f হিসাব করা হলো, কারণ বাস্তবে প্রতিটি f-হিসাব মানে পুরো ডেটার ওপর একবার loss.

import numpy as np
import matplotlib.pyplot as plt

calls = 0
def f(p):                                   # Rosenbrock, সর্বনিম্ন (1, 1)-এ, মান 0
    global calls; calls += 1
    return (1 - p[0])**2 + 100 * (p[1] - p[0]**2)**2
grad = lambda p: np.array([-2*(1 - p[0]) - 400*p[0]*(p[1] - p[0]**2), 200*(p[1] - p[0]**2)])
hess = lambda p: np.array([[2 - 400*p[1] + 1200*p[0]**2, -400*p[0]], [-400*p[0], 200.0]])

def armijo(p, d, g, t=1.0, c=1e-4, beta=0.5):   # backtracking: শর্ত না মেলা পর্যন্ত t অর্ধেক
    fp = f(p)
    while f(p + t * d) > fp + c * t * (g @ d) and t > 1e-10:
        t *= beta
    return t
def run(name, direction, use_ls, iters, eta=1e-3):
    global calls; calls, p = 0, np.array([-1.2, 1.0]); path = [p]
    for _ in range(iters):
        g = grad(p); d = direction(p, g)
        p = p + (armijo(p, d, g) if use_ls else eta) * d
        path.append(p)
    print(f"{name:<18} iters {iters:4d}  f-evals {calls:5d}  f = {f(p):.1e}  x = {np.round(p, 4)}")
    return np.array(path)
sd = lambda p, g: -g                                   # steepest descent
nt = lambda p, g: -np.linalg.solve(hess(p), g)          # Newton-এর দিক
P1 = run("GD, fixed 1e-3", sd, False, 2000)
P2 = run("GD + Armijo", sd, True, 2000)
P3 = run("Newton + Armijo", nt, True, 20)
X, Y = np.meshgrid(np.linspace(-1.5, 1.5, 200), np.linspace(-0.5, 1.6, 200))
plt.figure(figsize=(6.4, 4), layout="constrained")
plt.contour(X, Y, (1 - X)**2 + 100 * (Y - X**2)**2, levels=np.logspace(-1, 3, 14), colors="0.7")
for P, s, lab in [(P1, "C0", "GD fixed step"), (P2, "C1", "GD + Armijo"), (P3, "C3.-", "Newton + Armijo")]:
    plt.plot(P[:, 0], P[:, 1], s, lw=1.5, ms=4, label=lab)
plt.plot(1, 1, "k*", ms=12); plt.legend(fontsize=8); plt.title("Rosenbrock: paths from (-1.2, 1)"); plt.show()

GD, fixed 1e-3     iters 2000  f-evals     0  f = 7.8e-02  x = [0.7207 0.518 ]
GD + Armijo        iters 2000  f-evals 21742  f = 1.3e-04  x = [0.9887 0.9774]
Newton + Armijo    iters   20  f-evals    47  f = 8.5e-12  x = [1. 1.]

৬ · Quasi-Newton আর BFGS — gradient-এর পার্থক্য থেকে বাঁক শেখা

Newton-এর দামি অংশ Hessian. কিন্তু প্রতিটি ধাপে আমরা একটি সস্তা তথ্য পাই, যা এতক্ষণ ফেলে দিচ্ছিলাম: পরপর দুটি বিন্দুর gradient. এক চলকে দুটি বিন্দুর ঢালের পার্থক্যকে দূরত্ব দিয়ে ভাগ করলেই বাঁকের আন্দাজ: f″ ≈ (f′(x_k+1) − f′(x_k)) / (x_k+1 − x_k). এটি পুরোনো secant পদ্ধতির ভাবনা। বহু চলকে লিখুন

s_k = x_k+1 − x_k (কতটা সরলাম),   y_k = g_k+1 − g_k (gradient কতটা বদলাল)

Quadratic function-এ y = Hs হুবহু, অন্যথায় প্রায়। Quasi-Newton পদ্ধতিগুলো inverse Hessian-এর একটি আন্দাজ M রাখে (শুরুতে প্রায়ই M₀ = I), দিক নেয় d = −M·g, আর প্রতি ধাপে M-কে এমনভাবে হালনাগাদ করে যাতে সে সর্বশেষ জোড়াটি হুবহু মানে:

M_k+1 y_k = s_k   (secant condition)

এক চলকে এই শর্তই M পুরোপুরি ঠিক করে দেয়। বহু চলকে অসংখ্য M এই শর্ত মানে; কোনটি বাছবেন? BFGS-এর উত্তর: যেটি আগের M_k-এর সবচেয়ে কাছে (একটি নির্দিষ্ট ওজনদার মাপে), প্রতিসম আর positive definite. ফল একটি rank-২ হালনাগাদ, ρ = ১/(yᵀs) লিখে:

M_k+1 = (I − ρ s yᵀ) M_k (I − ρ y sᵀ) + ρ s sᵀ

নামটি চারজনের: Charles Broyden, Roger Fletcher, Donald Goldfarb আর David Shanno, যাঁরা ১৯৭০ সালে আলাদা আলাদা প্রবন্ধে প্রায় একই সূত্রে পৌঁছান। এর আগের DFP পদ্ধতি (William Davidon-এর ১৯৫৯ সালের একটি প্রতিবেদন, পরে Fletcher আর Powell ১৯৬৩ সালে জনপ্রিয় করেন) একই ধারণার আরেকটি রূপ; বাস্তবে BFGS বেশি স্থিতিশীল বলে প্রমাণিত হয়েছে। একটি জরুরি শর্ত: yᵀs > ০ হতে হবে, নইলে ρ ঋণাত্মক আর M positive definite থাকে না। Convex function-এ এটি সবসময় সত্য; non-convex-এ Wolfe-এর curvature শর্ত এটি নিশ্চিত করে, কারণ φ′(t) ≥ c₂φ′(০) থেকেই yᵀs > ০ আসে। এজন্যই BFGS-এর সঙ্গে Wolfe line search জোড়া হয়।

নিচের ব্লক পর্ব ১১-এর বাটি f = ½xᵀAx, A = [[২.৫, ১.৫], [১.৫, ২.৫]]-এ BFGS নিজে লিখে চালায়। A কোথাও M-এর হালনাগাদে ঢোকে না; সব বাঁকের তথ্য আসে শুধু gradient-এর পার্থক্য y থেকে। দুইবার: একবার নিখুঁত line search দিয়ে (quadratic-এ সেরা t বের করা যায় শুধু gradient দিয়ে, t = −gᵀd / dᵀ(∇f(x + d) − g)), একবার Armijo দিয়ে। প্রতি ধাপে ছাপা হয় M আসল A⁻¹ থেকে কত দূরে, আর secant শর্ত কতটা মানা হলো।

import numpy as np
np.set_printoptions(precision=4, suppress=True)

A = np.array([[2.5, 1.5], [1.5, 2.5]])      # আসল Hessian; BFGS এটি কখনো দেখে না
f = lambda x: 0.5 * x @ A @ x
grad = lambda x: A @ x
A_inv = np.linalg.inv(A)                    # শুধু তুলনার জন্য

def bfgs_update(M, s, y):                   # inverse Hessian-এর আন্দাজ M-এর হালনাগাদ
    rho, I = 1.0 / (y @ s), np.eye(len(s))
    return (I - rho * np.outer(s, y)) @ M @ (I - rho * np.outer(y, s)) + rho * np.outer(s, s)

def run(exact, steps):
    x, M = np.array([2.0, -0.5]), np.eye(2)  # M0 = I: প্রথম ধাপটি সাধারণ gradient descent
    for k in range(1, steps + 1):
        g = grad(x); d = -M @ g
        if exact:                             # quadratic-এ নিখুঁত line search, শুধু gradient দিয়ে
            t = -(g @ d) / (d @ (grad(x + d) - g))
        else:
            t = 1.0                           # Armijo backtracking
            while f(x + t * d) > f(x) + 1e-4 * t * (g @ d):
                t *= 0.5
        s = t * d; y = grad(x + s) - g        # বাঁকের তথ্য শুধু gradient-এর পার্থক্যে
        M = bfgs_update(M, s, y); x = x + s
        print(f"step {k}: t = {t:.3f}  |M - A^-1| = {np.linalg.norm(M - A_inv):.1e}"
              f"  |My - s| = {np.linalg.norm(M @ y - s):.0e}  f = {f(x):.1e}")
    return M

print("exact line search:"); M = run(True, 2)
print("M after 2 steps:\n", M, "\ntrue A^-1:\n", A_inv)
print("Armijo line search:"); run(False, 6)

exact line search:
step 1: t = 0.281  |M - A^-1| = 1.2e-01  |My - s| = 4e-16  f = 8.4e-01
step 2: t = 0.889  |M - A^-1| = 1.6e-16  |My - s| = 2e-16  f = 1.2e-32
M after 2 steps:
 [[ 0.625 -0.375]
 [-0.375  0.625]] 
true A^-1:
 [[ 0.625 -0.375]
 [-0.375  0.625]]
Armijo line search:
step 1: t = 0.500  |M - A^-1| = 1.2e-01  |My - s| = 4e-16  f = 2.6e+00
step 2: t = 1.000  |M - A^-1| = 6.7e-02  |My - s| = 7e-16  f = 1.3e-02
step 3: t = 1.000  |M - A^-1| = 7.2e-03  |My - s| = 3e-17  f = 5.1e-05
step 4: t = 1.000  |M - A^-1| = 6.9e-03  |My - s| = 0e+00  f = 1.2e-08
step 5: t = 1.000  |M - A^-1| = 5.8e-03  |My - s| = 1e-20  f = 2.3e-12
step 6: t = 1.000  |M - A^-1| = 3.3e-03  |My - s| = 0e+00  f = 2.7e-17

৭ · L-BFGS — শুধু শেষ m জোড়া মনে রাখা

BFGS Hessian হিসাবের খরচ বাঁচায়, কিন্তু মেমরি নয়: M তবুও d × d ঘন matrix, d = ১০ লাখে ৮ টেরাবাইট। Jorge Nocedal ১৯৮০ সালের একটি প্রবন্ধে সমাধান দেন: M জমা না রেখে শুধু শেষ m জোড়া (s_i, y_i) রাখুন, আর প্রতিবার দরকারি গুণফল M·g সরাসরি সেই জোড়াগুলো থেকে হিসাব করুন। এর নাম L-BFGS (L মানে limited memory, অর্থাৎ সীমিত স্মৃতি)। Liu আর Nocedal-এর ১৯৮৯ সালের প্রবন্ধ একে বড় সমস্যায় জনপ্রিয় করে। হিসাবের কৌশলটির নাম two-loop recursion: প্রথম loop নতুন থেকে পুরোনো জোড়ার দিকে চলে, মাঝে একটি সরল শুরুর আন্দাজ γI, তারপর দ্বিতীয় loop পুরোনো থেকে নতুনে। খরচ প্রায় 4md গুণ-যোগ, মেমরি 2md সংখ্যা।

নিচের stepper ছোট সংখ্যায় পুরো two-loop দেখায়, m = ২. বাটি f = ½(x² + ১০y²), Hessian diag(১, ১০), বর্তমান বিন্দু (২, ১), gradient g = (২, ১০). স্মৃতিতে দুটি পুরোনো জোড়া: s₁ = (১, −১), y₁ = (১, −১০) আর নতুনটি s₂ = (১, ০), y₂ = (১, ০). (Quadratic বলে y = Hs.)

বাস্তবে m সাধারণত ৩ থেকে ২০; SciPy-র L-BFGS-B-তে এর নাম maxcor, default ১০. m বাড়ালে প্রতি ধাপ একটু দামি হয় আর দিক একটু ভালো হয়, কিন্তু লাভ দ্রুত কমে যায়। L-BFGS-B-র “B” মানে bounds: প্রতিটি চলকের নিচের আর উপরের সীমা (যেমন ওজন ≥ ০) মেনে চলার ক্ষমতা, Byrd, Lu, Nocedal আর Zhu-র ১৯৯৫ সালের কাজ থেকে। লেখার সময়ে SciPy-র সোর্সে লেখা আছে, তারা L-BFGS-B v3.0-র মূল Fortran কোডের একটি C-তে অনূদিত ও পরিবর্তিত রূপ ব্যবহার করে। সীমাসহ optimization-এর তত্ত্ব পর্ব ১৬ (Constraint-সহ optimization — Lagrange multiplier আর KKT)-এর বিষয়।

৮ · SciPy আর scikit-learn-এ — চারটি পদ্ধতি পাশাপাশি

নিজে লেখা optimizer শেখার জন্য ভালো, কাজের জন্য library. scipy.optimize.minimize-এ এই পর্বের প্রতিটি ধারণার একটি তৈরি রূপ আছে:

Newton-CG-এর কৌশলটি বিশেষ করে মনে রাখার মতো। Conjugate gradient পদ্ধতি Hd = −g সমাধান করে H-কে কখনো পুরোটা না দেখে, শুধু বারবার “H গুণ একটি vector v” চেয়ে। এই Hessian-vector product Hv বের করতে পুরো Hessian লাগে না: logistic regression-এ Xᵀ(S(Xv)), দুটি matrix-vector গুণ; আর যেকোনো autodiff-যোগ্য model-এ gradient-এর দিকনির্দেশক derivative হিসেবে, খরচ প্রায় একটি gradient-এর কয়েক গুণ (Pearlmutter, ১৯৯৪)। CG-কে কয়েক ধাপ পরে থামিয়ে দিলে একটি “যথেষ্ট ভালো” Newton-এর দিক পাওয়া যায়।

নিচের ব্লক চারটি method দুটি সমস্যায় চালায়: Rosenbrock (SciPy-র নিজস্ব rosen, rosen_der, rosen_hess) আর বাজারঘরের ৫০০টি কাল্পনিক অর্ডারের ফেরত-সম্ভাবনার logistic regression, ছোট L2 শাস্তিসহ, সূত্রে দেওয়া gradient আর Hessian দিয়ে (পর্ব ১১-এর সূত্র, শুধু শাস্তির λI যোগ)। মাপকাঠি x* বের করা হয় ৩০ ধাপ বিশুদ্ধ Newton দিয়ে। শেষে একই logistic objective scikit-learn-এর তিনটি solver দিয়ে।

import numpy as np
from scipy.optimize import minimize, rosen, rosen_der, rosen_hess
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(42)          # কাল্পনিক বাজারঘর: ৫০০টি অর্ডার, ৫টি feature, ফেরত এল কি না
n, d = 500, 6
X = np.column_stack([np.ones(n), rng.normal(size=(n, d - 1))])
y = (rng.random(n) < 1 / (1 + np.exp(-(X @ rng.normal(size=d))))).astype(float)
sig = lambda z: 1 / (1 + np.exp(-z))
loss = lambda w: np.sum(np.logaddexp(0, X @ w) - y * (X @ w)) + 0.5 * w @ w   # L2, λ = 1
grad = lambda w: X.T @ (sig(X @ w) - y) + w
hess = lambda w: X.T @ (X * (sig(X @ w) * (1 - sig(X @ w)))[:, None]) + np.eye(d)

for name, (f, g, H, x) in {"Rosenbrock": (rosen, rosen_der, rosen_hess, np.array([-1.2, 1.0])),
                           "logistic": (loss, grad, hess, np.zeros(d))}.items():
    ref = x.copy()
    for _ in range(30):                  # তুলনার মাপকাঠি: বিশুদ্ধ Newton, ৩০ ধাপ
        ref = ref - np.linalg.solve(H(ref), g(ref))
    print(f"--- {name}, dimension {x.size}, x* = {np.round(ref, 3)}")
    print("method        nit  nfev  njev  nhev  |x - x*|")
    for m in ["BFGS", "L-BFGS-B", "Newton-CG", "trust-constr"]:
        kw = {"hess": H} if m in ("Newton-CG", "trust-constr") else {}
        r = minimize(f, x, jac=g, method=m, **kw)
        print(f"{m:<12} {r.nit:4d} {r.nfev:5d} {r.njev:5d} {getattr(r, 'nhev', 0):5d}  {np.linalg.norm(r.x - ref):.1e}")

for s in ["lbfgs", "newton-cg", "newton-cholesky"]:     # একই objective, scikit-learn দিয়ে
    clf = LogisticRegression(C=1.0, fit_intercept=False, solver=s, tol=1e-8).fit(X, y)
    print(f"sklearn {s:<16} n_iter {clf.n_iter_[0]:3d}  |w - w*| = {np.linalg.norm(clf.coef_[0] - ref):.1e}")

--- Rosenbrock, dimension 2, x* = [1. 1.]
method        nit  nfev  njev  nhev  |x - x*|
BFGS           32    39    39     0  6.1e-08
L-BFGS-B       36    44    44     0  2.2e-06
Newton-CG      83   105   105    83  3.9e-05
trust-constr   44    44    30    30  3.3e-14
--- logistic, dimension 6, x* = [1.313 0.208 1.189 0.25  1.2   0.693]
method        nit  nfev  njev  nhev  |x - x*|
BFGS           17    24    24     0  6.2e-08
L-BFGS-B        8     9     9     0  5.1e-06
Newton-CG       8     9     9     8  1.7e-09
trust-constr   22    14    13    13  1.2e-09
sklearn lbfgs            n_iter  10  |w - w*| = 3.8e-08
sklearn newton-cg        n_iter   6  |w - w*| = 5.6e-11
sklearn newton-cholesky  n_iter   5  |w - w*| = 1.0e-07

scikit-learn কেন lbfgs দিয়ে শুরু করে

scikit-learn-এর ডকুমেন্টেশন অনুযায়ী LogisticRegression-এর default solver ০.২২ version-এ liblinear থেকে বদলে lbfgs করা হয়েছে, আর docs-এর ভাষায় এটি “a good default solver because it works reasonably well for a wide class of problems”. কারণগুলো এই পর্ব থেকেই বোঝা যায়: L2 শাস্তিসহ log loss মসৃণ আর convex, তাই quasi-Newton-এর জন্য আদর্শ; মেমরি O(md), তাই feature অনেক হলেও চলে; Hessian লাগে না; আর multinomial (বহু শ্রেণির softmax) loss-ও সরাসরি সামলায়। একই docs বলে, নমুনা অনেক আর feature কম হলে (n ≫ d) newton-cholesky ভালো পছন্দ, যা ১.২ version-এ যোগ হয়েছে আর পুরো Hessian বানায়, তাই এর মেমরি d²-এর অনুপাতে বাড়ে। L1 শাস্তি মসৃণ নয় বলে lbfgs তা পারে না; সেখানে liblinear বা saga. ব্যবহারিক সতর্কতা: feature-এর scale খুব আলাদা হলে lbfgs-এর ধাপ বেশি লাগে, আর max_iter (default ১০০)-এর মধ্যে না থামলে scikit-learn একটি ConvergenceWarning দেয়; সমাধান সাধারণত আগে StandardScaler, তারপর প্রয়োজনে max_iter বাড়ানো।

পরিসংখ্যানের library-তে একই logistic regression প্রায়ই অন্য নামে চলে: IRLS (iteratively reweighted least squares). পর্ব ১১-এ দেখেছেন, logistic-এ Newton-এর প্রতিটি ধাপ আসলে S-এর ঘরগুলোকে ওজন ধরে একটি weighted least squares. statsmodels-এর GLM এভাবে fit করে, আর standard error-এর জন্য শেষ Hessian কাজে লাগায়। পুরো exponential family-তে এই ভাবনার সাধারণ রূপ পর্ব ২০ (Exponential family আর GLM — Poisson, binomial আর link function)-এ আসবে।

৯ · কোন পদ্ধতি কখন — সিদ্ধান্তের ছক

১০ · AI-এর যুগে — AI-এর লেখা optimizer কোড যাচাই

“এই loss কমানোর কোড লিখে দাও” বললে AI assistant প্রায়ই চলনসই কোড দেয়, কিন্তু optimizer-এর ভুল চুপচাপ থাকে: কোড চলে, একটি সংখ্যা ছাপে, আর সেটি সর্বনিম্ন নয়। নিচের ভুলগুলো বারবার দেখা যায়।

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

• প্রতি ধাপে loss ছাপিয়ে দেখেছি তা কমছে, আর শেষ gradient-এর মাপ প্রায় শূন্য
• Gradient (আর Hessian থাকলে Hessian) central difference দিয়ে float64-এ মিলিয়েছি
• d থেকে মেমরি হিসাব করেছি: d² না 2md, আর পদ্ধতি সেই অনুযায়ী বেছেছি
• success আর message পড়েছি, আর ConvergenceWarning চাপা না দিয়ে কারণ খুঁজেছি
• দুটি আলাদা method বা কয়েকটি শুরু থেকে একই উত্তর আসে কি না দেখেছি
• Library-র objective (শাস্তির রূপ, intercept, গড় না যোগফল) নিজের version-এর docs থেকে মিলিয়েছি
• ডেটা আগে ভাগ করেছি, test set তালা দেওয়া; optimizer-এর tolerance বা শাস্তি test set দেখে বাছিনি, আর কোনো ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি