SVD — low-rank approximation, pseudo-inverse আর recommender-এর ভিত্তি

AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ np.linalg.solve, lstsq আর condition number-এর প্রথম পরিচয় হয়েছিল, পর্ব ৬৫-এ SVD দিয়ে PCA, আর পর্ব ৬৯-এ truncated SVD দিয়ে LSA. সেখানে SVD ছিল একটি কাজের যন্ত্র। এই পর্বে প্রশ্ন: যন্ত্রটি ভেতরে কী করে, কেন সবসময় কাজ করে, আর কেন ML-এর এত জায়গায় (compression, least squares, text, recommender) ঘুরেফিরে আসে? এই মডিউলের আগের পর্বগুলো (পর্ব ০৬-এ rank, পর্ব ০৭-এ orthogonal matrix আর projection, পর্ব ০৮-এ eigen আর spectral theorem) এখানে একসঙ্গে কাজে লাগবে। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, matplotlib আর scikit-learn); scikit-learn প্রথমবার নামতে প্রায় ১০ সেকেন্ড লাগে।

১ · SVD-এর ছবি — ঘোরানো, টানা, আবার ঘোরানো

একটি m × n matrix A একটি n-মাত্রার vector x নিয়ে m-মাত্রার Ax দেয়। Singular value decomposition বলে, A-কে সবসময় লেখা যায়

A = U Σ Vᵀ

যেখানে U একটি m × m orthogonal matrix (কলামগুলো পরস্পর লম্ব আর একক দৈর্ঘ্যের, তাই UᵀU = I), V একটি n × n orthogonal matrix, আর Σ একটি m × n “কর্ণ” matrix, যার কর্ণে σ₁ ≥ σ₂ ≥ … ≥ 0 বসে আর বাকি সব ঘর শূন্য। Vᵀ মানে V-এর transpose. σ_i-গুলো singular value, V-এর কলাম v_i হলো right singular vector (ইনপুটের জগতে বিশেষ দিক), U-এর কলাম u_i left singular vector (আউটপুটের জগতে বিশেষ দিক)। Orthogonal matrix দৈর্ঘ্য আর কোণ বদলায় না, তাই সে কেবল ঘোরায় বা আয়নায় উল্টায় (পর্ব ০৭)।

এক লাইনে মূল সম্পর্ক: A v_i = σ_i u_i. অর্থাৎ ইনপুটের i-তম বিশেষ দিককে A আউটপুটের i-তম বিশেষ দিকে পাঠায়, দৈর্ঘ্য σ_i গুণ করে। Eigenvector-এর সঙ্গে তুলনা করুন (পর্ব ০৮): সেখানে Av = λv, একই দিক। SVD দিক বদলাতে দেয় (v থেকে u), বিনিময়ে সবসময় থাকে, যেকোনো আকারের যেকোনো matrix-এর জন্য, আর σ কখনো ঋণাত্মক বা জটিল সংখ্যা হয় না। Eigen-decomposition-এর এই নিশ্চয়তা নেই।

নিচের stepper A = [[3, 0], [4, 5]] নিয়ে একক বৃত্তকে তিন ধাপে ellipse বানায়। এখানে V আর U দুটিই ঘূর্ণন হিসেবে বেছে নেওয়া (অংশ ২-এ দেখবেন NumPy আয়নাসহ অন্য একটি সমান সঠিক রূপ দেয়)।

Matrix বর্গাকার না হলে Σ-ও বর্গাকার নয়। ML-এ ডেটা matrix প্রায় সবসময় “লম্বা” (n সারি উদাহরণ, d কলাম feature, n ≫ d), তাই full SVD-এর U-এর অনেক কলাম শূন্য σ-এর সঙ্গে গুণ হয়ে কাজে লাগে না। ব্যবহারিক রূপ thin SVD: শুধু প্রথম min(m, n)-টি কলাম রাখা। NumPy-তে np.linalg.svd(A, full_matrices=False) এটাই দেয়। আর শুধু প্রথম k-টি রাখলে truncated SVD, যা অংশ ৪ থেকে আমাদের মূল বিষয়।

২ · AᵀA আর AAᵀ থেকে SVD — হাতে একটি 2×2

SVD কোথা থেকে আসে, তা দুই লাইনে বের করা যায়। A = UΣVᵀ হলে

AᵀA = VΣᵀUᵀUΣVᵀ = V (ΣᵀΣ) Vᵀ,   AAᵀ = U (ΣΣᵀ) Uᵀ

কারণ UᵀU = I. ΣᵀΣ একটি কর্ণ matrix, কর্ণে σ_i². তাই V-এর কলাম হলো AᵀA-র eigenvector, আর σ_i² তার eigenvalue; একইভাবে U-এর কলাম AAᵀ-এর eigenvector. উল্টো দিক থেকে এটিই অস্তিত্বের প্রমাণের নকশা: AᵀA প্রতিসম (symmetric), তাই পর্ব ০৮-এর spectral theorem অনুযায়ী তার orthonormal eigenvector আছে, আর xᵀAᵀAx = ‖Ax‖² ≥ 0 বলে eigenvalue কখনো ঋণাত্মক নয়, তাই বর্গমূল নেওয়া যায়। σ_i = √λ_i ধরে, σ_i > 0 হলে u_i = Av_i/σ_i সংজ্ঞায়িত করুন। এগুলো একক দৈর্ঘ্যের আর পরস্পর লম্ব, কারণ u_iᵀu_j = v_iᵀAᵀAv_j/(σ_iσ_j) = λ_jv_iᵀv_j/(σ_iσ_j), যা i ≠ j হলে ০ আর i = j হলে ১. প্রয়োজনে আরও লম্ব vector যোগ করে U সম্পূর্ণ করা হয়। ব্যস, যেকোনো বাস্তব matrix-এর SVD আছে।

হাতে একটি পূর্ণ উদাহরণ (পরীক্ষায় প্রায়ই আসে এমন ধরনের)। A = [[3, 0], [4, 5]].

1. AᵀA = [[3·3 + 4·4, 3·0 + 4·5], [20, 0·0 + 5·5]] = [[25, 20], [20, 25]].
2. Characteristic polynomial: λ² − (trace)λ + det = λ² − 50λ + (625 − 400) = λ² − 50λ + 225 = (λ − 45)(λ − 5). তাই λ₁ = ৪৫, λ₂ = ৫, আর σ₁ = √45 = 3√5 ≈ ৬.৭০৮, σ₂ = √5 ≈ ২.২৩৬.
3. λ = ৪৫-এর eigenvector: (25 − 45)x + 20y = 0 → y = x, তাই v₁ = (1, 1)/√2. λ = ৫-এর: 20x + 20y = 0 → v₂ = (−1, 1)/√2.
4. u₁ = Av₁/σ₁ = (3, 9)/√2 ÷ √45 = (3, 9)/√90 = (1, 3)/√10 ≈ (0.316, 0.949). u₂ = Av₂/σ₂ = (−3, 1)/√2 ÷ √5 = (−3, 1)/√10.
5. যাচাই: σ₁σ₂ = √225 = ১৫ = |det A|, আর σ₁/σ₂ = ৩, যা A-এর condition number (অংশ ৭)।

নিচের ব্লক NumPy-র SVD হাতের হিসাবের সঙ্গে মেলায়, AᵀA-র eigen দিয়ে V যাচাই করে, তারপর stepper-এর চার ধাপ ছবিতে আঁকে।

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

A = np.array([[3., 0.], [4., 5.]])
U, s, Vt = np.linalg.svd(A)
print("singular values:", s)
print("U =\n", U, "\nVt =\n", Vt)
print("A == U Σ Vt ?", np.allclose(A, U @ np.diag(s) @ Vt))

lam, Q = np.linalg.eigh(A.T @ A)          # ছোট থেকে বড় ক্রমে আসে
print("eig(AᵀA):", lam, "→ √ =", np.sqrt(lam[::-1]))
print("V-র কলাম == AᵀA-র eigenvector (চিহ্ন বাদে)?",
      np.allclose(np.abs(Vt.T), np.abs(Q[:, ::-1])))
print("u1 == A v1 / σ1 ?", np.allclose(U[:, 0], A @ Vt[0] / s[0]))
print("det(A) =", round(np.linalg.det(A), 4), "| σ1·σ2 =", round(s.prod(), 4))

t = np.linspace(0, 2 * np.pi, 200)
C = np.vstack([np.cos(t), np.sin(t)])          # একক বৃত্তের ২০০টি বিন্দু
stages = [C, Vt @ C, np.diag(s) @ Vt @ C, U @ np.diag(s) @ Vt @ C]
names = ["unit circle", "after Vt (rotate)", "after Σ (stretch)", "after U = A·x"]
fig, axes = plt.subplots(1, 4, figsize=(9, 2.8), layout="constrained")
for ax, P, name in zip(axes, stages, names):
    ax.plot(P[0], P[1], "C0")
    ax.plot(P[0, 0], P[1, 0], "o", color="C3")     # x = (1, 0) বিন্দুটি অনুসরণ
    ax.set_xlim(-7, 7); ax.set_ylim(-7, 7); ax.set_aspect("equal")
    ax.axhline(0, color="0.8", lw=0.8); ax.axvline(0, color="0.8", lw=0.8)
    ax.set_title(name, fontsize=9)
plt.show()

singular values: [6.7082 2.2361]
U =
 [[-0.3162 -0.9487]
 [-0.9487  0.3162]]
Vt =
 [[-0.7071 -0.7071]
 [-0.7071  0.7071]]
A == U Σ Vt ? True
eig(AᵀA): [ 5. 45.] → √ = [6.7082 2.2361]
V-র কলাম == AᵀA-র eigenvector (চিহ্ন বাদে)? True
u1 == A v1 / σ1 ? True
det(A) = 15.0 | σ1·σ2 = 15.0

৩ · Singular value, rank আর matrix-এর মাপ

A = UΣVᵀ গুণ করে খুলে লিখলে SVD একটি যোগফল হয়ে যায়:

A = σ₁u₁v₁ᵀ + σ₂u₂v₂ᵀ + … + σ_ru_rv_rᵀ

প্রতিটি u_iv_iᵀ একটি outer product: m × n matrix যার প্রতিটি সারি v_iᵀ-এর গুণিতক, তাই এর rank ঠিক ১. অর্থাৎ যেকোনো matrix হলো rank-১ “স্তর”-এর ওজন দেওয়া যোগফল, আর ওজনগুলো বড় থেকে ছোট সাজানো। r হলো শূন্য নয় এমন σ-এর সংখ্যা। এখান থেকে কয়েকটি গুরুত্বপূর্ণ সত্য সরাসরি আসে:

• Rank = শূন্য নয় এমন singular value-এর সংখ্যা। Column space-এর orthonormal basis হলো u₁…u_r, row space-এর v₁…v_r, আর null space-এর v_r+1…v_n. পর্ব ০৬-এর চারটি মৌলিক subspace এক SVD-তেই পাওয়া যায়।
• Numerical rank: float-এর হিসাবে “ঠিক শূন্য” প্রায় আসে না। np.linalg.matrix_rank একটি tolerance ধরে, ডিফল্টে σ_max × max(m, n) × machine epsilon (প্রায় 2.2 × 10⁻¹⁶), আর তার নিচের σ-কে শূন্য গোনে। ডেটায় মাপার ভুল থাকলে তার মাত্রা অনুযায়ী tolerance বড় করাই যুক্তিসঙ্গত।
• Spectral norm ‖A‖₂ = σ₁: একক vector-কে A সর্বোচ্চ যতগুণ লম্বা করতে পারে, অর্থাৎ ellipse-এর লম্বা অর্ধ-অক্ষ।
• Frobenius norm: ‖A‖_F² = সব ঘরের বর্গের যোগফল = σ₁² + … + σ_r². কারণ orthogonal matrix দিয়ে গুণ করলে Frobenius norm বদলায় না, তাই ‖A‖_F = ‖Σ‖_F. σ_i²-কে তাই অনেকে i-তম স্তরের “শক্তি” (energy) বলেন।
• বর্গাকার A-তে |det A| = σ₁σ₂…σ_n, আর A উল্টানো যায় ঠিক তখনই যখন σ_n > 0.

৪ · Eckart–Young — সবচেয়ে ভালো low-rank approximation

যোগফলের রূপ দেখে একটি স্বাভাবিক ধারণা আসে: শুধু প্রথম k-টি স্তর রাখলে কেমন হয়?

A_k = σ₁u₁v₁ᵀ + … + σ_ku_kv_kᵀ = U_kΣ_kV_kᵀ

১৯৩৬ সালে Carl Eckart আর Gale Young (Psychometrika-তে) দেখান, এটি কেবল একটি সুবিধাজনক পছন্দ নয়, সবচেয়ে ভালো পছন্দ। Eckart–Young উপপাদ্য: rank সর্বোচ্চ k এমন সব matrix B-এর মধ্যে ‖A − B‖ সবচেয়ে কম হয় B = A_k-এ, Frobenius আর spectral দুই norm-এই, আর ভুলের মান সঠিকভাবে জানা:

‖A − A_k‖_F² = σ_k+1² + … + σ_r²,   ‖A − A_k‖₂ = σ_k+1

প্রথম সমতা সহজে দেখা যায়: A − A_k বাদ পড়া স্তরগুলোর যোগফল, যার SVD-তেই σ_k+1, …, σ_r, তাই Frobenius² তাদের বর্গের যোগফল। “এর চেয়ে ভালো কিছু নেই” অংশটি গভীরতর। Spectral norm-এর প্রমাণের নকশা: যেকোনো rank-k B-এর null space-এর মাত্রা অন্তত n − k, আর v₁, …, v_k+1-এর span-এর মাত্রা k + ১; দুটি মিলে n-এর বেশি, তাই তারা একটি একক vector w-তে ছেদ করে। সেখানে Bw = 0, তাই ‖(A − B)w‖ = ‖Aw‖ ≥ σ_k+1. কোনো B-ই σ_k+1-এর নিচে নামতে পারে না, আর A_k ঠিক সেখানে পৌঁছায়।

Low-rank approximation-এর অর্থ ML-এ সবচেয়ে স্পষ্ট দেখা যায় রেটিং-এ। বাজারঘরের পাঁচজন কাল্পনিক ক্রেতা চারটি পণ্যকে ১ থেকে ৫ রেটিং দিয়েছেন: P1 চাল, P2 ডাল, P3 হেডফোন, P4 চার্জার। নিচের stepper প্রথম k-টি σ রেখে matrix আবার গড়ে, আর ভুল মেলায় সূত্রের সঙ্গে। সংখ্যাগুলো seed ছাড়াই নির্দিষ্ট (হাতে বানানো matrix), NumPy দিয়ে হিসাব করা।

৫ · ছবিতে low-rank — rank বনাম ভুল, সূত্রের সঙ্গে মিলিয়ে

AI-এর ভিত্তি কোর্সের পর্ব ১৮-এ দেখেছেন, ধূসর ছবি মানে একটি 2-D array. তাই ছবিতে Eckart–Young চোখে দেখা যায়। scikit-learn তার package-এর ভেতরে দুটি নমুনা ছবি রাখে (load_sample_image; ইন্টারনেট লাগে না); এখানে “china.jpg”, একটি মন্দিরের ছবি, যা Creative Commons BY 2.0 লাইসেন্সে Flickr ব্যবহারকারী danielbuechele-এর। ব্লকটি ছবিকে ধূসর করে, দ্রুততার জন্য প্রতি দ্বিতীয় pixel নেয় (২১৪ × ৩২০), তারপর k = ৫, ২০, ৫০ রেখে ভুল মাপে আর সূত্রের সঙ্গে মেলায়।

import numpy as np
import matplotlib.pyplot as plt
import PIL                                   # load_sample_image ছবি খুলতে Pillow লাগে
from sklearn.datasets import load_sample_image

img = load_sample_image("china.jpg").mean(axis=2)[::2, ::2]   # ধূসর, প্রতি দ্বিতীয় pixel
m, n = img.shape
U, s, Vt = np.linalg.svd(img, full_matrices=False)
print("shape:", img.shape, "| rank:", np.linalg.matrix_rank(img))
print("σ1..σ5:", np.round(s[:5], 1))

total = np.sum(s**2)                         # = ‖img‖_F²
print(" k  ‖A−A_k‖²_F   বাদ পড়া σ²   আপেক্ষিক ভুল%  জায়গা%")
fig, axes = plt.subplots(1, 4, figsize=(9, 2.2), layout="constrained")
axes[0].imshow(img, cmap="gray"); axes[0].set_title(f"original (rank {min(m, n)})", fontsize=9)
for ax, k in zip(axes[1:], (5, 20, 50)):
    Ak = (U[:, :k] * s[:k]) @ Vt[:k]
    err = np.sum((img - Ak) ** 2)
    tail = np.sum(s[k:] ** 2)
    store = k * (m + n + 1) / (m * n)
    print(f"{k:2d}  {err:.4e}  {tail:.4e}  {100*np.sqrt(tail/total):6.2f}  {100*store:6.1f}")
    ax.imshow(Ak, cmap="gray"); ax.set_title(f"rank {k}", fontsize=9)
for ax in axes: ax.axis("off")
plt.show()

print("spectral norm ‖A−A_20‖₂ =", round(np.linalg.norm(img - (U[:, :20] * s[:20]) @ Vt[:20], 2), 2),
      "| σ21 =", round(s[20], 2))

shape: (214, 320) | rank: 214
σ1..σ5: [41777.5  7683.2  4879.   2897.7  2454.2]
 k  ‖A−A_k‖²_F   বাদ পড়া σ²   আপেক্ষিক ভুল%  জায়গা%
 5  6.4183e+07  6.4183e+07   18.35     3.9
20  3.3498e+07  3.3498e+07   13.25    15.6
50  1.6092e+07  1.6092e+07    9.19    39.1
spectral norm ‖A−A_20‖₂ = 976.73 | σ21 = 976.73

৬ · Pseudo-inverse — least squares আর minimum-norm সমাধান এক সূত্রে

Ax = b সমাধান করতে A⁻¹ লাগে, কিন্তু A বর্গাকার আর পূর্ণ rank না হলে inverse নেই। দুই রকম সমস্যা আসে। লম্বা A (সমীকরণ বেশি, অজানা কম, যেমন linear regression): সাধারণত কোনো x সব সমীকরণ মানে না, তাই চাই এমন x যাতে ‖Ax − b‖ সবচেয়ে কম (least squares). চওড়া A (অজানা বেশি, underdetermined system): অসীম সমাধান, কোনটি নেব? SVD দুটিরই একটি স্বাভাবিক উত্তর দেয়। Σ-এর শূন্য নয় এমন প্রতিটি σ_i-কে 1/σ_i দিয়ে বদলে, শূন্যগুলো শূন্যই রেখে, আকার উল্টে (m × n থেকে n × m) পাওয়া যায় Σ⁺. তখন

A⁺ = V Σ⁺ Uᵀ,   x⁺ = A⁺b

এটিই Moore–Penrose pseudo-inverse. E. H. Moore ১৯২০ সালের দিকে আর Roger Penrose ১৯৫৫ সালে আলাদাভাবে এর ধারণা দেন; Penrose চারটি শর্ত দিয়ে এটিকে একমাত্র matrix হিসেবে চিহ্নিত করেন: AA⁺A = A, A⁺AA⁺ = A⁺, আর AA⁺ ও A⁺A দুটিই প্রতিসম। x⁺-এর দুটি গুণ একসঙ্গে প্রমাণ করা যায়: (১) সব x-এর মধ্যে x⁺ ‖Ax − b‖ সবচেয়ে কম করে; (২) এমন সব x-এর মধ্যে x⁺-এর নিজের দৈর্ঘ্য ‖x‖ সবচেয়ে কম। দ্বিতীয়টির কারণ: x⁺ পুরোপুরি row space-এ থাকে (v₁…v_r-এর span), আর যেকোনো অন্য সমাধান = x⁺ + (null space-এর একটি vector), যা x⁺-এর সঙ্গে লম্ব, তাই পিথাগোরাস অনুযায়ী শুধু দৈর্ঘ্য বাড়ায়। এই সমাধানকে বলে minimum-norm solution.

A লম্বা আর পূর্ণ কলাম-rank হলে A⁺ = (AᵀA)⁻¹Aᵀ, যা normal equations-এর সূত্র (AI-এর ভিত্তি কোর্সের পর্ব ৫৮; projection হিসেবে ব্যাখ্যা পর্ব ০৭-এ)। চওড়া আর পূর্ণ সারি-rank হলে A⁺ = Aᵀ(AAᵀ)⁻¹. SVD-এর সূত্র দুটিকেই ধরে, আর rank কম হলেও কাজ করে। সবচেয়ে ছোট উদাহরণে হাতে দেখি: একটি কলামের matrix a = (3, 4)ᵀ.

৭ · Condition number = σ_max/σ_min — কেন কিছু সমাধান কাঁপে

A⁺-এর সূত্রে 1/σ_i আছে। কোনো σ_i খুব ছোট হলে 1/σ_i বিশাল, আর b-এর সামান্য পরিবর্তন x⁺-এ বিশাল পরিবর্তন আনে। x⁺ = Σ (u_iᵀb/σ_i) v_i লিখলে দেখা যায়: b-কে δb পরিমাণ নাড়ালে u_i দিকের অংশ 1/σ_i গুণ বড় হয়ে x-এ ঢোকে। সবচেয়ে খারাপ দিক u_min. তাই আপেক্ষিক পরিবর্তনের সর্বোচ্চ গুণক হলো

κ(A) = σ_max / σ_min,   ‖δx‖/‖x‖ ≤ κ(A) · ‖δb‖/‖b‖  (বর্গাকার, উল্টানো যায় এমন A)

এই condition number AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ np.linalg.cond হিসেবে এসেছিল; এখন দেখছি এর উৎস। ellipse-এর ভাষায় κ হলো লম্বা আর খাটো অক্ষের অনুপাত: বৃত্তের কাছাকাছি হলে ভালো, সুচের মতো সরু হলে ill-conditioned. দুটি পরিণাম মনে রাখুন। এক, κ(AᵀA) = κ(A)², কারণ AᵀA-র singular value σ_i²; তাই normal equations সরাসরি সমাধান করলে হারানো অঙ্ক দ্বিগুণ হয়, আর lstsq SVD বা QR দিয়ে কাজ করে। দুই, ML-এ ছোট σ সাধারণত আসে প্রায়-একই feature থেকে (collinearity), আর তখন coefficient অস্থির হয়, যদিও অনুমান (Xw) প্রায় ঠিক থাকে।

প্রতিকার দুটি, দুটিই SVD-এর ভাষায় সরল। Truncated pseudo-inverse: σ_i একটি সীমার নিচে হলে 1/σ_i-এর বদলে ০ বসানো (np.linalg.pinv-এর rtol, পুরোনো নাম rcond). Ridge: 1/σ_i-এর বদলে σ_i/(σ_i² + λ), যা বড় σ-তে প্রায় 1/σ_i আর ছোট σ-তে প্রায় ০; একটি মসৃণ ছাঁটাই। নিচের ব্লক তিনটি দাবি যাচাই করে: SVD থেকে বানানো pseudo-inverse np.linalg.pinv-এর সমান, এক সমীকরণ তিন অজানার সমস্যায় এটি minimum-norm সমাধান দেয়, আর প্রায় একই দুটি feature-এ (ডেলিভারির সরলরেখা দূরত্ব আর রাস্তার দূরত্ব, যা প্রায় ১.৩ গুণ) condition number আর ছাঁটাইয়ের প্রভাব।

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

def pinv_svd(A, rtol=1e-12):                 # A⁺ = V Σ⁺ Uᵀ, ছোট σ বাদ
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    inv = np.where(s > rtol * s[0], 1 / s, 0.0)
    return Vt.T @ (inv[:, None] * U.T)

# ১) এক সমীকরণ, তিন অজানা: x1 + 2·x2 + 2·x3 = 9
A = np.array([[1., 2., 2.]]); b = np.array([9.])
x = pinv_svd(A) @ b
print("pinv-এর সমাধান:", x, "| ‖x‖ =", np.linalg.norm(x))
print("np.linalg.pinv-এর সঙ্গে মেলে?", np.allclose(pinv_svd(A), np.linalg.pinv(A)))
for other in ([9., 0., 0.], [3., 1., 2.]):
    print("আরেকটি সমাধান", other, "→ A·x =", A @ other, "‖x‖ =", round(np.linalg.norm(other), 4))

# ২) প্রায় একই দুটি feature: দূরত্ব (km) আর রাস্তার দূরত্ব ≈ 1.3 × km
rng = np.random.default_rng(42)
km = rng.uniform(1, 15, 30)
X = np.column_stack([km, 1.3 * km + rng.normal(0, 0.01, 30)])
y = 10 + 3 * km + rng.normal(0, 1, 30)            # ডেলিভারির মিনিট
Xc, yc = X - X.mean(0), y - y.mean()
s = np.linalg.svd(Xc, compute_uv=False)
print("σ:", s, "| σmax/σmin =", round(s[0] / s[-1], 1),
      "| np.linalg.cond =", round(np.linalg.cond(Xc), 1))
y2 = yc + rng.normal(0, 0.2, 30)                  # label-এ সামান্য নড়চড়
for rt in (1e-12, 1e-2):
    w1, w2 = pinv_svd(Xc, rt) @ yc, pinv_svd(Xc, rt) @ y2
    print(f"rtol={rt:g}: w = {w1} → নড়চড়ের পর {w2}")

pinv-এর সমাধান: [1. 2. 2.] | ‖x‖ = 3.0
np.linalg.pinv-এর সঙ্গে মেলে? True
আরেকটি সমাধান [9.0, 0.0, 0.0] → A·x = [9.] ‖x‖ = 9.0
আরেকটি সমাধান [3.0, 1.0, 2.0] → A·x = [9.] ‖x‖ = 3.7417
σ: [36.5872  0.0264] | σmax/σmin = 1386.7 | np.linalg.cond = 1386.7
rtol=1e-12: w = [ 22.4617 -15.0066] → নড়চড়ের পর [ 21.7716 -14.4774]
rtol=0.01: w = [1.0985 1.428 ] → নড়চড়ের পর [1.0977 1.4269]

৮ · Sparse লেখায় truncated SVD — TF-IDF থেকে LSA

AI-এর ভিত্তি কোর্সের পর্ব ৬৭-এ TF-IDF matrix এসেছিল আর পর্ব ৬৯-এ তার SVD নিয়ে LSA (Deerwester ও সহকর্মীরা, ১৯৯০)। সেখানে কাজ হয়েছিল ঘন (dense) array-তে np.linalg.svd দিয়ে। বাস্তব সংগ্রহে লক্ষ লক্ষ নথি আর লক্ষ শব্দ; matrix-এর প্রায় সব ঘর শূন্য, তাই সেটি sparse matrix হিসেবে রাখা হয়। ঘন করলে memory ফেটে যায়, আর পুরো SVD-এর দরকারও নেই, প্রথম কয়েকশো σ-ই যথেষ্ট। scikit-learn-এর TruncatedSVD ঠিক এজন্য: sparse ইনপুট নেয়, ডিফল্টে randomized SVD চালায় (এলোমেলো vector দিয়ে column space-এর একটি ছোট basis খুঁজে তার ওপর ছোট একটি SVD; n_iter=5, n_oversamples=10), আর শুধু k-টি দিক দেয়।

একটি গুরুত্বপূর্ণ পার্থক্য: TruncatedSVD কলামের গড় বাদ দেয় না, PCA দেয়। গড় বাদ দিলে প্রতিটি শূন্য ঘর একটি ঋণাত্মক সংখ্যা হয়ে যায়, sparse আর থাকে না। তাই TruncatedSVD-এর প্রথম দিক প্রায়ই শুধু “গড় নথি”-র দিক, আর এর explained variance-এর ক্রম σ-এর ক্রমের সঙ্গে না-ও মিলতে পারে। নিচের ব্লক বাজারঘরের দশটি কাল্পনিক বাংলা রিভিউয়ে এটি দেখায়।

import numpy as np
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD

reviews = [                                   # কাল্পনিক বাজারঘর রিভিউ
    "ডেলিভারি খুব দ্রুত এসেছে", "ডেলিভারি দেরিতে এসেছে রাইডার ফোন ধরেনি",
    "রাইডার ভদ্র ডেলিভারি দ্রুত", "দুই দিন দেরিতে ডেলিভারি পেলাম",
    "শাড়ির কাপড় ভালো রং সুন্দর", "কাপড় পাতলা রং উঠে যায়",
    "রং ছবির মতো কাপড় নরম", "দাম কম কিন্তু কাপড় ভালো",
    "দাম বেশি তবে ডেলিভারি দ্রুত", "রাইডার দেরিতে এসেছে কাপড় ভালো",
]
vec = TfidfVectorizer(token_pattern=r"\S+")   # ফাঁকা জায়গায় ভাঙা: বাংলা যুক্তবর্ণ অক্ষত
X = vec.fit_transform(reviews)
print("X:", X.shape, "| শূন্য নয় এমন ঘর:", X.nnz, f"({100 * X.nnz / np.prod(X.shape):.0f}%)")

svd = TruncatedSVD(n_components=2, random_state=42)
Z = svd.fit_transform(X)                      # প্রতিটি রিভিউ এখন ২টি সংখ্যা
print("σ (TruncatedSVD):", np.round(svd.singular_values_, 4))
print("σ (পুরো SVD)    :", np.round(np.linalg.svd(X.toarray(), compute_uv=False)[:2], 4))
print("explained variance ratio:", np.round(svd.explained_variance_ratio_, 3))

terms = vec.get_feature_names_out()
for i, comp in enumerate(svd.components_):
    top = np.argsort(comp)[::-1][:4]
    print(f"দিক {i + 1}:", ", ".join(f"{terms[j]} {comp[j]:+.2f}" for j in top))
for r, z in zip(reviews, Z):
    print(f"{z[0]:+.2f} {z[1]:+.2f}  {r}")

X: (10, 28) | শূন্য নয় এমন ঘর: 49 (18%)
σ (TruncatedSVD): [1.5    1.2907]
σ (পুরো SVD)    : [1.5    1.2907]
explained variance ratio: [0.043 0.196]
দিক 1: ডেলিভারি +0.39, রাইডার +0.38, এসেছে +0.38, দ্রুত +0.32
দিক 2: কাপড় +0.46, রং +0.40, ভালো +0.33, সুন্দর +0.22
+0.59 -0.36  ডেলিভারি খুব দ্রুত এসেছে
+0.67 -0.22  ডেলিভারি দেরিতে এসেছে রাইডার ফোন ধরেনি
+0.59 -0.36  রাইডার ভদ্র ডেলিভারি দ্রুত
+0.36 -0.17  দুই দিন দেরিতে ডেলিভারি পেলাম
+0.32 +0.67  শাড়ির কাপড় ভালো রং সুন্দর
+0.19 +0.54  কাপড় পাতলা রং উঠে যায়
+0.19 +0.54  রং ছবির মতো কাপড় নরম
+0.33 +0.44  দাম কম কিন্তু কাপড় ভালো
+0.42 -0.25  দাম বেশি তবে ডেলিভারি দ্রুত
+0.72 +0.23  রাইডার দেরিতে এসেছে কাপড় ভালো

৯ · Recommender-এর ভিত্তি — matrix factorization

অংশ ৪-এর stepper আসলে একটি ছোট recommender system-এর ভিত্তি দেখিয়েছে। R ≈ U_kΣ_kV_kᵀ লিখে P = U_kΣ_k^1/2 আর Q = V_kΣ_k^1/2 নিলে R ≈ PQᵀ: প্রতিটি ক্রেতা একটি k-মাত্রার vector p_u, প্রতিটি পণ্য একটি q_i, আর আনুমানিক রেটিং তাদের dot product p_uᵀq_i. এই কাঠামোর নাম matrix factorization, আর এটি collaborative filtering-এর একটি প্রধান পদ্ধতি: অন্য ক্রেতাদের আচরণ থেকে একজনের পছন্দ আন্দাজ।

কিন্তু বাস্তবে একটি বড় বাধা: রেটিং matrix-এর বেশিরভাগ ঘর খালি। একজন ক্রেতা হাজার পণ্যের মধ্যে হয়তো দশটিতে রেটিং দেন। খালি ঘরকে ০ ধরলে SVD শেখে “এই ক্রেতা এই পণ্য অপছন্দ করেন”, যা ভুল; আর খালি ঘরের SVD সংজ্ঞায়িতই নয়। তাই recommender-এর “SVD” আসলে SVD নয়, শুধু যে ঘরগুলো জানা সেখানে ভুল কমানো একটি optimization:

min_{P, Q} ∑_{(u,i) জানা} (r_ui − p_uᵀq_i)² + λ(‖p_u‖² + ‖q_i‖²)

এখানে r_ui ক্রেতা u-এর পণ্য i-কে দেওয়া রেটিং, λ regularization-এর শক্তি। এর জনপ্রিয়তা আসে Netflix Prize (২০০৬–২০০৯) থেকে: Simon Funk ছদ্মনামে একজন প্রতিযোগী ২০০৬ সালে ব্লগে stochastic gradient descent দিয়ে এভাবে factor শেখার পদ্ধতি লেখেন, যা “Funk SVD” নামে পরিচিত হয়; পরে Koren, Bell আর Volinsky ২০০৯ সালে IEEE Computer-এ এর bias-সহ রূপ গুছিয়ে লেখেন। বিকল্প পথ ALS (alternating least squares): Q স্থির রেখে P-এর জন্য least squares, তারপর উল্টো, পালাক্রমে; প্রতিটি ধাপ এই পর্বের pseudo-inverse-এর মতো ছোট সমস্যা। নতুন ক্রেতা বা নতুন পণ্যের কোনো রেটিং না থাকলে তার vector শেখার কিছু নেই; এই cold start সমস্যা, implicit feedback (ক্লিক, কেনা) আর মূল্যায়নের পদ্ধতি পর্ব ৫৩ (Recommender system — collaborative filtering আর matrix factorization)-এ আসবে।

১০ · AI-এর যুগে — SVD-র হিসাব আর AI-এর কোড যাচাই

AI assistant linear algebra-র ধাপগুলো সুন্দর করে লেখে, আর ছোট চিহ্নের ভুল চোখে পড়ে না। SVD-তে কয়েকটি ভুল বিশেষভাবে বারবার দেখা যায়।

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

• UΣVᵀ গুণ করে A ফেরে কি না, আর UᵀU = I, VᵀV = I কি না দেখেছি
• σ বড় থেকে ছোট, অঋণাত্মক, আর σ₁σ₂… = |det A| (বর্গাকার হলে) মিলিয়েছি
• ছোট একটি 2×2 উদাহরণে নিজে হাতে করে AI-এর সংখ্যার সঙ্গে মিলিয়েছি
• Low-rank ভুল বাদ পড়া σ²-এর যোগফলের সমান কি না পরীক্ষা করেছি
• Recommender বা মডেলের ফল তালা দেওয়া test ডেটায়, একটি সরল baseline-এর পাশে, কয়েকটি seed-এ মেপেছি
• Library-র parameter নাম (rtol, full_matrices, n_iter) নিজের version-এর অফিসিয়াল docs-এ মিলিয়েছি
• আসল গ্রাহকের রেটিং বা রিভিউ AI-তে পেস্ট করিনি