Eigenvalue আর eigenvector — ধারণা থেকে PCA পর্যন্ত

AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ NumPy দিয়ে matrix গুণ, solve আর eigenvalue-র প্রথম ঝলক দেখেছেন, পর্ব ৬৫-এ SVD দিয়ে PCA চালিয়েছেন, আর পর্ব ৫৩-এ stationary distribution তিন পথে বের করেছেন। এই কোর্সের পর্ব ০৬ (Vector space — span, basis, linear independence আর rank) basis আর null space দিয়েছে, আর পর্ব ০৭ (Matrix-এর জ্যামিতি — projection, rotation আর least squares) matrix-কে একটি জ্যামিতিক রূপান্তর হিসেবে দেখিয়েছে। এই পর্বের প্রশ্ন: একটি matrix কোন দিকগুলোকে ঘোরায় না, আর সেই দিকগুলো জানলে PCA, PageRank বা ক্রেতার দীর্ঘমেয়াদি আচরণ কেন হঠাৎ সরল হয়ে যায়? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড, আর scikit-learn প্রথমবার প্রায় ১০ সেকেন্ড লাগতে পারে।

১ · Eigenvector কী — যে দিক matrix শুধু টানে, ঘোরায় না

পর্ব ০৭-এ দেখেছেন, একটি বর্গাকার matrix A হলো একটি linear map: প্রতিটি vector x-কে নতুন vector Ax-এ পাঠায়। সাধারণত Ax আর x আলাদা দিকে থাকে। কিন্তু কিছু বিশেষ দিক থাকতে পারে, যেখানে A শুধু টানে বা চাপে (অথবা উল্টে দেয়), ঘোরায় না:

A v = λ v,   যেখানে v ≠ 0

এখানে A একটি n×n matrix, v একটি শূন্য নয় এমন vector যাকে বলে eigenvector, আর λ (ল্যাম্বডা) একটি সংখ্যা, A-এর eigenvalue. λ = ৩ মানে সেই দিকে সবকিছু তিনগুণ লম্বা; λ = ০.৫ মানে অর্ধেক; λ ঋণাত্মক হলে দিক উল্টে যায় কিন্তু একই রেখায় থাকে; λ = ০ হলে সেই দিক পুরো চেপে শূন্যে যায়, অর্থাৎ v পড়েছে A-এর null space-এ (পর্ব ০৬)। “Eigen” জার্মান শব্দ, মোটামুটি অর্থ “নিজস্ব”: matrix-এর নিজস্ব দিক আর নিজস্ব গুণক।

দুটি ছোট কিন্তু জরুরি পর্যবেক্ষণ। প্রথমত, v eigenvector হলে 3v বা −v-ও eigenvector, কারণ A(cv) = cAv = λ(cv). তাই eigenvector আসলে একটি দিক বা রেখা, একটি নির্দিষ্ট তীর নয়। সফটওয়্যার সাধারণত দৈর্ঘ্য ১ করে দেয় (unit vector), কিন্তু চিহ্ন (+ না −) যেকোনোটি হতে পারে; অংশ ৬ আর ৮-এ এটি বারবার ফিরবে। দ্বিতীয়ত, একই λ-এর সব eigenvector আর শূন্য vector মিলে একটি subspace গড়ে, যাকে বলে সেই λ-এর eigenspace.

নিচের stepper-এ A = [[4, 1], [2, 3]] (প্রথম সারি 4, 1; দ্বিতীয় সারি 2, 3). প্রতিটি ধাপে একটি v (ধূসর তীর) আর তার ছবি Av (নীল তীর) দেখুন; লাল আর সবুজ ড্যাশ রেখা দুটি eigen-দিক।

শেষ ধাপটিই পুরো পর্বের চাবি। একটি জটিল দেখতে matrix, সঠিক অক্ষে তাকালে, শুধু কয়েকটি স্বাধীন গুণ। AI-এর ভিত্তি কোর্সের পর্ব ৩৩-এ যে dot product আর norm শিখেছেন, সেগুলো দিয়েই এখন “দিক ঘুরল কি না” মাপা যায়: u আর Au-এর মধ্যের কোণ শূন্য বা ১৮০° হলে u eigen-দিক।

২ · Characteristic polynomial — 2×2 matrix হাতে সমাধান

Eigenvalue খোঁজার যুক্তি তিন লাইনে। Av = λv-কে লিখুন Av − λIv = 0, অর্থাৎ (A − λI)v = 0, যেখানে I হলো identity matrix. আমরা শূন্য নয় এমন v চাই। পর্ব ০৬-এর ভাষায়, (A − λI)-এর null space-এ শূন্য ছাড়া কিছু থাকতে হবে, অর্থাৎ matrix-টি singular. আর singular মানে determinant শূন্য:

det(A − λI) = 0

বাঁদিকটি λ-এর একটি বহুপদী (polynomial), যার নাম characteristic polynomial; n×n matrix-এ তার ঘাত n. A = [[a, b], [c, d]] হলে

det [[a − λ, b], [c, d − λ]] = (a − λ)(d − λ) − bc = λ² − (a + d)λ + (ad − bc)

এখানে a + d হলো কর্ণের যোগফল, যাকে বলে trace (tr A), আর ad − bc হলো det A. তাই 2×2-এর সূত্র: λ² − (tr A)λ + det A = 0. দ্বিঘাত সমীকরণের মূলের নিয়ম থেকে সঙ্গে সঙ্গে দুটি যাচাই আসে: λ₁ + λ₂ = tr A আর λ₁λ₂ = det A. এই দুটি সব n×n matrix-এর জন্যও সত্য (সব eigenvalue-র যোগ = trace, গুণ = determinant).

সব 2×2 matrix এত ভদ্র নয়। তিনটি বিশেষ ক্ষেত্র মনে রাখার মতো, কারণ পরীক্ষা আর ইন্টারভিউতে এগুলোই আসে:

এবার কম্পিউটারে। নিচের ব্লক একই A-এর characteristic polynomial-এর মূল np.roots দিয়ে আর eigenvalue np.linalg.eig দিয়ে বের করে, eigenvector যাচাই করে, তারপর একক বৃত্তের ২৪টি দিক নিয়ে মাপে A প্রতিটিকে কত ডিগ্রি ঘোরায়, আর ছবি আঁকে।

import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=4, suppress=True)
A = np.array([[4.0, 1.0],
              [2.0, 3.0]])
tr, det = np.trace(A), np.linalg.det(A)       # λ² − (trace)·λ + det = 0
print(f"trace {tr:.0f}, det {det:.4f}")
print("np.roots  :", np.roots([1, -tr, det]))
w, V = np.linalg.eig(A)                       # কলামগুলো eigenvector
print("eig dtype:", w.dtype, "| max |imag|:", np.abs(np.imag(w)).max())
w, V = np.real(w), np.real(V)                 # imag 0, তাই নিরাপদ
print("np.linalg.eig:", w)
print("V (columns):\n", V)
for lam, v in zip(w, V.T):
    print(f"λ = {lam:.0f}: A@v = {A @ v}  λ·v = {lam * v}")
t = np.linspace(0, 2 * np.pi, 24, endpoint=False)
U = np.vstack([np.cos(t), np.sin(t)])          # ২৪টি unit vector
AU = A @ U
cos = np.sum(U * AU, 0) / np.linalg.norm(AU, axis=0)
turn = np.degrees(np.arccos(np.clip(cos, -1, 1)))     # u আর Au-এর কোণ
print("turn (deg), first 12 directions:", turn[:12].round(1))
fig, ax = plt.subplots(figsize=(5.2, 5.2), layout="constrained")
for (x, y), (p, q) in zip(U.T, AU.T):
    ax.plot([x, p], [y, q], color="0.75", lw=0.8)
ax.plot(*U, "o", color="0.45", ms=3, label="unit vector u")
ax.plot(*AU, "o", color="C0", ms=4, label="image A u")
for lam, v, c in zip(w, V.T, ("C3", "C2")):
    ax.axline((0, 0), v, color=c, ls="--", lw=1)
    ax.annotate("", lam * v, (0, 0), arrowprops=dict(color=c, width=2))
    ax.plot([], [], color=c, label=f"eigen-direction, λ = {lam:.0f}")
ax.set_aspect("equal"); ax.set_xlim(-5.5, 5.5); ax.set_ylim(-5.5, 5.5)
ax.set_title("A stretches eigen-directions, turns the rest")
ax.legend(fontsize=8, loc="lower right")
plt.show()

trace 7, det 10.0000
np.roots  : [5. 2.]
eig dtype: float64 | max |imag|: 0.0
np.linalg.eig: [5. 2.]
V (columns):
 [[ 0.7071 -0.4472]
 [ 0.7071  0.8944]]
λ = 5: A@v = [3.5355 3.5355]  λ·v = [3.5355 3.5355]
λ = 2: A@v = [-0.8944  1.7889]  λ·v = [-0.8944  1.7889]
turn (deg), first 12 directions: [26.6 18.3  9.2  0.   8.5 15.4 18.4 13.3  5.4 26.6 34.5 32.8]

৩ · Diagonalisation — eigenvector-এর ভাষায় matrix লেখা

A-এর n-টি eigenvector-কে কলাম করে একটি matrix V বানান, আর eigenvalue-গুলো কর্ণে বসিয়ে Λ (বড় হাতের ল্যাম্বডা)। প্রতিটি কলামে Av_i = λ_iv_i, তাই একসঙ্গে লিখলে AV = VΛ. V-এর কলামগুলো linearly independent হলে V-এর inverse আছে, আর

A = V Λ V⁻¹

একে বলে diagonalisation বা eigen-decomposition. ডান থেকে বাঁয়ে পড়ুন: V⁻¹x দেয় x-এর স্থানাঙ্ক eigenvector-এর basis-এ (পর্ব ০৬-এর “basis-এ স্থানাঙ্ক”), Λ প্রতিটি স্থানাঙ্ককে নিজের λ দিয়ে গুণ করে, আর V ফল আবার সাধারণ অক্ষে ফিরিয়ে আনে। Stepper-এর শেষ ধাপে ঠিক এটিই হাতে করেছিলেন: (2, −1)-এর eigen-স্থানাঙ্ক (1, 1), গুণের পর (5, 2), ফেরত (7, 1).

ঘাত সস্তা হয়ে যায়। A² = VΛV⁻¹VΛV⁻¹ = VΛ²V⁻¹, আর সাধারণভাবে A^k = VΛ^kV⁻¹; Λ^k মানে শুধু কর্ণের প্রতিটি সংখ্যার k-তম ঘাত। আমাদের A-তে V = [[1, 1], [1, −2]], det V = −3, V⁻¹ = (1/3)[[2, 1], [1, −1]]. গুণ করলে

A^k = (1/3) · [[2·5^k + 2^k, 5^k − 2^k], [2·5^k − 2·2^k, 5^k + 2·2^k]]

k = ১ বসিয়ে যাচাই: (10 + 2)/3 = 4, (5 − 2)/3 = 1, (10 − 4)/3 = 2, (5 + 4)/3 = 3 ✓. k = ১০-এ 5¹⁰ = 9,765,625 আর 2¹⁰ = 1,024, তাই প্রথম ঘর (19,531,250 + 1,024)/3 = 6,510,758; np.linalg.matrix_power(A, 10)-ও ঠিক এই সংখ্যা দেয়। আর একটি জিনিস চোখে পড়ে: k বড় হলে 5^k 2^k-কে পুরো ছাপিয়ে যায়, তাই A^kx প্রায় পুরোপুরি λ = ৫-এর দিকে ঘুরে যায়। এটিই অংশ ৫-এর power iteration আর অংশ ৭-এর Markov chain-এর মূল কারণ।

কখন diagonalise করা যায়? যখন n-টি linearly independent eigenvector আছে। একটি যথেষ্ট শর্ত: n-টি eigenvalue সব আলাদা। দুটির জন্য ছোট প্রমাণ: ধরুন c₁v₁ + c₂v₂ = 0. A দিয়ে গুণ করলে c₁λ₁v₁ + c₂λ₂v₂ = 0; প্রথম সমীকরণকে λ₂ দিয়ে গুণ করে বিয়োগ করলে c₁(λ₁ − λ₂)v₁ = 0. λ₁ ≠ λ₂ আর v₁ ≠ 0, তাই c₁ = 0, তারপর c₂ = 0. Eigenvalue পুনরাবৃত্ত হলে কখনো হয় (3I), কখনো হয় না (shear). Defective matrix-এর জন্য Jordan form নামে একটি সাধারণীকরণ আছে, কিন্তু সংখ্যাগতভাবে তা অস্থির, আর ML-এ প্রায় দরকার পড়ে না, কারণ আমাদের গুরুত্বপূর্ণ matrix-গুলো বেশির ভাগই symmetric.

৪ · Symmetric matrix আর spectral theorem — লম্ব অক্ষের নিশ্চয়তা

A^T = A হলে matrix-টি symmetric: ঘর (i, j) আর (j, i) সমান। ML-এ সবচেয়ে গুরুত্বপূর্ণ matrix-গুলো এই জাতের: covariance matrix, loss-এর Hessian (পর্ব ১১), kernel বা Gram matrix (পর্ব ২৯), graph-এর Laplacian. এদের জন্য একটি অসাধারণ শক্তিশালী ফল আছে।

(খ)-এর ছোট প্রমাণ (আলাদা eigenvalue-র ক্ষেত্রে): Av₁ = λ₁v₁, Av₂ = λ₂v₂, λ₁ ≠ λ₂. তাহলে λ₁(v₁^Tv₂) = (Av₁)^Tv₂ = v₁^TA^Tv₂ = v₁^TAv₂ = λ₂(v₁^Tv₂). মাঝের ধাপে A^T = A ব্যবহার হয়েছে। তাই (λ₁ − λ₂)v₁^Tv₂ = 0, আর যেহেতু λ₁ ≠ λ₂, dot product v₁^Tv₂ = 0: দুটি orthogonal. পুনরাবৃত্ত eigenvalue-তে eigenspace-এর ভেতরে Gram–Schmidt (পর্ব ০৭) দিয়ে লম্ব basis বেছে নেওয়া যায়। (ক)-এর যুক্তি “আরও গভীরে” বাক্সে।

Rank-১ টুকরোর যোগ। QΛQ^T গুণ করে লিখলে A = λ₁q₁q₁^T + λ₂q₂q₂^T + … + λ_nq_nq_n^T. প্রতিটি q_iq_i^T হলো q_i-এর রেখায় orthogonal projection (পর্ব ০৭)। অর্থাৎ symmetric matrix-এর কাজ: x-কে লম্ব অক্ষগুলোতে ছায়া ফেলা, প্রতিটি ছায়াকে λ_i দিয়ে গুণ, তারপর যোগ। S-এর জন্য: 3 · ½[[1, 1], [1, 1]] + 1 · ½[[1, −1], [−1, 1]] = [[1.5 + 0.5, 1.5 − 0.5], [1.5 − 0.5, 1.5 + 0.5]] = [[2, 1], [1, 2]] ✓. ছোট λ-এর টুকরো ফেলে দিলে matrix-এর সবচেয়ে ভালো low-rank আন্দাজ মেলে (symmetric positive semidefinite matrix-এ); সাধারণ matrix-এ এর রূপ Eckart–Young, পর্ব ০৯-এ।

সব eigenvalue ধনাত্মক হলে x^TAx > 0 সব x ≠ 0-এ; এমন matrix-কে positive definite বলে। Covariance matrix-এর সব eigenvalue ≥ ০ কেন, আর quadratic form x^TAx-এর বাটি, জিন বা উপত্যকার আকার, সেটি পর্ব ১০ (Positive definite matrix — quadratic form আর covariance)-এর বিষয়। এখানে শুধু PCA-র জন্য যতটুকু লাগে, ততটুকু।

৫ · Power iteration — নিজের হাতে একটি eigen-solver

অংশ ৩-এ দেখেছেন A^kx সবচেয়ে বড় eigenvalue-র দিকে ঘুরে যায়। এই পর্যবেক্ষণকে algorithm বানালে পাই power iteration, যা Richard von Mises আর Hilda Pollaczek-Geiringer-এর ১৯২৯ সালের কাজের সঙ্গে সাধারণত যুক্ত করা হয়:

y = A x_k,   x_k+1 = y / ‖y‖,   λ̂ = x_k+1^T A x_k+1

প্রতি ধাপে গুণ, তারপর দৈর্ঘ্য ১ করা (না করলে সংখ্যা 5^k-এর মতো ফুলে overflow হয়)। শেষ রাশিটি Rayleigh quotient: একটি unit x-এর জন্য x^TAx, যা x eigenvector হলে ঠিক λ দেয়।

কেন কাজ করে। A diagonalisable হলে শুরুর x₀ = c₁v₁ + c₂v₂ + … + c_nv_n. তাহলে

A^kx₀ = λ₁^k [ c₁v₁ + c₂(λ₂/λ₁)^kv₂ + … + c_n(λ_n/λ₁)^kv_n ]

|λ₁| > |λ₂| ≥ … হলে বন্ধনীর বাকি পদ প্রতি ধাপে অন্তত |λ₂/λ₁| গুণ ছোট হয়। তাই দিকের ভুল কমে মোটামুটি |λ₂/λ₁|^k হারে। Symmetric A-তে Rayleigh quotient আরও ভালো: তার ভুল কমে প্রায় |λ₂/λ₁|^2k হারে, কারণ eigenvector-এর কাছে x^TAx প্রথম ক্রমে বদলায় না। দুটি শর্ত: (১) শুরুর x₀-এ v₁-এর অংশ থাকতে হবে (c₁ ≠ 0); এলোমেলো x₀ নিলে কার্যত সবসময় থাকে। (২) সবচেয়ে বড় |λ| একটিই হতে হবে; |λ₁| = |λ₂| হলে (যেমন λ = ১ আর −১) x চিরকাল দোলে। 1 − |λ₂/λ₁|-এর মতো মাপকে সাধারণভাবে বলে eigenvalue gap: ফাঁক যত বড়, convergence তত দ্রুত।

এবার ৪×৪ matrix-এ, দুটি ভিন্ন gap নিয়ে। দুটি matrix-এর সবচেয়ে বড় eigenvalue ৫ আর eigenvector-ও এক (একই orthogonal Q দিয়ে বানানো, QR থেকে, যা পর্ব ০৭-এর বিষয়), পার্থক্য শুধু λ₂: একটিতে ৪.৫ (অনুপাত ০.৯), অন্যটিতে ১ (অনুপাত ০.২)। শুরু একই এলোমেলো x₀, আর শেষে np.linalg.eigh-এর উত্তরের সঙ্গে মেলানো।

import numpy as np

rng = np.random.default_rng(42)
Q, _ = np.linalg.qr(rng.standard_normal((4, 4)))   # একটি এলোমেলো orthogonal matrix

def make(lams):                                     # জানা eigenvalue-এর symmetric matrix
    return Q @ np.diag(lams) @ Q.T

def power_iteration(A, steps, x):
    rows = []
    for k in range(1, steps + 1):
        y = A @ x                                   # ১. গুণ
        x = y / np.linalg.norm(y)                   # ২. দৈর্ঘ্য 1 করা
        lam = x @ A @ x                             # ৩. Rayleigh quotient
        rows.append((k, lam, x))
    return rows

x0 = rng.standard_normal(4)
for lams in ([5, 4.5, 2, 1], [5, 1, 0.5, 0.2]):
    A = make(lams)
    w, V = np.linalg.eigh(A)                        # উত্তর মেলানোর জন্য
    v1 = V[:, -1]                                   # eigh: ছোট থেকে বড়, শেষটি সবচেয়ে বড়
    print(f"eigenvalues {lams}, ratio λ2/λ1 = {lams[1] / lams[0]}")
    for k, lam, x in power_iteration(A, 30, x0):
        if k in (1, 2, 3, 5, 10, 20, 30):
            ang = np.degrees(np.arccos(min(1.0, abs(x @ v1))))
            print(f"  step {k:>2}: λ estimate {lam:.6f} | error {5 - lam:.1e}"
                  f" | angle to v1 {ang:.4f}°")
    print("  eigh top eigenvalue:", round(w[-1], 6))

eigenvalues [5, 4.5, 2, 1], ratio λ2/λ1 = 0.9
  step  1: λ estimate 4.775344 | error 2.2e-01 | angle to v1 30.3040°
  step  2: λ estimate 4.888443 | error 1.1e-01 | angle to v1 26.2928°
  step  3: λ estimate 4.916856 | error 8.3e-02 | angle to v1 23.7404°
  step  5: λ estimate 4.943883 | error 5.6e-02 | angle to v1 19.5632°
  step 10: λ estimate 4.978917 | error 2.1e-02 | angle to v1 11.8497°
  step 20: λ estimate 4.997338 | error 2.7e-03 | angle to v1 4.1842°
  step 30: λ estimate 4.999675 | error 3.3e-04 | angle to v1 1.4612°
  eigh top eigenvalue: 5.0
eigenvalues [5, 1, 0.5, 0.2], ratio λ2/λ1 = 0.2
  step  1: λ estimate 4.930570 | error 6.9e-02 | angle to v1 7.4884°
  step  2: λ estimate 4.997573 | error 2.4e-03 | angle to v1 1.4078°
  step  3: λ estimate 4.999906 | error 9.4e-05 | angle to v1 0.2772°
  step  5: λ estimate 5.000000 | error 1.5e-07 | angle to v1 0.0110°
  step 10: λ estimate 5.000000 | error 1.7e-14 | angle to v1 0.0000°
  step 20: λ estimate 5.000000 | error 2.7e-15 | angle to v1 0.0000°
  step 30: λ estimate 5.000000 | error 2.7e-15 | angle to v1 0.0000°
  eigh top eigenvalue: 5.0

বাস্তবে কী চলে। শুদ্ধ power iteration শুধু একটি eigenvector দেয়, আর gap ছোট হলে ধীর। বাস্তব সফটওয়্যার এর উন্নত রূপ ব্যবহার করে: বাকিগুলোর জন্য deflation বা একসঙ্গে কয়েকটি vector (subspace iteration); কোনো নির্দিষ্ট মানের কাছের eigenvalue-র জন্য inverse iteration ((A − σI)⁻¹ দিয়ে গুণ); আর বড় sparse matrix-এ Lanczos পদ্ধতি (১৯৫০), যা SciPy-র eigsh ARPACK-এর মাধ্যমে চালায়। ছোট dense matrix-এর সব eigenvalue-র জন্য LAPACK-এর QR-ভিত্তিক routine (John Francis আর Vera Kublanovskaya, ১৯৬১ সালের কাছাকাছি স্বাধীনভাবে)। নিচে একটি ২০০০×২০০০ covariance matrix-এর সবচেয়ে বড় তিনটি eigenvalue তিন পথে, ডেস্কটপে সত্যিকারের run (Python 3.14.4, NumPy 2.5.4, SciPy 1.18.1; সময় একটি ল্যাপটপের একটি run, সাধারণ ফল নয়):

~/bazarghar-ml$ python eig_check.py
…
numpy eigh (all 2000) : [4.719  4.6754 4.6016]  0.54 s
scipy eigh subset (3) : [4.719  4.6754 4.6016]  0.29 s
eigsh Lanczos (3)     : [4.719  4.6754 4.6016]  0.07 s

তিনটির উত্তর এক; সব ২০০০টি eigenvalue না চাইলে scipy.linalg.eigh(..., subset_by_index=...) বা eigsh অনেক কম কাজ করে। Matrix যত বড় আর sparse, পার্থক্য তত বেশি।

৬ · Covariance-এর eigen থেকে PCA — variance সর্বোচ্চ করার derivation

AI-এর ভিত্তি কোর্সের পর্ব ৬৫-এ PCA চালিয়েছেন SVD দিয়ে, আর explained variance মেপেছেন। সেখানে “সবচেয়ে বেশি variance-এর দিক”-এর উত্তর হিসেবে eigenvector এসেছিল। এখানে প্রশ্নটি উল্টো দিক থেকে: কেন সেই দিক eigenvector-ই হবে?

সেটআপ। X একটি n×d ডেটা matrix (n গ্রাহক, d feature). প্রতিটি কলাম থেকে গড় বাদ দিলে পাই centred X_c. তার covariance matrix

Σ = X_c^TX_c / (n − 1)

একটি d×d symmetric matrix: কর্ণে প্রতিটি feature-এর variance, বাইরে প্রতি জোড়ার covariance. এখন একটি unit দিক u (‖u‖ = ১) নিন, আর প্রতিটি গ্রাহককে সেই দিকে ছায়া ফেলুন: z = X_cu, n-টি সংখ্যা। তাদের variance

Var(z) = z^Tz / (n − 1) = u^TX_c^TX_cu / (n − 1) = u^TΣu

অর্থাৎ যেকোনো দিকে ছড়ানো = একটি quadratic form. PCA-র প্রথম অক্ষের সমস্যা তাই: ‖u‖ = ১ শর্তে u^TΣu সর্বোচ্চ করুন।

প্রমাণ ২, Lagrange দিয়ে (পর্ব ১৬-এর হাতিয়ারের প্রথম ঝলক): L(u, λ) = u^TΣu − λ(u^Tu − 1). u-এর সাপেক্ষে gradient 2Σu − 2λu (quadratic form-এর gradient পর্ব ১১-এ derive হবে); শূন্য বসালে Σu = λu. অর্থাৎ শর্তসাপেক্ষ সর্বোচ্চ বা সর্বনিম্নের প্রতিটি সম্ভাব্য বিন্দুই eigenvector, আর সেখানে u^TΣu = λu^Tu = λ, তাই সবচেয়ে বড়টি বেছে নিন। এখানে λ হলো Lagrange multiplier, যা শেষে eigenvalue হয়ে দাঁড়ায়। দ্বিতীয় PC: u ⟂ q₁ শর্ত যোগ করলে a₁ = ০, আর একই যুক্তিতে সর্বোচ্চ λ₂, q₂-এ। এভাবে সব PC পরস্পর লম্ব, আর j-তম PC-র ভাগ explained variance ratio = λ_j / (λ₁ + … + λ_d); হর হলো tr Σ, মোট variance.

SVD-এর সঙ্গে সম্পর্ক। X_c = USV^T হলে Σ = X_c^TX_c/(n − 1) = V(S²/(n − 1))V^T. এটি হুবহু spectral decomposition: SVD-এর ডান singular vector = Σ-এর eigenvector, আর λ_j = s_j²/(n − 1). তাহলে কোনটি ব্যবহার করবেন? Covariance বানালে condition number বর্গ হয়ে যায়, তাই খুব ছোট-বড় মেশানো singular value থাকলে SVD বেশি নির্ভুল। আবার n অনেক বড় আর d ছোট হলে d×d covariance-এর eigh অনেক সস্তা। scikit-learn 1.5 থেকে PCA-তে svd_solver="covariance_eigh" যোগ হয়েছে, আর docs অনুযায়ী "auto" নীতিতে feature ১০০০-এর কম আর sample তার ১০ গুণের বেশি হলে এটিই বাছা হয়। অর্থাৎ নিচের উদাহরণে scikit-learn নিজেও ভেতরে covariance-এর eigen-decomposition চালাচ্ছে।

নিচের ব্লকে বাজারঘরের ৫০০ জন কাল্পনিক গ্রাহক, চারটি feature: মাসে অর্ডার, গড় বাস্কেট (৳), শেষ অর্ডারের পর কত দিন (recency), আর রিভিউ সংখ্যা। চারটিই একটি লুকানো “আগ্রহ” থেকে কিছুটা আসে, তাই পরস্পর সম্পর্কিত। এককগুলো আলাদা বলে আগে standardize করা হয়েছে (AI-এর ভিত্তি কোর্সের পর্ব ৬৫-এর scaling-এর ফাঁদ মনে করুন)। তারপর তিন পথে PCA: নিজের eigh, scikit-learn, আর SVD.

import numpy as np
from sklearn.decomposition import PCA
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(42)                     # কাল্পনিক বাজারঘর: ৫০০ গ্রাহক
n = 500
engage = rng.normal(0, 1, n)                        # লুকানো "আগ্রহ"
orders = 4 + 1.5 * engage + rng.normal(0, 0.8, n)              # মাসে অর্ডার
basket = 900 + 150 * engage + rng.normal(0, 200, n)            # গড় বাস্কেট (৳)
recency = 20 - 6 * engage + rng.normal(0, 5, n)                # শেষ অর্ডারের পর দিন
reviews = 1 + 0.2 * engage + rng.normal(0, 1, n)               # রিভিউ সংখ্যা
X = np.column_stack([orders, basket, recency, reviews])
Z = (X - X.mean(0)) / X.std(0)                      # standardize (একক মুছে ফেলা)
S = np.cov(Z, rowvar=False)                         # Σ, ভাগ n − 1
lam, V = np.linalg.eigh(S)                          # ছোট থেকে বড়
lam, V = lam[::-1], V[:, ::-1]                      # বড় থেকে ছোট
print("eigenvalues of Σ :", lam)
print("explained ratio  :", lam / lam.sum())
pca = PCA().fit(Z)
print("sklearn variance :", pca.explained_variance_)
_, s, Vt = np.linalg.svd(Z - Z.mean(0), full_matrices=False)
print("SVD s²/(n−1)     :", s**2 / (n - 1))
sign = np.sign(np.sum(V * pca.components_.T, 0))   # চিহ্ন মেলানো
print("same directions (signs aside):",
      np.allclose(V * sign, pca.components_.T), np.allclose(np.abs(Vt), np.abs(V.T)))
print("PC1 weights (orders, basket, recency, reviews):", V[:, 0] * sign[0])
v1 = V[:, 0]
print("v1ᵀΣv1 =", round(v1 @ S @ v1, 4), "| var of scores =", round((Z @ v1).var(ddof=1), 4))
for name, u in [("orders only", np.eye(4)[0]), ("random unit", rng.normal(size=4))]:
    u = u / np.linalg.norm(u)
    print(f"uᵀΣu for {name:<11}: {u @ S @ u:.4f}  (≤ λ1 = {lam[0]:.4f})")

eigenvalues of Σ : [2.2007 0.9335 0.5661 0.3077]
explained ratio  : [0.5491 0.2329 0.1413 0.0768]
sklearn variance : [2.2007 0.9335 0.5661 0.3077]
SVD s²/(n−1)     : [2.2007 0.9335 0.5661 0.3077]
same directions (signs aside): True True
PC1 weights (orders, basket, recency, reviews): [ 0.5934  0.5041 -0.5775  0.2457]
v1ᵀΣv1 = 2.2007 | var of scores = 2.2007
uᵀΣu for orders only: 1.0020  (≤ λ1 = 2.2007)
uᵀΣu for random unit: 0.4485  (≤ λ1 = 2.2007)

৭ · Markov chain-এর stationary distribution — eigenvalue ১-এর eigenvector

AI-এর ভিত্তি কোর্সের পর্ব ৫৩-এ বাজারঘরের ক্রেতার চারটি অবস্থার (Active, Occasional, Dormant, Churned) একটি Markov chain ছিল, আর transition matrix P-এর stationary distribution π (যেখানে πP = π) তিন পথে বের হয়েছিল, একটি পথ eig(P.T). সেখানে যা “কেন” প্রশ্ন বাকি ছিল, এখানে তার উত্তর।

১. eigenvalue ১ কেন সবসময় আছে? P-এর প্রতিটি সারির যোগ ১. তাই সব-১ vector 1-এর জন্য P1 = 1: ১ হলো P-এর eigenvalue, ডান eigenvector 1. আর det(P^T − λI) = det((P − λI)^T) = det(P − λI), তাই P আর P^T-এর eigenvalue একই; P^T-এরও eigenvalue ১, আর তার eigenvector-ই π^T (πP = π-কে transpose করলে P^Tπ^T = π^T).

২. কোনো eigenvalue ১-এর চেয়ে বড় নয় কেন? Px-এর প্রতিটি ঘর x-এর ঘরগুলোর একটি ওজনদার গড় (ওজন অঋণাত্মক, যোগ ১), তাই max|(Px)_i| ≤ max|x_i|. Px = λx হলে |λ| · max|x_i| ≤ max|x_i|, অর্থাৎ |λ| ≤ ১.

৩. π একটিই আর ধনাত্মক কেন? এর উত্তর Perron (১৯০৭) আর Frobenius (১৯১২)-এর নামে পরিচিত Perron–Frobenius উপপাদ্য: সব ঘর ধনাত্মক এমন matrix-এ সবচেয়ে বড় eigenvalue একক (পুনরাবৃত্ত নয়), বাকি সবার |λ| তার চেয়ে কঠোরভাবে ছোট, আর তার eigenvector-এর সব ঘর ধনাত্মক নেওয়া যায়। AI-এর ভিত্তি কোর্সের পর্ব ৫৩-এর irreducible আর aperiodic শর্ত ঠিক এই নিশ্চয়তা পাওয়ার শর্ত।

৪. কত দ্রুত? অংশ ৩-এর মতোই: শুরুর বণ্টন x₀-কে বাঁ eigenvector-এ ভাঙলে x_n = x₀Pⁿ = π + (λ₂ⁿ-এর অনুপাতে কমতে থাকা পদ) + …. তাই π থেকে দূরত্ব মোটামুটি |λ₂|ⁿ হারে কমে, আর 1 − |λ₂| হলো chain-এর spectral gap. বাজারঘরের chain-এ |λ₂| = ০.৭২৫১: প্রতি মাসে দূরত্ব প্রায় ২৭% কমে। এভাবে দেখলে x ← xP পুনরাবৃত্তিটি আসলে P^T-এর ওপর power iteration, যার “সবচেয়ে বড়” eigenvalue ১. PageRank-এর ক্ষেত্রেও একই কথা (AI-এর ভিত্তি কোর্সের পর্ব ৫৩-এ দেখেছেন): damping factor c (প্রায়ই ০.৮৫) যোগ করলে Google matrix-এর দ্বিতীয় eigenvalue-র মান c-এর বেশি হয় না বলে Haveliwala আর Kamvar ২০০৩ সালে দেখান, তাই কোটি পাতার web graph-এও power iteration নিশ্চিতভাবে এগোয়।

নিচের ব্লকের প্রথম অংশ এই চারটি কথা সংখ্যায় যাচাই করে। দ্বিতীয় অংশ পরের অংশের জন্য সতর্কবার্তা: symmetric-এর জন্য বানানো eigh-কে non-symmetric P^T দিলে কী হয়, আর প্রায়-symmetric matrix-এ eig কী করে।

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

P = np.array([[0.70, 0.20, 0.08, 0.02],      # Active, Occasional, Dormant, Churned
              [0.30, 0.40, 0.25, 0.05],      # সারি = এই মাস, কলাম = পরের মাস
              [0.10, 0.20, 0.50, 0.20],
              [0.10, 0.05, 0.05, 0.80]])
print("P @ ones =", P @ np.ones(4))           # সারির যোগ 1 → eigenvalue 1
w, V = np.linalg.eig(P.T)                     # πP = π  ⇔  Pᵀπᵀ = πᵀ
order = np.argsort(-np.abs(w))
w, V = np.real(w[order]), np.real(V[:, order])
pi = V[:, 0] / V[:, 0].sum()
print("eigenvalues:", w, "| π:", pi)
x = np.array([1.0, 0, 0, 0])                  # Active থেকে শুরু
for n in range(1, 31):
    x = x @ P
    if n in (5, 10, 20, 30):
        print(f"n={n:>2}: ‖x−π‖₁ = {np.abs(x - pi).sum():.2e} | |λ2|^n = {abs(w[1])**n:.2e}")
print("eigh on non-symmetric Pᵀ (wrong!):", np.linalg.eigh(P.T)[0])

S = np.array([[2.0, 1e-9], [-1e-9, 2.0]])     # প্রায় symmetric, eigenvalue জোড়া
we = np.linalg.eig(S)[0]
print("eig :", we.real, "| imag:", [f"{v:.0e}" for v in we.imag], we.dtype)
print("eigh:", np.linalg.eigh(S)[0], np.linalg.eigh(S)[0].dtype)
C = np.array([[4.0, 2.0], [2.0 + 1e-6, 3.0]])  # রাউন্ডিংয়ে সামান্য অসমান
L = np.tril(C) + np.tril(C, -1).T              # নিচের ত্রিভুজ থেকে বানানো
print("eig − eigh:", f"{np.abs(np.sort(np.linalg.eig(C)[0].real) - np.linalg.eigh(C)[0]).max():.1e}",
      "| eigh(C) == eigh(lower):", np.array_equal(np.linalg.eigh(C)[0], np.linalg.eigh(L)[0]))

P @ ones = [1. 1. 1. 1.]
eigenvalues: [1.     0.7251 0.4979 0.177 ] | π: [0.3499 0.1998 0.1827 0.2676]
n= 5: ‖x−π‖₁ = 2.00e-01 | |λ2|^n = 2.00e-01
n=10: ‖x−π‖₁ = 4.33e-02 | |λ2|^n = 4.02e-02
n=20: ‖x−π‖₁ = 1.76e-03 | |λ2|^n = 1.61e-03
n=30: ‖x−π‖₁ = 7.09e-05 | |λ2|^n = 6.48e-05
eigh on non-symmetric Pᵀ (wrong!): [0.1615 0.4676 0.7614 1.0095]
eig : [2. 2.] | imag: ['1e-09', '-1e-09'] complex128
eigh: [2. 2.] float64
eig − eigh: 4.9e-07 | eigh(C) == eigh(lower): True

৮ · সংখ্যার সতর্কতা — eig বনাম eigh, ক্রম, চিহ্ন আর version

হাতে eigenvalue বের করা আর কম্পিউটারে বের করার মধ্যে কয়েকটি ফাঁদ আছে, যেগুলো কোনো error দেয় না, শুধু চুপচাপ ভুল বা অপ্রত্যাশিত উত্তর দেয়।

NumPy 2.5-এর পরিবর্তন। NumPy 2.5.0 (জুন ২০২৬)-এর release notes-এ compatibility note: linalg.eig আর linalg.eigvals এখন সবসময় complex array দেয়, আগে ফলের dtype মানের ওপর নির্ভর করত (gh-30411); notes নিজেই symmetric matrix-এর জন্য eigh/eigvalsh-এর পরামর্শ দেয়। এই পাতার ব্রাউজার-Python-এ NumPy 2.4.6, তাই অংশ ২-এর ব্লকে dtype float64 এসেছে। ডেস্কটপে একই A, সত্যিকারের run:

~/bazarghar-ml$ python eig_check.py
Python 3.14.4 | NumPy 2.5.4 | SciPy 1.18.1
eig(A): [5.+0.j 2.+0.j] complex128
eigvalsh(A @ A.T) dtype: float64
…

পুরোনো কোডে w.dtype == float ধরে নেওয়া থাকলে, বা ফল সরাসরি pandas-এ বা JSON-এ লেখা হলে, NumPy আপগ্রেডের পর তা ভাঙতে পারে। নিরাপদ অভ্যাস: symmetric হলে eigh; না হলে imag অংশ সত্যিই শূন্যের কাছে কি না যাচাই করে তারপর np.real, যেমন অংশ ২ আর ৭-এর কোডে।

Non-symmetric matrix-এর সংবেদনশীলতা। J = [[1, 1], [ε, 1]] নিন। tr = ২, det = 1 − ε, তাই λ = 1 ± √ε. ε = 10⁻¹⁰, প্রায় রাউন্ডিংয়ের মাপের একটি বদল, eigenvalue-কে 10⁻⁵ সরিয়ে দেয়: পরিবর্তনের চেয়ে লক্ষ গুণ বড় ফল। Symmetric matrix-এ এমন হয় না: Weyl-এর অসমতা অনুযায়ী E যোগ করলে প্রতিটি eigenvalue সর্বোচ্চ ‖E‖ (spectral norm) সরে। অংশ ৭-এর শেষ লাইনে C-এর অসমতা ১০^−৬, আর eig ও eigh-এর পার্থক্য তারও কম।

চিহ্ন আর পুনরাবৃত্ত eigenvalue: v আর −v দুটিই বৈধ উত্তর, তাই দুটি লাইব্রেরি বা দুটি version আলাদা চিহ্ন দিতে পারে; PCA-র loading-এ “recency-র ওজন ঋণাত্মক” এমন ব্যাখ্যা দেওয়ার আগে চিহ্ন একটি নিয়মে বেঁধে নিন (যেমন সবচেয়ে বড় |ওজন| ধনাত্মক)। আর λ₁ ≈ λ₂ হলে সেই দুটি eigenvector আলাদাভাবে প্রায় অর্থহীন: তাদের span স্থির, কিন্তু span-এর ভেতরের দিক ডেটার সামান্য বদলেই ঘুরে যায় (Davis–Kahan). ভাবনার প্রশ্ন ৫-এ এর ব্যবহারিক প্রভাব।

Covariance নিজেই একটু অসমান হলে। np.cov বা X^TX থেকে বানানো matrix সাধারণত symmetric, কিন্তু নিজে লেখা কোডে (যেমন দুটি আলাদা গণনা থেকে উপরের আর নিচের ত্রিভুজ) রাউন্ডিংয়ের অসমতা ঢুকতে পারে। তখন (C + C.T) / 2 দিয়ে symmetric করে নিয়ে eigh চালান, আর ডেটা থেকে PCA দরকার হলে সরাসরি SVD বা scikit-learn ব্যবহার করুন।

৯ · AI-এর যুগে — eigen আর PCA-র কোড ও derivation যাচাই

AI assistant linear algebra-র derivation খুব গুছিয়ে লেখে, আর চিহ্নগুলো নিখুঁত দেখায়। কিন্তু eigen-এর কোডে তার ভুলগুলো প্রায়ই নীরব: কোনো error নেই, শুধু ভুল দিক।

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

• প্রতিটি eigenpair-এ np.allclose(A @ v, lam * v) চালিয়েছি
• Matrix symmetric কি না যাচাই করে তবেই eigh নিয়েছি, আর ফলের ক্রম (ascending) মাথায় রেখেছি
• eigenvalue-র যোগ trace-এর সঙ্গে আর গুণ determinant-এর সঙ্গে মিলিয়েছি (ছোট matrix-এ)
• নিজের PCA-র explained variance scikit-learn আর SVD-এর সঙ্গে মিলিয়েছি, চিহ্ন বাদে
• ডেটা আগে ভাগ করেছি, scaler আর PCA শুধু train-এ fit করেছি, test set তালা দেওয়া
• PCA-র পরে model-এর ফল কয়েকটি seed বা fold-এ, PCA ছাড়া একটি সরল baseline-এর সঙ্গে তুলনা করেছি
• নিজের NumPy আর scikit-learn version-এর অফিসিয়াল docs পড়েছি; কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি