Constraint-সহ optimization — Lagrange multiplier আর KKT

AI-এর ভিত্তি কোর্সের পর্ব ৩৪-এ gradient descent শিখেছেন, যেখানে x যেকোনো দিকে যেতে পারত; আর পর্ব ৪৪-এ scipy.optimize-এর কয়েকটি global optimizer দেখেছেন। এই কোর্সের পর্ব ১১-এ gradient আর Hessian matrix আকারে এসেছে, আর পর্ব ০৮-এ PCA-র ভেতরে একবার Lagrange multiplier উঁকি দিয়েছিল। এই পর্বের প্রশ্ন: x যদি স্বাধীন না হয়, যদি বাজেট, সীমা বা “যোগফল ১” জাতীয় শর্ত মানতেই হয়, তাহলে সর্বনিম্ন কোথায় আর কীভাবে চিনবেন? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy আর matplotlib); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে।

১ · Constraint কেন — শর্ত মেনে সেরা সমাধান

Optimization-এর সাধারণ রূপ আগের পর্বগুলোতে ছিল “f(x)-কে যত ছোট সম্ভব করুন”, যেখানে f হলো objective function (ML-এ সাধারণত loss). এখন তার সঙ্গে যোগ হলো constraint, অর্থাৎ শর্ত। প্রচলিত লেখার ধরন (standard form):

সর্বনিম্ন করুন f(x), যেখানে g_i(x) ≤ 0 (i = 1…m) আর h_j(x) = 0 (j = 1…p)

এখানে x ∈ ℝⁿ চলকের vector; প্রতিটি g_i একটি inequality constraint (“≤” শর্ত), আর প্রতিটি h_j একটি equality constraint. “x ≤ ৫”-কে লেখা হয় g(x) = x − 5 ≤ 0, আর “x ≥ ০”-কে g(x) = −x ≤ 0; তাই সব অসমতা এক দিকে আনা যায়। যে সব x সব শর্ত মানে, তাদের মিলিত সেটই feasible set; সমাধান খোঁজা হয় শুধু তার ভেতরে। “সর্বোচ্চ করুন f” মানে “সর্বনিম্ন করুন −f”, তাই একটি রূপই যথেষ্ট।

ML-এ এ ধরনের শর্ত যতটা মনে হয়, তার চেয়ে অনেক বেশি জায়গায় আছে:

সমস্যার পরিবারগুলোর নামও জেনে রাখুন, কারণ টুল বাছাই এর ওপর নির্ভর করে। f আর সব শর্ত রৈখিক হলে linear programming (LP); f দ্বিঘাত আর শর্ত রৈখিক হলে quadratic programming (QP), যেমন SVM; বাকিগুলো সাধারণ nonlinear programming (NLP). f convex আর feasible set-ও convex হলে সমস্যাটি convex; তখন যেকোনো local সমাধানই global. Convexity-র পুরো আলোচনা পর্ব ১২ (Convexity — কেন কিছু loss সহজ আর কিছু কঠিন)-এর বিষয়; এখানে শুধু ফলটি ব্যবহার করব।

২ · সমতার শর্ত আর Lagrange multiplier — রেখাকে ছুঁয়ে যাওয়া contour

সবচেয়ে ছোট উদাহরণ দিয়ে শুরু করা যাক: f(x, y) = x² + 2y²-কে সর্বনিম্ন করুন, শর্ত x + y = 3. শর্তটিকে লিখি h(x, y) = x + y − 3 = 0. শর্ত না থাকলে উত্তর (০, ০), কিন্তু সেটি রেখার ওপর নেই। f-এর contour, অর্থাৎ সমান মানের বক্ররেখাগুলো, এখানে মূলবিন্দু-কেন্দ্রিক উপবৃত্ত। রেখা ধরে হাঁটলে f-এর মান কখনো বাড়ে, কখনো কমে। সর্বনিম্ন সেখানে, যেখানে রেখাটি কোনো একটি উপবৃত্তকে কেটে না গিয়ে শুধু ছুঁয়ে যায়। কেটে গেলে রেখার এক দিকে f আরও ছোট, তাই সেটি সর্বনিম্ন হতে পারে না।

কেন ছোঁয়ার বিন্দুতে দুটি gradient সমান্তরাল? পর্ব ১১-এ দেখেছেন, gradient সবসময় নিজের contour-এর সঙ্গে লম্ব। ∇h লম্ব শর্তের রেখার সঙ্গে, আর ∇f লম্ব f-এর contour-এর সঙ্গে। ছোঁয়ার বিন্দুতে দুটি বক্ররেখার স্পর্শক একই, তাই দুটি লম্বও একই রেখায়। ভিন্নভাবে বললে: রেখা বরাবর একক দিক d-তে f-এর বদলের হার ∇fᵀd. সর্বনিম্নে এটি শূন্য হতেই হবে, নইলে কোনো এক দিকে হাঁটলে f কমত। অর্থাৎ ∇f-এর রেখা বরাবর কোনো অংশ নেই, পুরোটা ∇h-এর দিকে:

∇f(x*) = λ ∇h(x*),   h(x*) = 0

এই সংখ্যা λ-ই Lagrange multiplier. সমীকরণ দুটি একসঙ্গে লেখার একটি কৌশল আছে: Lagrangian নামের নতুন function L(x, λ) = f(x) − λh(x) বানান, আর তার সব চলকে (x আর λ দুটিতেই) ঢাল শূন্য বসান। ∂L/∂x = 0 দেয় ∇f = λ∇h, আর ∂L/∂λ = 0 দেয় h(x) = 0. অর্থাৎ শর্তসহ একটি সমস্যা হয়ে গেল শর্তহীন একটি stationary point খোঁজার সমস্যা, একটি বাড়তি চলক নিয়ে। (সমতার শর্তে λ-এর চিহ্ন যা খুশি হতে পারে, তাই অনেক বই f + λh লেখে; উত্তরে শুধু λ-এর চিহ্ন উল্টায়।)

নিচের stepper রেখা বরাবর চারটি বিন্দুতে দুটি gradient পাশাপাশি রেখে দেখায় কোথায় তারা সমান্তরাল হয়। “রেখা বরাবর অংশ” মানে ∇f-এর সঙ্গে রেখার একক দিক d = (1, −1)/√2-এর dot product; সেটি শূন্য হলেই সর্বনিম্ন।

একাধিক সমতা থাকলে প্রতিটির জন্য আলাদা λ_j, আর শর্ত হয় ∇f = ∑λ_j∇h_j. একটি সূক্ষ্ম কথা: এই নিয়ম কাজ করে যদি সমাধানে শর্তগুলোর gradient-গুলো রৈখিকভাবে স্বাধীন হয় (এর নাম LICQ, একটি constraint qualification). এটি না মানলে সর্বনিম্ন বিন্দুতেও কোনো λ না-ও থাকতে পারে; ভাবনার প্রশ্ন ১-এ একটি উদাহরণ আছে।

৩ · বৃত্তের ওপর quadratic — হাতে Lagrange, তারপর SciPy দিয়ে যাচাই

এবার একটি পরিচিত প্রশ্ন: প্রতিসম A = [[৩, ১], [১, ৩]]-এর জন্য f(x) = xᵀAx-কে সর্বোচ্চ করুন, শর্ত ‖x‖² = xᵀx = 1, অর্থাৎ একক বৃত্তের ওপর। পর্ব ১০-এ এর contour দেখেছেন: A positive definite, তাই contour তির্যক উপবৃত্ত। শর্ত h(x) = xᵀx − 1. পর্ব ১১-এর সূত্রে ∇f = 2Ax (A প্রতিসম বলে) আর ∇h = 2x. Lagrange-এর শর্ত:

2Ax = λ · 2x  ⟹  Ax = λx

অর্থাৎ শর্তসহ stationary point মানেই A-এর eigenvector, আর λ তার eigenvalue. সেখানে f = xᵀAx = λxᵀx = λ, তাই সর্বোচ্চ মান সবচেয়ে বড় eigenvalue. A-এর eigenvalue ৪ (দিক (১, ১)/√২) আর ২ (দিক (১, −১)/√২)। এটিই পর্ব ০৮-এর Rayleigh quotient-এর ফল, আর PCA-র প্রথম principal component-এর যুক্তি: variance uᵀΣu সর্বোচ্চ করুন, শর্ত ‖u‖ = 1. নিচের ব্লক একই প্রশ্ন SciPy-র SLSQP দিয়ে সংখ্যায় সমাধান করে, তারপর তিনটি জিনিস মেলায়: eigh-এর উত্তর, দুটি gradient সমান্তরাল কি না, আর λ. শেষে contour-এর ছবি আঁকে।

SLSQP (Sequential Least SQuares Programming) প্রতিটি ধাপে f-এর একটি দ্বিঘাত আর শর্তের একটি রৈখিক আন্দাজ বানিয়ে ছোট একটি QP সমাধান করে; Dieter Kraft-এর ১৯৮৮ সালের Fortran কোড থেকে এসেছে বলে SciPy-র docs-এ উল্লেখ আছে। SciPy-তে minimize সবসময় ছোট করে, তাই সর্বোচ্চের জন্য −f দেওয়া হয়েছে; আর "eq" ধরনের শর্তের function শূন্য হতে হয়।

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

A = np.array([[3.0, 1.0], [1.0, 3.0]])
f = lambda v: v @ A @ v                         # এটিকে সর্বোচ্চ করতে চাই
circle = {"type": "eq", "fun": lambda v: v @ v - 1}   # h(x) = xᵀx − 1 = 0

res = minimize(lambda v: -f(v), x0=[1.0, 0.2], method="SLSQP",
               constraints=[circle], options={"ftol": 1e-12})
x = res.x
print("SLSQP x* =", x.round(4), " f(x*) =", round(f(x), 4), " nit =", res.nit)
vals, vecs = np.linalg.eigh(A)
print("eigh: values", vals.round(4), " top vector", vecs[:, -1].round(4))

g_f, g_h = 2 * A @ x, 2 * x                      # ∇f আর ∇h
lam = (g_f @ g_h) / (g_h @ g_h)                  # ∇f = λ∇h থেকে λ
cross = g_f[0] * g_h[1] - g_f[1] * g_h[0]        # সমান্তরাল হলে 0
print("grad f =", g_f.round(4), " grad h =", g_h.round(4))
print(f"lambda = {lam:.4f}   cross(grad f, grad h) = {cross:.1e}")
print("SLSQP multiplier for -f:", res.multipliers.round(4))

t = np.linspace(-1.6, 1.6, 201)
X, Y = np.meshgrid(t, t)
F = 3 * X**2 + 2 * X * Y + 3 * Y**2
fig, ax = plt.subplots(figsize=(5, 5))
cs = ax.contour(X, Y, F, levels=[1, 2, 4, 6, 9], cmap="viridis")
ax.clabel(cs, fontsize=8)
th = np.linspace(0, 2 * np.pi, 300)
ax.plot(np.cos(th), np.sin(th), "k--", label="constraint: x^2 + y^2 = 1")
ax.plot(*x, "ro", label=f"max on circle, f = {f(x):.2f}")
ax.quiver(*x, *(g_f / 8), color="red", angles="xy", scale_units="xy", scale=1)
ax.set(xlabel="x", ylabel="y", aspect="equal",
       title="Level curve f = 4 touches the circle")
ax.legend(loc="lower left", fontsize=8)

SLSQP x* = [0.7071 0.7071]  f(x*) = 4.0  nit = 9
eigh: values [2. 4.]  top vector [0.7071 0.7071]
grad f = [5.6569 5.6569]  grad h = [1.4142 1.4142]
lambda = 4.0000   cross(grad f, grad h) = -5.7e-08
SLSQP multiplier for -f: [-4.]

৪ · অসমতার শর্ত আর KKT-এর চারটি শর্ত

বাস্তবের বেশিরভাগ শর্ত সমতা নয়, অসমতা: “বাজেটের বেশি নয়”, “ওজন ঋণাত্মক নয়”। অসমতার শর্তের একটি বিশেষ আচরণ আছে: সমাধানে সে হয় দেয়ালে লেগে আছে (g_i(x*) = 0), নয়তো দূরে, ঢিলা (g_i(x*) < 0). প্রথম ক্ষেত্রে শর্তটিকে বলে active, দ্বিতীয় ক্ষেত্রে inactive. Inactive শর্ত সমাধানে কোনো ভূমিকা রাখে না: তাকে মুছে দিলেও উত্তর একই থাকে। Active শর্ত সমতার মতো আচরণ করে, তবে একটি পার্থক্যসহ: দেয়াল শুধু এক দিকে ঠেলতে পারে।

সাধারণ সমস্যার জন্য Lagrangian L(x, μ, λ) = f(x) + ∑μ_ig_i(x) + ∑λ_jh_j(x). এখানে μ_i অসমতার multiplier, λ_j সমতার। x* একটি local সর্বনিম্ন হলে (আর একটি constraint qualification মানলে) এমন μ আর λ থাকে, যাতে নিচের চারটি শর্ত একসঙ্গে মেলে। এদের বলে KKT conditions, Karush, Kuhn আর Tucker-এর নামে:

চতুর্থ শর্তটি সবচেয়ে কাজের: প্রতিটি অসমতার জন্য দুটির একটি শূন্য, হয় μ_i, নয় g_i. তাই ছোট সমস্যা হাতে সমাধানের কৌশল হলো active set অনুমান করা: কোন শর্তগুলো দেয়ালে আছে ধরে নিন, সেগুলোকে সমতা বানিয়ে Lagrange-এর মতো সমাধান করুন, তারপর বাকি তিনটি শর্ত (feasibility আর μ ≥ ০) যাচাই করুন। না মিললে অন্য অনুমান। একটি ঋণাত্মক μ বলে “এই দেয়াল আসলে বিন্দুটিকে টেনে রাখছে”, অর্থাৎ শর্তটি ছেড়ে দিলে f আরও কমত না; সেটিকে active তালিকা থেকে বাদ দিন। নিচের stepper একটি দুই-চলকের সমস্যা পুরো হিসাব করে।

লক্ষ করুন, μ₂ = ২ > μ₁ = ১: এই বিন্দুতে x ≤ ১.৫ সীমাটি বেশি “দামি”, ঢিলা করলে বেশি লাভ। এই তুলনাই অংশ ৮-এর shadow price-এর মূল কথা। আর KKT-এর একটি সীমাবদ্ধতা মনে রাখুন: সাধারণ (non-convex) সমস্যায় KKT শুধু প্রয়োজনীয় শর্ত; একটি KKT বিন্দু local সর্বনিম্ন, local সর্বোচ্চ বা saddle-ও হতে পারে। Convex সমস্যায় তা যথেষ্টও, তাই সেখানে KKT মিলল মানেই global সমাধান।

৫ · KKT সংখ্যায় যাচাই — একটি solver-এর উত্তর বিশ্বাস না করে পরীক্ষা

বড় সমস্যা হাতে নয়, solver দিয়ে সমাধান হয়। কিন্তু solver-এর “success: True” মানে শুধু তার নিজের থামার নিয়ম মিলেছে। উত্তরটি সত্যিই সমাধান কি না, তা KKT দিয়ে আলাদা করে যাচাই করা যায়, আর এতে active শর্তগুলোর তালিকাও পাওয়া যায়। নিচের ব্লকে তিন চলকের একটি quadratic বাটি f(x) = ½xᵀQx − cᵀx (Q positive definite, তাই convex), সঙ্গে পাঁচটি রৈখিক অসমতা: মোট x₁ + x₂ + x₃ ≤ ১.৫, x₁ ≤ ০.২, আর তিনটি x ≥ ০. সব শর্ত একটি matrix আকারে Gx − h ≤ 0, তাই ∇g_i হলো G-এর সারি i, আর stationarity হয় ∇f + Gᵀμ = 0.

SciPy 1.16 থেকে SLSQP ফলের সঙ্গে multipliers ফেরত দেয়। ব্লকটি সেগুলো পড়ে চারটি শর্ত যাচাই করে, তারপর শুধু active সারিগুলো নিয়ে least squares-এ μ আবার নিজে বের করে, যাতে solver-এর ওপর নির্ভর না করেও μ মেলানো যায়। SciPy-তে "ineq" মানে function ≥ 0, তাই h − Gx দেওয়া হয়েছে।

import numpy as np
from scipy.optimize import minimize

np.set_printoptions(precision=4, suppress=True)
Q = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
c = np.array([3.0, 4.0, 3.0])
f = lambda x: 0.5 * x @ Q @ x - c @ x
grad = lambda x: Q @ x - c
print("unconstrained min:", np.linalg.solve(Q, c))

# প্রতিটি constraint g_i(x) = G[i] @ x - h[i] <= 0
G = np.array([[1, 1, 1], [1, 0, 0], [-1, 0, 0], [0, -1, 0], [0, 0, -1.0]])
h = np.array([1.5, 0.2, 0.0, 0.0, 0.0])
names = ["x1+x2+x3 <= 1.5", "x1 <= 0.2", "x1 >= 0", "x2 >= 0", "x3 >= 0"]
cons = [{"type": "ineq", "fun": lambda x: h - G @ x, "jac": lambda x: -G}]
res = minimize(f, np.zeros(3), jac=grad, method="SLSQP",
               constraints=cons, options={"ftol": 1e-12})
x, mu = res.x, res.multipliers
g = G @ x - h
print("x* =", x, " f* =", round(res.fun, 4))
print(f"{'constraint':<17}{'g_i(x*)':>9}{'mu_i':>8}{'mu_i*g_i':>10}  active")
for nm, gi, mi in zip(names, g, mu):
    print(f"{nm:<17}{gi:9.1e}{mi:8.4f}{mi * gi:10.1e}  {abs(gi) < 1e-8}")

print("stationarity |grad f + G.T mu| =", f"{np.linalg.norm(grad(x) + G.T @ mu):.1e}")
print("primal feasible:", bool(np.all(g <= 1e-9)), " dual feasible:", bool(np.all(mu >= -1e-9)))
act = np.abs(g) < 1e-8                 # নিজে আলাদা করে: active সারিতে lstsq
mu_ls = np.linalg.lstsq(G[act].T, -grad(x), rcond=None)[0]
print("mu from lstsq on active rows:", mu_ls)

unconstrained min: [0.5556 0.7778 1.1111]
x* = [0.2 0.7 0.6]  f* = -3.465
constraint         g_i(x*)    mu_i  mu_i*g_i  active
x1+x2+x3 <= 1.5   -1.1e-15  1.1000  -1.2e-15  True
x1 <= 0.2         -1.1e-16  0.4000  -4.4e-17  True
x1 >= 0           -2.0e-01  0.0000  -0.0e+00  False
x2 >= 0           -7.0e-01  0.0000  -0.0e+00  False
x3 >= 0           -6.0e-01  0.0000  -0.0e+00  False
stationarity |grad f + G.T mu| = 1.7e-15
primal feasible: True  dual feasible: True
mu from lstsq on active rows: [1.1 0.4]

৬ · Lagrangian duality — নিচ থেকে সীমা, weak আর strong

Lagrangian-কে আরেকভাবেও দেখা যায়। ধরুন শর্ত মানা বাধ্যতামূলক নয়, কিন্তু ভাঙলে জরিমানা দিতে হয়, প্রতি একক ভাঙায় λ টাকা। প্রতিটি নির্দিষ্ট জরিমানা-হারের জন্য আপনি শর্তহীনভাবে L(x, λ) সর্বনিম্ন করেন। ফলাফল একটি সংখ্যা, যা λ-এর ওপর নির্ভর করে; একে বলে dual function:

q(λ, μ) = min_x L(x, λ, μ),   μ ≥ 0

অংশ ২-এর উদাহরণে L = x² + 2y² − λ(x + y − 3). x আর y-তে ঢাল শূন্য বসালে x = λ/2, y = λ/4, আর বসিয়ে q(λ) = 3λ − 3λ²/8. কয়েকটি মান:

দুটি ঘটনা চোখে পড়ে। প্রথমত, প্রতিটি q(λ) ≤ f* = ৬. এটি কাকতালীয় নয়: যেকোনো feasible x-এ শর্তের পদ শূন্য (বা অসমতায় μg ≤ ০), তাই L(x, λ, μ) ≤ f(x); সব x-এর ওপর সর্বনিম্ন নিলে q আরও ছোট। একেই বলে weak duality, আর এটি সব সমস্যায় খাটে, convex হোক বা না হোক। দ্বিতীয়ত, q-এর সর্বোচ্চ মান λ = ৪-এ, ঠিক ৬, অর্থাৎ ফাঁক শূন্য। একে বলে strong duality, আর ফাঁকটির নাম duality gap. q-কে সর্বোচ্চ করার সমস্যাটিই dual problem, আর মূল সমস্যাটি primal.

Dual কেন দরকার? তিনটি কারণ। (১) Dual-এর চলক শর্তের সংখ্যার সমান, primal-এর চলক feature-এর সংখ্যার সমান। SVM-এ dual-এর চলক প্রতিটি training উদাহরণের জন্য একটি α_i, আর উদাহরণগুলো সেখানে শুধু dot product x_iᵀx_j হিসেবে আসে। ওই dot product-কে kernel দিয়ে বদলানোই kernel trick. (২) Complementary slackness বলে যে উদাহরণ margin থেকে দূরে, তার α_i = ০; শুধু সীমানার উদাহরণগুলো, অর্থাৎ support vector, সমাধান ঠিক করে। (৩) Dual-এর মান সবসময় একটি নিম্নসীমা, তাই primal আর dual-এর মানের পার্থক্য দিয়ে একটি solver জানতে পারে সে সমাধান থেকে কত দূরে; অনেক solver এই ফাঁক দিয়েই থামার নিয়ম বানায়। এই তিনটির পুরো ব্যবহার পর্ব ২৭ (Support Vector Machine — margin, support vector আর hard margin) আর পর্ব ২৮ (Soft margin — hinge loss, C আর dual problem)-এ আসবে।

৭ · Projected gradient descent আর box constraint — ধাপ নিন, তারপর ফিরে আসুন

Lagrange আর KKT শর্ত লেখার ভাষা; কিন্তু হাজার হাজার চলকের সমস্যায় সমাধানের সবচেয়ে সরল algorithm প্রায়ই এটি: সাধারণ gradient descent-এর একটি ধাপ নিন, আর ধাপের পর বিন্দুটি feasible set-এর বাইরে চলে গেলে সবচেয়ে কাছের বৈধ বিন্দুতে ফিরিয়ে আনুন। এর নাম projected gradient descent:

x_k+1 = P_C( x_k − η ∇f(x_k) )

এখানে C হলো feasible set, η ধাপের মাপ (learning rate), আর P_C(v) মানে C-এর যে বিন্দু v-এর সবচেয়ে কাছে, অর্থাৎ ‖w − v‖² সর্বনিম্ন করে এমন w ∈ C. C convex হলে এই বিন্দু একটিই। ধারণাটি পুরোনো: ১৯৬০-এর দশকের মাঝামাঝি Goldstein আর Levitin ও Polyak-এর কাজে এর বিশ্লেষণ আছে বলে সাধারণত উল্লেখ করা হয়। কৌশলটির শক্তি নির্ভর করে projection কত সস্তা তার ওপর:

• Box l ≤ x ≤ u: projection মানে প্রতিটি ঘর আলাদা করে কেটে দেওয়া, np.clip(v, l, u). খরচ n.
• Ball ‖x‖ ≤ r: v বাইরে থাকলে r·v/‖v‖-এ ছোট করা। Gradient clipping-এর সঙ্গে মিল লক্ষ করুন।
• Probability simplex {w ≥ 0, ∑w = 1}: সব ঘর থেকে একই সংখ্যা θ বাদ দিয়ে শূন্যে কাটা; θ খুঁজতে sort লাগে, খরচ n log n. নিচের stepper হাতে একবার করে।
• সাধারণ রৈখিক শর্ত Ax ≤ b: projection নিজেই একটি QP; তখন SLSQP বা interior-point solver ভালো।

Box শর্তের জন্য SciPy-তে আলাদা একটি পরিচিত টুল আছে: L-BFGS-B, যা পর্ব ১৫ (Second-order পদ্ধতি — Newton, L-BFGS আর line search)-এর L-BFGS-কে bounds মানতে শেখায় (Byrd, Lu, Nocedal আর Zhu, ১৯৯৫)। এটি শুধু box নেয়, সাধারণ শর্ত নয়। নিচের ব্লক অংশ ৫-এর একই বাটিতে দুটি সমস্যা চালায়: box [০, ০.৮]³-এ projected GD বনাম L-BFGS-B, আর simplex-এ projected GD বনাম SLSQP. ধাপের মাপ η = 1/L, যেখানে L হলো Q-এর সবচেয়ে বড় eigenvalue (পর্ব ১৩-এ এর কারণ আসবে)।

import numpy as np
from scipy.optimize import minimize

np.set_printoptions(precision=4, suppress=True)
Q = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
c = np.array([3.0, 4.0, 3.0])
f = lambda x: 0.5 * x @ Q @ x - c @ x
grad = lambda x: Q @ x - c
eta = 1 / np.linalg.eigvalsh(Q).max()            # ধাপ 1/L

def proj_box(v):                                  # বাক্স [0, 0.8]³: শুধু কেটে দেওয়া
    return np.clip(v, 0.0, 0.8)
def proj_simplex(v):                              # সব ঘর ≥ 0, যোগফল 1
    u = np.sort(v)[::-1]
    css = np.cumsum(u)
    k = np.arange(1, v.size + 1)
    rho = np.nonzero(u - (css - 1) / k > 0)[0][-1]
    theta = (css[rho] - 1) / (rho + 1)
    return np.maximum(v - theta, 0.0)
def pgd(proj, x, tol=1e-10, max_it=500):
    for it in range(1, max_it + 1):
        x_new = proj(x - eta * grad(x))           # gradient ধাপ, তারপর projection
        if np.linalg.norm(x_new - x) < tol:
            return x_new, it
        x = x_new
    return x, max_it
xb, nb = pgd(proj_box, np.zeros(3))
lb = minimize(f, np.zeros(3), jac=grad, method="L-BFGS-B", bounds=[(0, 0.8)] * 3)
print("box     PGD:", xb, f"iters {nb}", "| L-BFGS-B:", lb.x, f"iters {lb.nit}")
xs, ns = pgd(proj_simplex, np.ones(3) / 3)
sl = minimize(f, np.ones(3) / 3, jac=grad, method="SLSQP", bounds=[(0, None)] * 3,
              constraints=[{"type": "eq", "fun": lambda x: x.sum() - 1}],
              options={"ftol": 1e-12})
print("simplex PGD:", xs, f"iters {ns}", "| SLSQP:", sl.x.round(4), f"iters {sl.nit}")
print("projection of (0.5, 1.2, -0.3):", proj_simplex(np.array([0.5, 1.2, -0.3])))

box     PGD: [0.55 0.8  0.8 ] iters 13 | L-BFGS-B: [0.55 0.8  0.8 ] iters 3
simplex PGD: [0.1429 0.5714 0.2857] iters 57 | SLSQP: [0.1429 0.5714 0.2857] iters 7
projection of (0.5, 1.2, -0.3): [0.15 0.85 0.  ]

৮ · বাজারঘরের বিজ্ঞাপন বাজেট — SLSQP, linprog আর shadow price

বাজারঘর আগামী মাসে বিজ্ঞাপনে মোট ১ লাখ টাকা (১০০ হাজার) খরচ করবে, চারটি চ্যানেলে: সোশ্যাল মিডিয়া, সার্চ বিজ্ঞাপন, ভিডিও আর SMS. প্রতিটি চ্যানেলে খরচের একটি ঊর্ধ্বসীমা আছে (চুক্তি বা দর্শকের সীমা)। এখানের সব সংখ্যা কাল্পনিক, শুধু পদ্ধতি দেখানোর জন্য। প্রতিটি চ্যানেলে ফেরত diminishing returns মানে: প্রথম কয়েক হাজার টাকায় অনেক বিক্রি, পরের টাকায় কম। একে মডেল করা হয়েছে r_i(x) = a_i log(1 + x/b_i) দিয়ে, যার প্রান্তিক ফেরত r_i′(x) = a_i/(b_i + x) খরচ বাড়লে কমে। বাস্তবে এই function-গুলো আসে অতীতের ডেটায় fit করা একটি model থেকে, আর তার নিজস্ব অনিশ্চয়তা আছে।

KKT এখানে একটি সুন্দর ব্যবসায়িক নিয়ম দেয়। বাজেট শর্তের multiplier ν হলে, যে চ্যানেল সীমায় বা শূন্যে আটকে নেই, তার জন্য stationarity বলে r_i′(x_i) = ν: সব “মুক্ত” চ্যানেলে শেষ হাজার টাকার ফেরত সমান। কোনো চ্যানেলে এটি বেশি হলে সেখানে টাকা সরানো লাভজনক হতো। যে চ্যানেল নিজের সীমায় আটকে আছে, তার প্রান্তিক ফেরত ν-এর চেয়ে বেশি, পার্থক্যটি তার সীমার multiplier. দ্বিতীয় অংশে একই প্রশ্নের একটি LP সংস্করণ: প্রতি হাজার টাকায় স্থির ফেরত, যা linprog (HiGHS solver) দিয়ে সমাধান হয়, আর solver নিজেই প্রতিটি শর্তের dual value, অর্থাৎ shadow price দেয়।

import numpy as np
from scipy.optimize import minimize, linprog

np.set_printoptions(precision=2, suppress=True)
# ক্রম: social, search, video, sms — সব সংখ্যা কাল্পনিক, হাজার টাকায়
a = np.array([120.0, 90.0, 150.0, 50.0])         # ফেরতের উচ্চতা
b = np.array([20.0, 10.0, 60.0, 5.0])            # কত দ্রুত ফেরত কমে
cap = np.array([50.0, 30.0, 60.0, 15.0])         # প্রতি চ্যানেলের সীমা
B = 100.0                                        # মোট বাজেট

rev = lambda x: np.sum(a * np.log1p(x / b))      # কমতে-থাকা ফেরত
d_rev = lambda x: a / (b + x)                    # প্রতি হাজার টাকায় প্রান্তিক ফেরত
budget = {"type": "ineq", "fun": lambda x: B - x.sum(), "jac": lambda x: -np.ones(4)}
res = minimize(lambda x: -rev(x), np.full(4, 25.0), jac=lambda x: -d_rev(x),
               method="SLSQP", bounds=list(zip(np.zeros(4), cap)),
               constraints=[budget], options={"ftol": 1e-12})
x = res.x
print("SLSQP spend :", x, " total", round(x.sum(), 2))
print("marginal    :", d_rev(x), " revenue", round(rev(x), 2))
print("budget multiplier:", res.multipliers.round(4))

r = np.array([2.4, 3.1, 1.6, 2.8])              # LP: প্রতি হাজার টাকায় স্থির ফেরত
lp = linprog(-r, A_ub=np.ones((1, 4)), b_ub=[B],
             bounds=list(zip(np.zeros(4), cap)), method="highs")
print("LP spend    :", lp.x, " revenue", round(-lp.fun, 2))
print("shadow price of budget  :", 0.0 - lp.ineqlin.marginals)
print("shadow price of each cap:", 0.0 - lp.upper.marginals)
lp2 = linprog(-r, A_ub=np.ones((1, 4)), b_ub=[B + 1],
              bounds=list(zip(np.zeros(4), cap)), method="highs")
print("revenue with budget 101 :", round(-lp2.fun, 2))

SLSQP spend : [40. 30. 15. 15.]  total 100.0
marginal    : [2.   2.25 2.   2.5 ]  revenue 359.39
budget multiplier: [2.]
LP spend    : [50. 30.  5. 15.]  revenue 263.0
shadow price of budget  : [1.6]
shadow price of each cap: [0.8 1.5 0.  1.2]
revenue with budget 101 : 264.6

৯ · ML-এ constraint কোথায় দাঁড়ায় — শাস্তি না শর্ত, আর কোন টুল

ML-এর বেশিরভাগ training loop-এ সরাসরি শর্ত দেখা যায় না, কারণ শর্তকে প্রায়ই শাস্তিতে রূপান্তর করা হয়। Ridge-এর দুটি রূপ ভাবুন: “‖y − Xw‖² সর্বনিম্ন করুন, শর্ত ‖w‖² ≤ t” আর “‖y − Xw‖² + α‖w‖² সর্বনিম্ন করুন”। প্রথমটির KKT লিখলে stationarity হয় ∇‖y − Xw‖² + μ · 2w = 0, যা দ্বিতীয়টির gradient-এর সঙ্গে হুবহু এক যখন α = μ. অর্থাৎ convex সমস্যায় প্রতিটি t-এর জন্য এমন একটি α আছে যার সমাধান একই, আর α হলো শর্তের Lagrange multiplier. তবে t আর α-এর সম্পর্ক ডেটার ওপর নির্ভর করে, আগে থেকে জানা যায় না। AI-এর ভিত্তি কোর্সের পর্ব ৬১-এ L1 আর L2-এর জ্যামিতির ছবিতে যে হীরা আর বৃত্ত দেখেছেন, সেগুলোই এই constrained রূপের feasible set; পুরো তত্ত্ব পর্ব ১৭ (Regularization-এর তত্ত্ব — L1/L2-এর জ্যামিতি, MAP আর early stopping)-এ।

কখন শাস্তি যথেষ্ট নয়, সত্যিকারের শর্ত লাগে? যখন শর্ত ভাঙলে ফলের কোনো মানে থাকে না: সম্ভাবনা ঋণাত্মক হতে পারে না, মজুদ শূন্যের নিচে যেতে পারে না, বাজেট পেরোনো যায় না। Fairness-এর কিছু পদ্ধতিও শর্ত হিসেবে লেখা হয়, যদিও কোন মাপটি ঠিক, তা নিয়ে বিতর্ক আছে।

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

AI assistant constrained optimization-এর কোড দ্রুত লিখে দেয়, আর বেশিরভাগ সময় কোডটি চলেও। সমস্যা হয় চুপচাপ ভুলে: শর্তের চিহ্ন উল্টো, ভুল solver, বা solver-এর “success” দেখেই উত্তর বিশ্বাস। এই পর্বের KKT যাচাই ঠিক এই ভুলগুলো ধরার যন্ত্র।

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

• সমস্যাটি নিজে standard form-এ লিখেছি, আর প্রতিটি শর্তের চিহ্ন solver-এর নিয়মের সঙ্গে মিলিয়েছি
• সমাধানে প্রতিটি শর্তের মান ছেপে দেখেছি কোনগুলো active
• Stationarity residual, μ ≥ ০ আর μ_ig_i ≈ ০ সংখ্যায় পরীক্ষা করেছি
• ছোট একটি সংস্করণ হাতে বা আরেকটি solver দিয়ে সমাধান করে মিলিয়েছি
• Non-convex হলে কয়েকটি আলাদা শুরুর বিন্দু থেকে চালিয়েছি
• Shadow price শর্তের মান বদলে আবার solve করে যাচাই করেছি
• ফেরতের model যে ডেটা থেকে এসেছে, তা আলাদা validation ডেটায় পরীক্ষা করেছি; optimizer ভুল model-কেও নিখুঁতভাবে optimize করে
• কোনো আসল গ্রাহক বা ব্যবসার গোপন ডেটা AI-তে পেস্ট করিনি