Positive definite matrix — quadratic form আর covariance

AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ matrix গুণ, solve আর condition number এসেছিল; পর্ব ৩৭-এ দুটি কলামের covariance আর correlation; আর পর্ব ৬৫-এ covariance-এর বড় দিক খুঁজে PCA. এই মডিউলের পর্ব ০৮ (Eigenvalue আর eigenvector — ধারণা থেকে PCA পর্যন্ত) প্রতিসম matrix-এর spectral theorem দেখায়। এই পর্বের প্রশ্ন: একটি প্রতিসম matrix কখন “ভালো আচরণের”, অর্থাৎ বাটির মতো, যার একটিই তলা আছে? উত্তরটি ML-এর প্রায় সব জায়গায় দরকার হয়: covariance, Gaussian, ridge, Newton পদ্ধতি, kernel. চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy আর matplotlib); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে।

১ · Quadratic form xᵀAx — একটি matrix থেকে একটি উচ্চতার মানচিত্র

একটি বর্গাকার matrix A (আকার n × n) আর একটি vector x (দৈর্ঘ্য n) দিলে xᵀAx একটি সংখ্যা: প্রথমে Ax (একটি vector), তারপর x-এর সঙ্গে তার dot product. ᵀ মানে transpose, এখানে x-কে কলাম থেকে সারি বানানো। এই function-কে বলে quadratic form, কারণ এর প্রতিটি পদ x-এর দুটি উপাদানের গুণফল, অর্থাৎ ঘাত ২. ২ × ২ ক্ষেত্রে A = [[a, b], [c, d]] হলে:

xᵀAx = a·x₁² + (b + c)·x₁x₂ + d·x₂²

এখানে x₁, x₂ হলো x-এর দুটি উপাদান, a আর d কর্ণের দুটি ঘর, b আর c কর্ণের বাইরের দুটি। লক্ষ করুন, b আর c আলাদা করে আসে না, আসে শুধু তাদের যোগফল। তাই [[3, 1], [1, 2]] আর [[3, 2], [0, 2]] একই quadratic form দেয়। সাধারণ প্রমাণ এক লাইনের: xᵀAx একটি সংখ্যা, আর একটি সংখ্যার transpose সে নিজেই, তাই xᵀAx = (xᵀAx)ᵀ = xᵀAᵀx. দুটির গড় নিলে xᵀAx = xᵀ S x, যেখানে S = (A + Aᵀ)/2 হলো A-এর প্রতিসম অংশ। এজন্য এই পর্বে, আর সাধারণভাবে definiteness-এর আলোচনায়, A সবসময় symmetric matrix (A = Aᵀ) ধরা হয়।

শেষ ধাপের ঘটনাটিই এই পর্বের মূল যন্ত্র। পর্ব ০৮-এর spectral theorem বলে, প্রতিসম A-কে লেখা যায় A = VΛVᵀ, যেখানে V-এর কলামগুলো পরস্পর লম্ব একক eigenvector, আর Λ কর্ণে eigenvalue-গুলো। তাহলে y = Vᵀx বসালে (V শুধু ঘোরায়, দৈর্ঘ্য বদলায় না):

xᵀAx = xᵀVΛVᵀx = yᵀΛy = λ₁y₁² + λ₂y₂² + … + λ_ny_n²

প্রতিটি y_i² ≥ ০, তাই xᵀAx-এর চিহ্ন পুরোপুরি λ-গুলোর চিহ্নের ওপর নির্ভর করে। আর একক দৈর্ঘ্যের x-এ xᵀAx-এর সর্বনিম্ন মান ঠিক ক্ষুদ্রতম eigenvalue, সর্বোচ্চ মান বৃহত্তম (এই অনুপাত xᵀAx / xᵀx-কে বলে Rayleigh quotient).

২ · তিন আকার — বাটি, জিন আর উপত্যকা

উপরের সূত্র থেকে সরাসরি পাঁচটি সংজ্ঞা আসে। প্রতিসম A-এর জন্য, শূন্য নয় এমন সব x-এ:

PD হলে PSD-ও, কিন্তু উল্টোটা নয়। চিহ্নগুলো লেখা হয় A ≻ 0 (PD) আর A ⪰ 0 (PSD). নিচের ছবিতে তিনটি ২ × ২ matrix-এর contour, অর্থাৎ সমান উচ্চতার রেখা, যেমন মানচিত্রে পাহাড়ের উচ্চতা-রেখা।

উপবৃত্তের জ্যামিতি সরাসরি eigenvalue থেকে আসে। xᵀAx = c রেখায়, eigen-অক্ষে λ₁y₁² + λ₂y₂² = c, তাই অর্ধ-অক্ষের দৈর্ঘ্য √(c/λ₁) আর √(c/λ₂). ছোট λ মানে লম্বা অক্ষ। দুটি λ খুব আলাদা হলে উপবৃত্ত সরু আর লম্বা, যা gradient descent-এর জন্য কঠিন উপত্যকা; এর মাপ হলো condition number λ_max/λ_min, এখানে ৩.৬১৮/১.৩৮২ = ২.৬২. পর্ব ১৩ (Gradient descent-এর পরিবার) এই সংখ্যাটি দিয়েই অভিসরণের গতি ব্যাখ্যা করবে।

নিচের ব্লক তিনটি matrix-এর contour আঁকে, eigenvalue ছাপে, আর প্রতিটি প্যানেলে ক্ষুদ্রতম eigenvalue-এর দিকে একটি কালো তীর দেয়।

import numpy as np
import matplotlib.pyplot as plt

mats = {"bowl: positive definite": np.array([[3.0, 1.0], [1.0, 2.0]]),
        "saddle: indefinite":      np.array([[1.0, 2.0], [2.0, 1.0]]),
        "valley: semidefinite":    np.array([[1.0, 1.0], [1.0, 1.0]])}

g = np.linspace(-2, 2, 201)
X1, X2 = np.meshgrid(g, g)
P = np.stack([X1, X2], axis=-1)              # প্রতিটি গ্রিড-বিন্দু x = (x1, x2)
levels = [-4, -2, -1, -0.5, 0.5, 1, 2, 4, 8]

fig, axes = plt.subplots(1, 3, figsize=(9, 3.3), layout="constrained")
for ax, (name, A) in zip(axes, mats.items()):
    Q = np.einsum("...i,ij,...j->...", P, A, P)   # প্রতিটি বিন্দুতে xᵀAx
    lam, V = np.linalg.eigh(A)                    # ছোট থেকে বড় eigenvalue
    print(f"{name:24s} λ = {np.round(lam, 3)}  min xᵀAx = {Q.min():6.2f}")
    print(f"{'':24s} ছোট λ-এর দিক v = {np.round(V[:, 0], 3)}")
    ax.contour(X1, X2, Q, levels=levels, cmap="RdBu_r")
    ax.quiver(0, 0, *V[:, 0], color="k", angles="xy", scale_units="xy", scale=1)
    ax.set_title(name, fontsize=10)
    ax.set_xlabel("x1"); ax.set_aspect("equal")
axes[0].set_ylabel("x2")
plt.show()

bowl: positive definite  λ = [1.382 3.618]  min xᵀAx =   0.00
                         ছোট λ-এর দিক v = [ 0.526 -0.851]
saddle: indefinite       λ = [-1.  3.]  min xᵀAx =  -8.00
                         ছোট λ-এর দিক v = [-0.707  0.707]
valley: semidefinite     λ = [0. 2.]  min xᵀAx =  -0.00
                         ছোট λ-এর দিক v = [-0.707  0.707]

৩ · PD-র তিন পরীক্ষা — eigenvalue, leading minor আর Cholesky

সংজ্ঞা “সব x-এ xᵀAx > ০” সরাসরি যাচাই করা যায় না, কারণ x অসীম। তাই তিনটি সমতুল্য পরীক্ষা ব্যবহার হয়, প্রতিসম A-এর জন্য:

২ × ২-এ Sylvester-এর নিয়ম হয় খুব সরল: [[a, b], [b, d]] PD হবে যদি এবং কেবল যদি a > ০ আর ad − b² > ০. কেন? বর্গ পূরণ করলে a·x₁² + 2b·x₁x₂ + d·x₂² = a(x₁ + (b/a)x₂)² + ((ad − b²)/a)·x₂². দুটি সহগ a আর (ad − b²)/a; দুটিই ধনাত্মক হলে যোগফল ধনাত্মক। অংশ ১-এর stepper-এ এই দুটিই ছিল ৩ আর ৫/৩. এই সহগগুলোকে বলে pivot; k-তম pivot = (k-তম minor)/(আগের minor), তাই সব minor ধনাত্মক ⇔ সব pivot ধনাত্মক ⇔ Cholesky সফল (কারণ Cholesky-র কর্ণ ঠিক pivot-গুলোর বর্গমূল)। তিনটি পরীক্ষা একই সত্যের তিন রূপ।

এবার তিন পরীক্ষা একসঙ্গে, সাতটি matrix-এ। শেষ দুটি ইচ্ছাকৃত ফাঁদ: F-এর leading minor দুটিই ০ (≥ ০), তবু সে PSD নয়; আর G প্রতিসম নয়।

import numpy as np
from scipy import linalg

def eig_test(A):                      # সব eigenvalue > 0?
    return f"{np.linalg.eigvalsh(A).min():+.3f}"

def minors_test(A):                   # Sylvester: উপরের-বাঁ কোণের k×k det
    return " ".join(f"{np.linalg.det(A[:k, :k]):+.0f}" for k in range(1, len(A) + 1))

def chol_test(A):                     # Cholesky সফল হলো, নাকি LinAlgError?
    try:
        linalg.cholesky(A, lower=True)
        return "ok"
    except linalg.LinAlgError as err:
        return f"LinAlgError: {err}"

cases = {
    "A bowl":   [[3, 1], [1, 2]],
    "B saddle": [[1, 2], [2, 1]],
    "C valley": [[1, 1], [1, 1]],
    "D 3x3":    [[4, 2, 0], [2, 3, 1], [0, 1, 2]],
    "E 3x3":    [[2, 3, 0], [3, 2, 1], [0, 1, 2]],
    "F trap":   [[0, 0], [0, -1]],
    "G nonsym": [[1, 4], [0, 1]],
}
print(f"{'matrix':9} {'min λ':>7} {'leading minors':>15}  Cholesky")
for name, rows in cases.items():
    A = np.array(rows, dtype=float)
    print(f"{name:9} {eig_test(A):>7} {minors_test(A):>15}  {chol_test(A)}")

G = np.array(cases["G nonsym"], dtype=float)
x = np.array([1.0, -1.0])
S = (G + G.T) / 2                     # quadratic form দেখে শুধু এই প্রতিসম অংশ
print("G: xᵀGx =", x @ G @ x, "| xᵀSx =", x @ S @ x)
print("S =", S.tolist(), "| eig(S) =", np.linalg.eigvalsh(S))

matrix      min λ  leading minors  Cholesky
A bowl     +1.382           +3 +5  ok
B saddle   -1.000           +1 -3  LinAlgError: Internal potrf return info = [2] for slices [0].
C valley   +0.000           +1 +0  LinAlgError: Internal potrf return info = [2] for slices [0].
D 3x3      +0.855       +4 +8 +12  ok
E 3x3      -1.162       +2 -5 -12  LinAlgError: Internal potrf return info = [2] for slices [0].
F trap     -1.000           +0 +0  LinAlgError: Internal potrf return info = [1] for slices [0].
G nonsym   +1.000           +1 +1  ok
G: xᵀGx = -2.0 | xᵀSx = -2.0
S = [[1.0, 2.0], [2.0, 1.0]] | eig(S) = [-1.  3.]

৪ · Cholesky decomposition — হাতে একটি ২ × ২, তারপর কেন এত জনপ্রিয়

PD matrix A-কে লেখা যায় A = LLᵀ, যেখানে L একটি lower triangular matrix (কর্ণের ওপরে সব শূন্য) আর কর্ণের সব ঘর ধনাত্মক। এই L একটিই (unique). এটি যেন matrix-এর “বর্গমূল”: সংখ্যায় ৯ = ৩ · ৩, matrix-এ A = L · Lᵀ. পদ্ধতিটির নাম ফরাসি সামরিক প্রকৌশলী André-Louis Cholesky-র নামে, যিনি জরিপের হিসাবের জন্য এটি বানান; ভাবনার প্রশ্ন ১-এ তাঁর গল্প আছে।

সাধারণ n × n-এ একই নিয়ম কলাম ধরে চলে: l_jj = √(a_jj − ∑_k<j l_jk²), আর j-এর নিচের সারিতে l_ij = (a_ij − ∑_k<j l_ikl_jk) / l_jj. খরচ প্রায় n³/3টি গুণ-যোগ, যেখানে সাধারণ LU decomposition-এর প্রায় 2n³/3; প্রতিসমতার কারণে অর্ধেক কাজ বাঁচে। আর PD matrix-এ সারি অদলবদল (pivoting) ছাড়াই পদ্ধতিটি সাংখ্যিকভাবে স্থিতিশীল। তত্ত্ব বলে প্রায় দ্বিগুণ দ্রুত; বাস্তবে কতটা? একটি ৩,০০০ × ৩,০০০ PD system-এ চারটি পথ তুলনা করা হলো (ডেস্কটপে সত্যিকারের run; Python 3.14.4, NumPy 2.5.4, SciPy 1.18.1; প্রতিটির ৫ বারের সেরা সময়, একটি local run, সাধারণ ফল নয়):

~/bazarghar-ml$ python chol_speed.py
numpy 2.5.4 | scipy 1.18.1 | n = 3000
LU   solve(assume_a='gen')        202.4 ms   residual 1.4e-13
Cholesky cho_factor+cho_solve     171.6 ms   residual 9.4e-14
solve(assume_a='pos')             171.2 ms   residual 9.4e-14
inv(A) @ b                        390.5 ms   residual 1.3e-13

Matrix-টি বানানো হয়েছিল XᵀX/m + ০.১·I হিসেবে, ঠিক ridge regression-এর মতো (অংশ ৮)। Cholesky এখানে LU-এর চেয়ে প্রায় ১৫% দ্রুত, দ্বিগুণ নয়: আধুনিক BLAS একাধিক core আর memory-র গতিতে আটকে যায়, তাই গণনার সংখ্যার অনুপাত সময়ে হুবহু আসে না। Inverse বানিয়ে গুণ করা সবচেয়ে ধীর, প্রায় দ্বিগুণ সময়; AI-এর ভিত্তি কোর্সের পর্ব ১৭-এর “inverse নয়, solve” নিয়ম এখানেও খাটে। SciPy-তে solve(A, b, assume_a="pos") ভেতরে Cholesky-ই চালায়, তাই সময় প্রায় এক। আপনার যন্ত্রে সংখ্যা আলাদা হবে।

৫ · Covariance matrix সবসময় PSD — প্রমাণ, rank কমে যাওয়া আর jitter

d-টি feature-এর covariance matrix Σ একটি d × d প্রতিসম matrix: কর্ণে প্রতিটি feature-এর variance, বাইরে প্রতিটি জোড়ার covariance. এর সবচেয়ে গুরুত্বপূর্ণ গুণ এক লাইনে প্রমাণ হয়। যেকোনো ধ্রুবক vector v নিন, আর নতুন একটি সংখ্যা বানান s = vᵀx, অর্থাৎ feature-গুলোর একটি ওজন-দেওয়া যোগফল। তাহলে

Var(vᵀx) = E[(vᵀ(x − μ))²] = vᵀ E[(x − μ)(x − μ)ᵀ] v = vᵀΣv ≥ 0

কারণ variance একটি বর্গের গড়, কখনো ঋণাত্মক নয়। এখানে μ হলো গড় vector, E প্রত্যাশিত মান। হাতের ডেটায় একই যুক্তি: কেন্দ্রীকৃত design matrix X_c (প্রতিটি কলাম থেকে তার গড় বাদ) নিয়ে sample covariance S = X_cᵀX_c/(n − 1), আর vᵀSv = ‖X_cv‖²/(n − 1) ≥ ০. দৈর্ঘ্যের বর্গ ঋণাত্মক হয় না। সমান চিহ্ন কখন? যখন X_cv = ০, অর্থাৎ feature-গুলোর কোনো একটি মিশ্রণ সব সারিতে ধ্রুবক। তখন S PSD কিন্তু PD নয়, আর তার rank d-এর চেয়ে কম। এটি ঘটে দুই ক্ষেত্রে: (১) একটি কলাম অন্যগুলোর রৈখিক ফল (পর্ব ০৬-এর linear dependence), (২) n − 1 < d, অর্থাৎ feature-এর চেয়ে সারি কম, কারণ S-এর rank সর্বোচ্চ n − 1.

নিচের ব্লক কাল্পনিক বাজারঘরের ২০০টি অর্ডারে (দূরত্ব, ডেলিভারির সময়, পণ্যের সংখ্যা) প্রমাণটি সংখ্যায় দেখায়, তারপর দুটি নির্ভরশীল কলাম যোগ করে: ডেলিভারি চার্জ, যা দূরত্বের হুবহু রৈখিক ফল (২০ + ৫ × km), আর একই ডেলিভারির সময় ঘণ্টায়, যা ডেটা জোড়া লাগানোর সময় প্রায়ই ঘটে।

import numpy as np
from scipy import linalg
np.set_printoptions(precision=3, linewidth=70)

def chol(S):                          # Cholesky: সফল হলে L-এর ক্ষুদ্রতম কর্ণ
    try:
        return f"ok, min L diag {linalg.cholesky(S, lower=True).diagonal().min():.1e}"
    except linalg.LinAlgError as err:
        return f"LinAlgError: {err}"

rng = np.random.default_rng(42)       # কাল্পনিক বাজারঘর: ২০০টি অর্ডার
n = 200
dist = rng.uniform(1, 15, n)                          # দূরত্ব, km
time = 15 + 2.5 * dist + rng.normal(0, 6, n)          # ডেলিভারির সময়, মিনিট
items = rng.integers(1, 8, n).astype(float)           # পণ্যের সংখ্যা
X = np.column_stack([dist, time, items])
S = np.cov(X, rowvar=False)           # 3×3 sample covariance
print("eig(S):", np.linalg.eigvalsh(S))
v = rng.normal(size=3)                # যেকোনো একটি দিক v
print(f"vᵀSv = {v @ S @ v:.4f} | var(Xv) = {np.var(X @ v, ddof=1):.4f}")

extras = {"charge = 20 + 5*dist": 20 + 5 * dist,      # দূরত্বের হুবহু রৈখিক ফল
          "time in hours":        time / 60}          # একই তথ্য, অন্য এককে
for name, col in extras.items():
    S2 = np.cov(np.column_stack([X, col]), rowvar=False)
    lam = np.linalg.eigvalsh(S2)
    print(f"+ {name}: rank {np.linalg.matrix_rank(S2)}, "
          f"min λ {lam[0]:.1e}, max λ {lam[-1]:.1e}")
    print(f"  Cholesky: {chol(S2)}")

S2 = np.cov(np.column_stack([X, time / 60]), rowvar=False)
for eps in (1e-10, 1e-8, 1e-6):       # jitter: S2 + εI
    print(f"eps={eps:g}: {chol(S2 + eps * np.eye(4))}")

eig(S): [  3.31    3.951 169.278]
vᵀSv = 6.7850 | var(Xv) = 6.7850
+ charge = 20 + 5*dist: rank 3, min λ -2.1e-15, max λ 5.4e+02
  Cholesky: ok, min L diag 3.4e-07
+ time in hours: rank 3, min λ -3.9e-17, max λ 1.7e+02
  Cholesky: LinAlgError: Internal potrf return info = [4] for slices [0].
eps=1e-10: ok, min L diag 1.0e-05
eps=1e-08: ok, min L diag 1.0e-04
eps=1e-06: ok, min L diag 1.0e-03

Jitter কী সারায়, কী লুকায়। কর্ণে ছোট ε যোগ করা Gaussian process, kernel পদ্ধতি আর বড় covariance-এর হিসাবে একটি প্রচলিত কৌশল। এটি rounding-এর ক্ষুদ্র ঋণাত্মক eigenvalue সারায় আর Cholesky চালু রাখে। কিন্তু উপরের ক্ষেত্রে আসল সমস্যা ডেটার: একই তথ্য দুবার আছে। সঠিক সমাধান একটি কলাম ফেলে দেওয়া; jitter শুধু উপসর্গ ঢাকে। ε বাছাইয়েও সাবধান: এটি matrix-এর মাপের সঙ্গে তুলনীয় হতে হয় (যেমন ε = ১০^−৮ × গড় কর্ণ), কারণ কর্ণে ১০০-এর মতো সংখ্যা থাকলে ১০^−১০ rounding-এর নিচে হারিয়ে যেতে পারে, আর বড় ε model-কে সত্যিই বদলে দেয় (ridge-এর λ-এর মতো)। আরেকটি সাধারণ ফাঁদ: হাতে বসানো বা আংশিক-অনুপস্থিত ডেটা থেকে জোড়ায় জোড়ায় হিসাব করা correlation matrix PSD না-ও হতে পারে; ভাবনার প্রশ্ন ৩ দেখুন।

একই যুক্তি আরেকটি বড় শ্রেণিতে খাটে। যেকোনো X-এর জন্য XᵀX (d × d) আর XXᵀ (n × n) দুটিই PSD, কারণ vᵀXᵀXv = ‖Xv‖². দ্বিতীয়টি হলো সারিগুলোর জোড়ায় জোড়ায় dot product-এর টেবিল, যাকে বলে Gram matrix. পর্ব ২৯ (Kernel trick)-এ dot product-এর জায়গায় একটি kernel function k(x_i, x_j) বসবে, আর “বৈধ kernel” মানেই এমন function যার Gram matrix সব ডেটায় PSD.

৬ · Cholesky দিয়ে correlated Gaussian — z থেকে x = μ + Lz

z ~ N(0, I) মানে z-এর উপাদানগুলো স্বাধীন standard normal, গড় ০, variance ১, তাই Cov(z) = I (identity matrix). এখন x = μ + Lz নিন, যেখানে Σ = LLᵀ. তাহলে

E[x] = μ, Cov(x) = E[(Lz)(Lz)ᵀ] = L E[zzᵀ] Lᵀ = L I Lᵀ = LLᵀ = Σ

আর Gaussian-এর রৈখিক রূপান্তর Gaussian-ই থাকে, তাই x ~ N(μ, Σ), একটি multivariate Gaussian. লক্ষ করুন, Σ = LLᵀ-এর যেকোনো বর্গমূল কাজ করত (যেমন eigen থেকে VΛ^½), কিন্তু Cholesky সবচেয়ে সস্তা, আর নিচের-ত্রিভুজ হওয়ায় একটি সুন্দর অর্থও আছে: x₁ = μ₁ + l₁₁z₁ শুধু প্রথম ছক্কা থেকে, x₂ = μ₂ + l₂₁z₁ + l₂₂z₂ প্রথমটির কিছু অংশ ধার নেয়। Σ PD না হলে Cholesky হয় না, আর PSD-ও না হলে এমন কোনো Gaussian-ই নেই: ঋণাত্মক variance-এর দিক বাস্তবে অসম্ভব।

বাজারঘরের কাল্পনিক ডেলিভারি: সময়ের গড় ৩৫ মিনিট, standard deviation ১০; দূরত্বের গড় ৬ km, standard deviation ৩; correlation ρ = ০.৭. তাহলে Σ = [[১০০, ২১], [২১, ৯]] (কারণ ০.৭ × ১০ × ৩ = ২১)। হাতে Cholesky: l₁₁ = ১০, l₂₁ = ২১/১০ = ২.১, l₂₂ = √(৯ − ৪.৪১) = √৪.৫৯ ≈ ২.১৪২.

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

mu = np.array([35.0, 6.0])            # গড়: ৩৫ মিনিট, ৬ km
sd = np.array([10.0, 3.0]); rho = 0.7
Sigma = np.array([[sd[0]**2,            rho * sd[0] * sd[1]],
                  [rho * sd[0] * sd[1], sd[1]**2]])
L = np.linalg.cholesky(Sigma)         # Sigma = L Lᵀ, L নিচের-ত্রিভুজ
print("Sigma =\n", Sigma, "\nL =\n", L)

rng = np.random.default_rng(42)
Z = rng.standard_normal((5000, 2))    # z ~ N(0, I): পরস্পর-স্বাধীন
Xs = mu + Z @ L.T                     # প্রতিটি সারিতে x = μ + L z
print("sample mean:", Xs.mean(axis=0))
print("sample cov:\n", np.cov(Xs, rowvar=False))
print("sample corr:", round(np.corrcoef(Xs, rowvar=False)[0, 1], 3))

orders = np.array([[55.0, 12.0], [55.0, 3.0]])   # দুটি নতুন অর্ডার: (মিনিট, km)
for x in orders:
    y = np.linalg.solve(L, x - mu)    # L y = x − μ, তারপর d² = ‖y‖²
    z = (x - mu) / sd                 # শুধু আলাদা আলাদা standardize
    print(f"order {x}: Mahalanobis {np.sqrt(y @ y):.2f} | "
          f"per-feature z {np.sqrt(z @ z):.2f}")

fig, ax = plt.subplots(1, 2, figsize=(8, 3.4), layout="constrained")
ax[0].scatter(Z[:800, 0], Z[:800, 1], s=4, alpha=0.4)
ax[0].set_title("z ~ N(0, I)"); ax[0].set_aspect("equal")
ax[1].scatter(Xs[:800, 1], Xs[:800, 0], s=4, alpha=0.4)
ax[1].scatter(orders[:, 1], orders[:, 0], c="red", marker="x", s=60)
ax[1].set_title("x = mu + L z"); ax[1].set_xlabel("distance (km)")
ax[1].set_ylabel("delivery time (min)")
plt.show()

Sigma =
 [[100.  21.]
 [ 21.   9.]]
L =
 [[10.     0.   ]
 [ 2.1    2.142]]
sample mean: [34.903  5.956]
sample cov:
 [[97.87  20.927]
 [20.927  9.278]]
sample corr: 0.694
order [55. 12.]: Mahalanobis 2.17 | per-feature z 2.83
order [55.  3.]: Mahalanobis 3.91 | per-feature z 2.24

৭ · Mahalanobis দূরত্ব — covariance মেনে “কতটা অস্বাভাবিক”

দুটি নতুন অর্ডার, দুটিতেই ৫৫ মিনিট লেগেছে: একটি ১২ km দূরে, অন্যটি মাত্র ৩ km. কোনটি বেশি সন্দেহজনক? প্রতিটি feature আলাদা করে standardize করলে (z-score, AI-এর ভিত্তি কোর্সের পর্ব ৩৭) প্রথমটি দূরে দেখায় (২.৮৩ বনাম ২.২৪), কারণ দুটি feature-ই গড় থেকে দূরে। কিন্তু সাধারণ বুদ্ধি বলে উল্টো: দূরের অর্ডারে বেশি সময় স্বাভাবিক, কাছের অর্ডারে ৫৫ মিনিট অস্বাভাবিক। Mahalanobis distance এই সম্পর্ক মানে:

d_M(x)² = (x − μ)ᵀ Σ⁻¹ (x − μ)

এটি আবার একটি quadratic form, এবার Σ⁻¹ দিয়ে। Σ PD হলে Σ⁻¹-ও PD (eigenvalue 1/λ_i, সব ধনাত্মক), তাই d_M ≥ ০ আর শুধু x = μ-তে শূন্য, যা দূরত্বের জন্য জরুরি। Cholesky দিয়ে হিসাব সহজ: Σ⁻¹ = L^−ᵀL⁻¹, তাই d_M² = ‖y‖², যেখানে Ly = x − μ. inverse বানাতেই হয় না; ব্লকে ঠিক তাই করা হয়েছে। হাতে প্রথম অর্ডার: x − μ = (২০, ৬); y₁ = ২০/১০ = ২, y₂ = (৬ − ২.১ × ২)/২.১৪২ = ০.৮৪০; d = √(৪ + ০.৭০৬) = ২.১৭. দ্বিতীয়টি: x − μ = (২০, −৩); y₂ = (−৩ − ৪.২)/২.১৪২ = −৩.৩৬১; d = √(৪ + ১১.২৯) = ৩.৯১. y হলো “সাদা করা” (whitened) স্থানাঙ্ক: সেখানে মেঘটি আবার গোল, তাই সাধারণ দৈর্ঘ্যই ঠিক মাপ। Gaussian ধরলে d_M²-এর বণ্টন chi-square, যা দিয়ে outlier-এর সীমা টানা যায়; এ সবই পর্ব ১৯ (Multivariate Gaussian — covariance, Mahalanobis দূরত্ব আর শর্তসাপেক্ষ বণ্টন)-এ বিস্তারিত আসবে।

৮ · ML-এ কোথায় কোথায় — Hessian, ridge, kernel আর Newton

Hessian আর convexity. একটি loss function f(θ)-এর দ্বিতীয় অবকলনগুলোর matrix হলো Hessian H, যার (i, j) ঘর ∂²f/∂θ_i∂θ_j. একটি বিন্দু θ₀-এর কাছে second-order Taylor approximation: f(θ₀ + Δ) ≈ f(θ₀) + ∇fᵀΔ + ½ΔᵀHΔ. শেষ পদটি একটি quadratic form! তাই gradient শূন্য এমন বিন্দুতে (stationary point) H PD হলে সেখানে কঠোর local minimum (বাটির তলা), H indefinite হলে saddle point, আর H শুধু PSD হলে পরীক্ষা অমীমাংসিত (উপত্যকা, আরও উঁচু ক্রমের পদ দেখতে হয়)। আর যদি সব বিন্দুতে H PSD হয়, তাহলে f একটি convex function, যার যেকোনো local minimum-ই global. Gradient আর Hessian কীভাবে বের করতে হয়, তা পর্ব ১১ (Matrix calculus — gradient, Jacobian আর Hessian)-এ; convexity কেন optimization সহজ করে, তা পর্ব ১২ (Convexity — কেন কিছু loss সহজ আর কিছু কঠিন)-এ।

Least squares আর ridge. f(w) = ‖Xw − y‖²-এর Hessian 2XᵀX, যা অংশ ৫-এর যুক্তিতে সবসময় PSD, তাই least squares সবসময় convex. X-এর কলামগুলো স্বাধীন হলে XᵀX PD, আর minimum একটিই (AI-এর ভিত্তি কোর্সের পর্ব ৫৮-এর normal equations, যার সমাধানে Cholesky-ই প্রচলিত)। কলাম নির্ভরশীল হলে, যেমন অংশ ৫-এর মিনিট-ঘণ্টা জোড়া, XᵀX singular আর অসংখ্য সমাধান। Ridge যোগ করে λ‖w‖², Hessian হয় 2(XᵀX + λI), আর vᵀ(XᵀX + λI)v = ‖Xv‖² + λ‖v‖² > ০: λ > ০ হলে সবসময় PD. Ridge তাই আসলে নীতিগত একটি “jitter”, যার একটি পরিসংখ্যানিক অর্থও আছে (পর্ব ১৭-এ MAP হিসেবে)।

Logistic regression. এর log loss-এর Hessian XᵀWX, যেখানে W কর্ণে p_i(1 − p_i) > ০, তাই PSD আর loss convex. Newton পদ্ধতি প্রতিটি ধাপে HΔ = −∇f সমাধান করে; H PD হলে Cholesky দিয়ে, আর না হলে ধাপটি নিচে না নেমে ওপরে বা saddle-এর দিকে যেতে পারে। এজন্য অনেক বাস্তব optimizer H-এর কর্ণে কিছু যোগ করে (পর্ব ১৫, Second-order পদ্ধতি)।

Kernel আর Gaussian process. SVM-এর dual problem আর kernel ridge-এ Gram matrix K PSD হলেই সমস্যাটি convex থাকে (পর্ব ২৮ আর ২৯)। Gaussian process regression-এ (পর্ব ৪৬) K + σ²I-এর Cholesky দিয়েই অনুমান আর likelihood হিসাব হয়, আর σ² খুব ছোট হলে অংশ ৫-এর jitter ঠিক এখানে লাগে।

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

AI assistant এই বিষয়ে সূত্র ঠিকঠাক লেখে, কিন্তু শর্তগুলো প্রায়ই বাদ দেয়, আর ছোট সংখ্যার ভুল ধরা কঠিন কারণ লেখা আত্মবিশ্বাসী।

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

• Matrix প্রতিসম কি না assert করেছি, আর definiteness-এর দাবি eigvalsh দিয়ে মিলিয়েছি
• সবচেয়ে ছোট উদাহরণে (২ × ২) হাতে হিসাব করে AI-এর সংখ্যার সঙ্গে মিলিয়েছি, আর একটি পাল্টা উদাহরণ খুঁজেছি
• Sampling কোডে নমুনার গড় আর covariance লক্ষ্যের সঙ্গে তুলনা করেছি, কয়েকটি seed-এ
• Jitter থাকলে তার মান, কারণ আর ফলের ওপর প্রভাব লিখে রেখেছি
• Model-এর ফল দাবি করার আগে ডেটা ভাগ করেছি, test set তালা দেওয়া, আর একটি সরল baseline-এর সঙ্গে তুলনা করেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি