Matrix-এর জ্যামিতি — projection, rotation আর least squares

AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ NumPy-তে matrix গুণ, solve, lstsq আর condition number শিখেছেন, পর্ব ৩৩-এ norm, dot product আর কোণ, আর পর্ব ৫৮-এ normal equations আর “y-এর ছায়া” হিসেবে linear regression-এর প্রথম ছবি। পর্ব ০৬ (Vector space — span, basis, linear independence আর rank) দিয়েছে column space আর rank-এর ভাষা। এই পর্বের প্রশ্ন: একটি matrix জ্যামিতিকভাবে কী করে, আর least squares সমাধানকারী software ভেতরে কেন QR ব্যবহার করে? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy আর matplotlib); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে।

১ · Matrix মানে একটি map — বর্গাকার জাল কীভাবে বাঁকে

একটি function f যদি দুটি নিয়ম মানে, f(u + v) = f(u) + f(v) আর f(cu) = c·f(u) (যেকোনো vector u, v আর সংখ্যা c-এর জন্য), তাকে বলে linear map বা linear transformation. নিয়ম দুটির সরাসরি ফল: origin নড়ে না (f(0) = 0), সোজা রেখা সোজা থাকে, আর সমান ফাঁকের বিন্দু সমান ফাঁকেই থাকে। প্রতিটি linear map-কে একটি matrix দিয়ে লেখা যায়, আর প্রতিটি matrix একটি linear map.

কেন? দুই মাত্রায় যেকোনো vector x = (x₁, x₂) লেখা যায় x₁e₁ + x₂e₂, যেখানে e₁ = (1, 0) আর e₂ = (0, 1) হলো standard basis. Linearity বলে f(x) = x₁f(e₁) + x₂f(e₂). তাই f(e₁) আর f(e₂)-কে পাশাপাশি কলাম হিসেবে বসালেই matrix A পাওয়া যায়, আর Ax ঠিক এই weighted যোগফল:

Ax = x₁ · (A-এর ১ম কলাম) + x₂ · (A-এর ২য় কলাম)

AI-এর ভিত্তি কোর্সের পর্ব ১৭-এ Ax পড়েছিলেন “সারি ধরে ধরে”, প্রতিটি সারির সঙ্গে x-এর dot product. এটি “কলাম ধরে” পড়া: একই হিসাব, কিন্তু জ্যামিতির জন্য এটিই বেশি কাজের। নিচের ছবিতে A = [[2, 1], [1, 2]] (সারি ধরে লেখা; NumPy-র মতো)। প্রথম কলাম (2, 1): e₁ সেখানে যায়; দ্বিতীয় কলাম (1, 2): e₂ সেখানে। একক বর্গটি একটি parallelogram হয়ে যায়, আর বর্গের কোণ (1, 1) যায় (2, 1) + (1, 2) = (3, 3)-এ।

ML-এ এই ছবি সব জায়গায়। একটি design matrix X আর weight w-এর গুণ Xw হলো X-এর কলামগুলোর (feature-গুলোর) weighted যোগফল; neural network-এর একটি layer Wx + b হলো একটি linear map আর তারপর একটি সরণ, যাকে বলে affine map. সরণ b origin সরায় বলে affine map নিজে linear নয়, কিন্তু তার ভেতরের জ্যামিতি পুরোটাই W-এর।

২ · পাঁচটি মৌলিক map — rotation, scaling, shear, reflection, projection

যেকোনো ২×২ matrix কয়েকটি সরল map জোড়া দিয়ে তৈরি হতে পারে (পর্ব ০৯ (SVD — low-rank approximation, pseudo-inverse আর recommender-এর ভিত্তি) দেখাবে যে ঠিক তিনটিই যথেষ্ট: ঘোরানো, টানা, আবার ঘোরানো)। তাই মৌলিক পাঁচটি চিনে রাখা দরকার। প্রতিটির কলাম দেখে বলুন e₁ আর e₂ কোথায় গেল।

Rotation-এর matrix মুখস্থ না করে কলাম থেকে বের করুন: e₁ = (1, 0)-কে θ কোণে ঘোরালে একক বৃত্তের ওপর (cos θ, sin θ); e₂ = (0, 1)-কে ঘোরালে (−sin θ, cos θ). θ = 30° হলে cos θ ≈ 0.866, sin θ = 0.5. নিচের ব্লক একক বর্গ আর ২৫টি বিন্দুর একটি জালকে পাঁচটি map-এ পাঠায়, প্রতিটির ছবি আঁকে, আর determinant-এর পাশে ছবির আসল ক্ষেত্রফল (বহুভুজের shoelace সূত্রে মাপা) ছাপে। সবুজ তীর প্রতিটি ছবিতে e₁-এর গন্তব্য, অর্থাৎ matrix-এর প্রথম কলাম।

import numpy as np
import matplotlib.pyplot as plt

t = np.pi / 6                                       # 30 ডিগ্রি
maps = {
    "rotate 30°": np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]]),
    "scale (2, 0.75)": np.array([[2.0, 0.0], [0.0, 0.75]]),
    "shear k = 1": np.array([[1.0, 1.0], [0.0, 1.0]]),
    "reflect in y = x": np.array([[0.0, 1.0], [1.0, 0.0]]),
    "project onto x-axis": np.array([[1.0, 0.0], [0.0, 0.0]]),
}
sq = np.array([[0, 1, 1, 0], [0, 0, 1, 1]])          # একক বর্গের ৪টি কোণ, কলামে
g = np.linspace(0, 1, 5)
grid = np.array([[a, b] for a in g for b in g]).T    # 2 x 25 বিন্দুর জাল

def area(P):                                        # shoelace সূত্রে বহুভুজের ক্ষেত্রফল
    x, y = P
    return 0.5 * abs(x @ np.roll(y, -1) - y @ np.roll(x, -1))

fig, axes = plt.subplots(2, 3, figsize=(7.2, 5.2), layout="constrained")
for ax, (name, A) in zip(axes.flat, maps.items()):
    img, d = A @ sq, np.linalg.det(A)
    ax.fill(*sq, color="0.85", label="unit square")
    ax.fill(*img, color="C0", alpha=0.45, label="image")
    ax.scatter(*(A @ grid), s=7, color="C3")
    ax.arrow(0, 0, *A[:, 0], color="C2", width=0.02, length_includes_head=True)
    ax.set_title(f"{name}\ndet = {d:.2f}", fontsize=10)
    ax.set_aspect("equal"); ax.set_xlim(-0.8, 2.2); ax.set_ylim(-0.6, 1.6)
    side = "উল্টে যায়" if d < 0 else ("চ্যাপ্টা" if np.isclose(d, 0) else "ঠিক থাকে")
    print(f"{name:20s} det = {d:6.3f}  ক্ষেত্রফল = {area(img):.3f}  দিক: {side}")
axes.flat[5].axis("off")
axes.flat[0].legend(fontsize=8, loc="upper left")
plt.show()

rotate 30°           det =  1.000  ক্ষেত্রফল = 1.000  দিক: ঠিক থাকে
scale (2, 0.75)      det =  1.500  ক্ষেত্রফল = 1.500  দিক: ঠিক থাকে
shear k = 1          det =  1.000  ক্ষেত্রফল = 1.000  দিক: ঠিক থাকে
reflect in y = x     det = -1.000  ক্ষেত্রফল = 1.000  দিক: উল্টে যায়
project onto x-axis  det =  0.000  ক্ষেত্রফল = 0.000  দিক: চ্যাপ্টা

৩ · পরপর দুটি map মানে matrix গুণ — আর ক্রম কেন জরুরি

প্রথমে R প্রয়োগ, তারপর S: ফল S(Rx). Matrix গুণের সংজ্ঞাই এমনভাবে বানানো যাতে (SR)x = S(Rx) হয়; অর্থাৎ matrix গুণ মানে map-এর composition. লেখার ক্রম আর প্রয়োগের ক্রম উল্টো: SR-এ ডানেরটি আগে। কলামের ভাষায়: SR-এর প্রথম কলাম হলো S(R-এর প্রথম কলাম), অর্থাৎ e₁ আগে R-এ, তারপর S-এ।

শেষ ধাপের নিয়মটি রোজকার: আগে মোজা, তারপর জুতা পরলে খোলার সময় আগে জুতা, তারপর মোজা। এবার ML-এর একটি গুরুত্বপূর্ণ ফল: দুটি linear layer পরপর বসালে W₂(W₁x) = (W₂W₁)x, অর্থাৎ একটিই linear map. একশটি linear layer-ও একটি matrix-এর বেশি কিছু বলতে পারে না। এজন্যই neural network-এ প্রতিটি layer-এর পর ReLU-এর মতো একটি non-linear function বসে; AI-এর ভিত্তি কোর্সের পর্ব ৬৮-এ XOR-এর উদাহরণে এটি দেখেছেন। আরেকটি ব্যবহারিক কথা: বাজারঘরের ডেলিভারি অ্যাপে মানচিত্রের দিক ঘোরানো, জুম আর হেলানো পরপর হলে প্রতিটি বিন্দুতে তিনবার গুণ না করে আগে তিনটি matrix গুণ করে একটি M বানিয়ে নেওয়া সস্তা।

৪ · Determinant — ক্ষেত্রফল কত গুণ, আর দিক উল্টায় কি না

২×২ matrix A = [[a, b], [c, d]]-এর determinant det A = ad − bc. জ্যামিতিক অর্থ: A-এর কলাম দুটি (a, c) আর (b, d) যে parallelogram গড়ে, তার চিহ্নসহ ক্ষেত্রফল। একক বর্গের ক্ষেত্রফল ১, তাই det A বলে A যেকোনো আকৃতির ক্ষেত্রফল কত গুণ করে, কারণ যেকোনো আকৃতিকে ছোট ছোট বর্গে ভাঙা যায় আর linear map প্রতিটি বর্গকে একই অনুপাতে বদলায়। তিন মাত্রায় det হলো একক ঘনকের আয়তন কত গুণ হলো।

ad − bc কেন? একটি ছোট প্রমাণ। ধরুন a ≠ 0. দ্বিতীয় কলাম থেকে প্রথম কলামের (b/a) গুণ বাদ দিলে, অর্থাৎ একটি shear করলে, parallelogram-এর ভূমি একই থাকে আর উচ্চতাও (shear ক্ষেত্রফল বদলায় না, অংশ ২-এর আউটপুট)। নতুন দ্বিতীয় কলাম (0, d − bc/a). এখন parallelogram-এর একটি বাহু খাড়া, দৈর্ঘ্য d − bc/a, আর অন্য বাহুর অনুভূমিক বিস্তার a. ক্ষেত্রফল = a × (d − bc/a) = ad − bc. চিহ্নের অর্থ: প্রথম কলাম থেকে দ্বিতীয় কলামে যেতে ঘড়ির কাঁটার উল্টো দিকে ঘুরতে হলে ধনাত্মক, ঘড়ির কাঁটার দিকে হলে ঋণাত্মক।

তিনটি নিয়ম, প্রতিটির জ্যামিতিক কারণসহ:

• det(AB) = det A · det B. আগে B ক্ষেত্রফল det B গুণ করে, তারপর A আরও det A গুণ। অংশ ৩-এ det(SR) = ১ × ১ = ১.
• det A = 0 ⇔ A-এর inverse নেই। ক্ষেত্রফল শূন্য মানে plane একটি রেখায় (বা বিন্দুতে) চ্যাপ্টা হয়েছে; অনেক বিন্দু একই জায়গায় পড়ে, তাই কোনটি থেকে এসেছিল আর বলা যায় না। এমন matrix-কে বলে singular; এর কলামগুলো linearly dependent, আর rank মাত্রার চেয়ে কম (পর্ব ০৬-এর ভাষা)।
• det(A⁻¹) = 1 / det A. Inverse map-কে undo করে, তাই ক্ষেত্রফলের গুণকও উল্টায়।

২×২-এর inverse: A⁻¹ = (1 / det A) · [[d, −b], [−c, a]]. পাঁচটি map-এর inverse জ্যামিতি দিয়েই বলা যায়: rotation θ-এর inverse rotation −θ; scaling (s₁, s₂)-এর inverse (1/s₁, 1/s₂); shear k-এর inverse shear −k; reflection নিজেই নিজের inverse (দুবার আয়নায় উল্টালে আগের মতো); আর projection-এর inverse নেই, কারণ মোছা উচ্চতা ফেরানোর উপায় নেই।

একটি পূর্ণ numerical. A = [[3, 1], [2, 4]]. det A = 3·4 − 1·2 = 10: যেকোনো আকৃতির ক্ষেত্রফল ১০ গুণ, দিক অক্ষত। A⁻¹ = (1/10)·[[4, −1], [−2, 3]] = [[0.4, −0.1], [−0.2, 0.3]]. যাচাই: A·A⁻¹-এর প্রথম সারি (3·0.4 + 1·(−0.2), 3·(−0.1) + 1·0.3) = (1, 0). b = (5, 10) হলে Ax = b-এর সমাধান x = A⁻¹b = (0.4·5 − 0.1·10, −0.2·5 + 0.3·10) = (1, 2); যাচাই: 3·1 + 1·2 = 5, 2·1 + 4·2 = 10.

৫ · Orthogonal matrix — দৈর্ঘ্য আর কোণ অক্ষত

পাঁচটি map-এর মধ্যে rotation আর reflection বিশেষ: এরা আকৃতি বদলায় না, শুধু ঘোরায় বা উল্টায়। এমন matrix Q-এর কলামগুলো একে অপরের সঙ্গে লম্ব আর প্রতিটির দৈর্ঘ্য ১, অর্থাৎ orthonormal. কলামের dot product-গুলো একসঙ্গে লিখলে শর্তটি এক লাইনে আসে:

QᵀQ = I   ⇔   Q⁻¹ = Qᵀ

এখানে Qᵀ হলো Q-এর transpose (সারি আর কলাম অদলবদল), আর I হলো identity matrix. QᵀQ-এর (i, j) ঘর = i-তম কলাম · j-তম কলাম; i = j হলে ১ (দৈর্ঘ্য ১), না হলে ০ (লম্ব)। এমন বর্গ matrix-কে বলে orthogonal matrix. তিনটি ফল, প্রতিটি এক লাইনে প্রমাণ হয়:

• দৈর্ঘ্য অক্ষত: ‖Qx‖² = (Qx)ᵀ(Qx) = xᵀQᵀQx = xᵀx = ‖x‖².
• কোণ অক্ষত: (Qx)·(Qy) = xᵀQᵀQy = x·y; দৈর্ঘ্য আর dot product দুটিই না বদলালে cosine-ও বদলায় না।
• det = ±1: det(QᵀQ) = (det Q)² = det I = 1. +১ হলে rotation, −১ হলে তাতে একটি reflection মেশানো।

সংখ্যার হিসাবের দিক থেকে orthogonal matrix সবচেয়ে নিরাপদ: inverse বের করতে কোনো ভাগ লাগে না (শুধু transpose), আর condition number ঠিক ১, তাই float-এর ছোট ভুল একটুও বড় হয় না। এজন্যই ভালো algorithm চেষ্টা করে কাজের বড় অংশ orthogonal matrix দিয়ে করতে। একটি গুরুত্বপূর্ণ উদাহরণ Householder reflection: একক vector u-এর জন্য H = I − 2uuᵀ. এটি u-এর লম্ব সমতলে আয়না: u-এর দিকের অংশ উল্টে দেয়, বাকি সব অক্ষত। যাচাই করুন: Hᵀ = H, আর H² = I − 4uuᵀ + 4u(uᵀu)uᵀ = I. অংশ ৮-এ দেখবেন QR-এর software এই আয়নাগুলো দিয়েই চলে।

৬ · Orthogonal projection — একটি রেখায়, তারপর একটি subspace-এ

রেখায় projection. রেখাটি origin দিয়ে vector a-এর দিকে। খুঁজছি p = c·a, এমন c যাতে ভুলের অংশ r = b − ca রেখার সঙ্গে লম্ব: aᵀ(b − ca) = 0 ⇒ aᵀb − c·aᵀa = 0 ⇒ c = aᵀb / aᵀa. তাই p = a(aᵀb)/(aᵀa) = [aaᵀ / (aᵀa)] b, আর projection matrix P = aaᵀ / (aᵀa). aaᵀ একটি কলাম আর একটি সারির গুণ, একটি পূর্ণ matrix (একে বলে outer product).

Subspace-এ projection. এবার A একটি m×k matrix, যার k-টি কলাম linearly independent, আর তাদের column space col(A) হলো একটি k-মাত্রার subspace (তিন মাত্রায় k = 2 হলে একটি সমতল)। col(A)-এর প্রতিটি বিন্দু Ax̂ আকারের। একই যুক্তি: ভুলের অংশ r = b − Ax̂ যেন A-এর প্রতিটি কলামের সঙ্গে লম্ব হয়, অর্থাৎ Aᵀ(b − Ax̂) = 0:

AᵀA x̂ = Aᵀb  ⇒  x̂ = (AᵀA)⁻¹Aᵀb  ⇒  p = Ax̂ = A(AᵀA)⁻¹Aᵀ b

P = A(AᵀA)⁻¹Aᵀ

k = 1 হলে A = a, আর এটি ঠিক আগের aaᵀ / aᵀa. AᵀA-এর inverse থাকে ঠিক তখনই, যখন A-এর কলামগুলো independent (পর্ব ০৬: full column rank). এবার দুটি মূল ধর্মের প্রমাণ:

• P² = P (idempotent): P·P = A(AᵀA)⁻¹(AᵀA)(AᵀA)⁻¹Aᵀ = A(AᵀA)⁻¹Aᵀ = P. মাঝের AᵀA তার inverse-এর সঙ্গে কেটে যায়।
• Pᵀ = P (symmetric): (AᵀA) symmetric, তাই তার inverse-ও; Pᵀ = A((AᵀA)⁻¹)ᵀAᵀ = P. এই দ্বিতীয় ধর্মটিই projection-কে “orthogonal” করে; শুধু P² = P হলে তির্যক (oblique) projection-ও হতে পারত।
• সবচেয়ে কাছে: col(A)-এর যেকোনো q-এর জন্য b − q = r + (p − q), আর p − q col(A)-এ বলে r-এর সঙ্গে লম্ব। পিথাগোরাসে ‖b − q‖² = ‖r‖² + ‖p − q‖² ≥ ‖r‖², সমান শুধু q = p হলে।

আরও দুটি দরকারি তথ্য। (১) I − P-ও একটি projection, col(A)-এর লম্ব দিকে: (I − P)b = r. (২) Q-এর কলাম orthonormal হলে QᵀQ = I, তাই সূত্রটি ছোট হয়ে P = QQᵀ: কোনো inverse লাগে না। অংশ ৮-এর QR এই সুবিধাটিই কাজে লাগায়। নিচের ব্লক তিন মাত্রায় একটি সমতলে projection করে প্রতিটি দাবি সংখ্যায় যাচাই করে।

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

A = np.array([[1.0, 3.0],
              [2.0, 3.0],
              [2.0, 0.0]])         # দুটি কলাম: R³-এর ভেতর একটি সমতল
b = np.array([3.0, 0.0, 9.0])     # সমতলের বাইরের একটি vector

P = A @ np.linalg.inv(A.T @ A) @ A.T     # শেখার জন্য সূত্রটি হুবহু
p = P @ b                                # সমতলে b-এর ছায়া
r = b - p                                # ছায়া থেকে b পর্যন্ত লম্ব অংশ
print("P =\n", P)
print("p =", p, "  r =", r)
print("P @ P == P ?", np.allclose(P @ P, P))
print("P.T == P ?", np.allclose(P.T, P))
print("A.T @ r =", A.T @ r, "← দুটি কলামের সঙ্গেই লম্ব")
print("rank(P) =", np.linalg.matrix_rank(P), "  trace(P) =", round(np.trace(P), 6))

rng = np.random.default_rng(42)          # সমতলের আরও ৫টি এলোমেলো বিন্দু
others = A @ rng.normal(size=(2, 5))
dist = np.linalg.norm(b[:, None] - others, axis=0)
print("‖r‖ =", round(np.linalg.norm(r), 3), "  অন্য বিন্দুগুলোর দূরত্ব:", dist)

M = np.eye(3) - P                        # লম্ব দিকের projection
print("(I − P) b =", M @ b, "  (I − P)² == I − P ?", np.allclose(M @ M, M))
print("P সমতলের vector বদলায় না:", np.allclose(P @ A, A))

P =
 [[ 0.556  0.444 -0.222]
 [ 0.444  0.556  0.222]
 [-0.222  0.222  0.889]]
p = [-0.333  3.333  7.333]   r = [ 3.333 -3.333  1.667]
P @ P == P ? True
P.T == P ? True
A.T @ r = [0. 0.] ← দুটি কলামের সঙ্গেই লম্ব
rank(P) = 2   trace(P) = 2.0
‖r‖ = 5.0   অন্য বিন্দুগুলোর দূরত্ব: [11.174 11.79   8.171  7.647 16.267]
(I − P) b = [ 3.333 -3.333  1.667]   (I − P)² == I − P ? True
P সমতলের vector বদলায় না: True

৭ · Least squares মানে projection — y-এর ছায়া X-এর column space-এ

AI-এর ভিত্তি কোর্সের পর্ব ৫৮-এ MSE-এর gradient শূন্য বসিয়ে normal equations XᵀXw = Xᵀy পেয়েছিলেন, আর দেখেছিলেন inverse-এর পথে condition number বর্গ হয়ে যায়। এখানে একই সমীকরণ আসে কোনো derivative ছাড়াই, শুধু জ্যামিতি থেকে, আর তা থেকে অনেক কিছু একসঙ্গে বোঝা যায়।

দৃষ্টিভঙ্গি বদলান: n-টি অর্ডারের ডেলিভারি সময় y একটি n-মাত্রার একটি vector. X-এর প্রতিটি কলাম (সব-১-এর কলাম, দূরত্বের কলাম, ওজনের কলাম) একটি করে n-মাত্রার vector. যেকোনো w-এর জন্য অনুমান Xw এই কলামগুলোর weighted যোগফল, অর্থাৎ col(X)-এর একটি বিন্দু। Least squares চায় ‖y − Xw‖ সবচেয়ে ছোট, অর্থাৎ col(X)-এর যে বিন্দু y-এর সবচেয়ে কাছে। অংশ ৬ বলে সেটি y-এর orthogonal projection, আর শর্ত Xᵀ(y − Xŵ) = 0: ঠিক normal equations.

Projection matrix-টির এখানে একটি নিজস্ব নাম আছে, hat matrix: H = X(XᵀX)⁻¹Xᵀ, কারণ সে y-কে “টুপি পরায়”, ŷ = Hy. অংশ ৬-এর সব ধর্ম সরাসরি regression-এর তথ্যে বদলে যায়:

• Residual সব কলামের সঙ্গে লম্ব (Xᵀr = 0). X-এ সব-১-এর কলাম থাকলে তার সঙ্গে লম্ব মানে ∑r_i = 0: intercept-সহ model-এ residual-এর যোগফল সবসময় শূন্য।
• ŷ আর r লম্ব, তাই পিথাগোরাস: ‖y‖² = ‖ŷ‖² + ‖r‖². y থেকে গড় বাদ দিয়ে একই কথা লিখলে পরিচিত ভাঙন আসে: মোট বর্গ = ব্যাখ্যা করা বর্গ + residual বর্গ (SST = SSR + SSE), যা থেকে R² = 1 − SSE/SST.
• trace(H) = p, parameter-এর সংখ্যা (intercept-সহ)। r থাকে col(X)-এর লম্ব, n − p মাত্রার একটি subspace-এ; এজন্যই গোলমালের variance-এর নিরপেক্ষ আন্দাজে n দিয়ে নয়, n − p দিয়ে ভাগ করা হয়।
• H-এর কর্ণের ঘর h_ii-কে বলে i-তম বিন্দুর leverage: ŷ_i তার নিজের y_i-এর ওপর কতটা নির্ভর করে। ০ থেকে ১-এর মধ্যে থাকে, আর গড়ে p/n.

পূর্ণ numerical (পর্ব ৫৮-এর তিনটি বিন্দু, এবার projection-এর চোখে). দূরত্ব x = 2, 4, 6 km, সময় y = 20, 26, 38 মিনিট; X = [[1, 2], [1, 4], [1, 6]]. XᵀX = [[3, 12], [12, 56]], det = 168 − 144 = 24, (XᵀX)⁻¹ = (1/24)·[[56, −12], [−12, 3]]. H-এর (i, j) ঘর = (56 − 12x_i − 12x_j + 3x_ix_j) / 24. বসালে:

H = [[5/6, 1/3, −1/6], [1/3, 1/3, 1/3], [−1/6, 1/3, 5/6]]

যাচাই: Hy-এর প্রথম ঘর (5/6)·20 + (1/3)·26 − (1/6)·38 = 16.667 + 8.667 − 6.333 = 19; মাঝের ঘর (20 + 26 + 38)/3 = 28; শেষ ঘর 37. তাই ŷ = (19, 28, 37), r = (1, −2, 1), আর Xᵀr = (1 − 2 + 1, 2 − 8 + 6) = (0, 0). trace(H) = 5/6 + 1/3 + 5/6 = 2 = parameter-এর সংখ্যা (b আর w). Leverage দুই প্রান্তে ৫/৬, মাঝে ১/৩: প্রান্তের বিন্দু রেখাটিকে বেশি টানে, যেমন লম্বা লাঠির প্রান্ত ধরে ঘোরানো সহজ। মাঝের বিন্দুর অনুমান ঠিক তিনটি y-এর গড়, কারণ x = 4 হলো x-গুলোর গড়, আর least squares রেখা সবসময় (x̄, ȳ) দিয়ে যায়। উচ্চ leverage-এর outlier কীভাবে fit নষ্ট করে আর তা সামলানোর পথ পর্ব ২৪ (Robust আর quantile regression — outlier আর অনুমানের পরিসর)-এ আসবে।

৮ · Gram–Schmidt আর QR — হাতে, তারপর NumPy-তে

অংশ ৬ বলেছিল: কলাম orthonormal হলে projection সহজ, P = QQᵀ. তাহলে X-এর কলাম থেকে একই column space-এর একটি orthonormal basis বানাতে পারলে least squares সহজ হয়ে যায়। Gram–Schmidt প্রক্রিয়া ঠিক এটি করে, আর তার হিসাবের খাতা থেকেই আসে QR decomposition.

নিয়মটি এক বাক্যে: প্রতিটি নতুন কলাম থেকে আগের সব q-এর দিকের ছায়া বাদ দিন, যা রইল তার দৈর্ঘ্য ১ করুন। ছায়া বাদ দেওয়া মানে ঠিক অংশ ৬-এর রেখায় projection-এর r অংশ রাখা।

কেন solver এই পথ বেছে নেয়? X = QR বসালে normal equations হয় RᵀQᵀQRw = RᵀQᵀy ⇒ RᵀRw = RᵀQᵀy ⇒ Rw = Qᵀy (R invertible হলে)। XᵀX কখনো বানাতেই হয় না। Q দৈর্ঘ্য বদলায় না, তাই cond(R) = cond(X), অথচ cond(XᵀX) = cond(X)². X-এ ৭টি দশমিক অঙ্ক হারানোর ঝুঁকি থাকলে normal equations-এ ১৪টি, QR-এ ৭টিই। কাজের দাম: n-টি সারি আর p-টি কলামে normal equations মোটামুটি np² গুণ-যোগ নেয়, Householder QR প্রায় তার দ্বিগুণ; বেশিরভাগ ক্ষেত্রে স্থিরতার জন্য এই দামটুকু দেওয়া হয়। নিচের ব্লক বাজারঘরের ৪০টি কাল্পনিক অর্ডারে তিনটি পথ পাশাপাশি চালায়, তারপর প্রায় একই দুটি কলামের একটি কঠিন ক্ষেত্রে।

import numpy as np
from scipy.linalg import solve_triangular

rng = np.random.default_rng(42)
n = 40                                       # বাজারঘরের ৪০টি কাল্পনিক অর্ডার
km = rng.uniform(1, 25, n)                   # দূরত্ব
kg = rng.uniform(0.2, 15, n)                 # ওজন
X = np.column_stack([np.ones(n), km, kg])
y = X @ np.array([12, 2.5, 0.8]) + rng.normal(0, 4, n)   # ডেলিভারির মিনিট

def three_ways(X, y):
    w_ne = np.linalg.solve(X.T @ X, X.T @ y)   # projection-এর সূত্র: XᵀXw = Xᵀy
    w_ls = np.linalg.lstsq(X, y)[0]            # LAPACK-এর SVD-ভিত্তিক least squares
    Q, R = np.linalg.qr(X)                     # reduced QR: Q 40x3, R 3x3
    w_qr = solve_triangular(R, Q.T @ y)        # Rw = Qᵀy, নিচ থেকে ওপরে সমাধান
    return w_ne, w_ls, w_qr, R

w_ne, w_ls, w_qr, R = three_ways(X, y)
for name, w in [("normal", w_ne), ("lstsq", w_ls), ("QR", w_qr)]:
    print(f"{name:6s} w = {np.array2string(w, precision=4)}")
print(f"সর্বোচ্চ পার্থক্য: {max(abs(w_ne - w_ls).max(), abs(w_qr - w_ls).max()):.1e}")
print(f"cond(X) = {np.linalg.cond(X):.1f}  cond(XᵀX) = {np.linalg.cond(X.T @ X):.1f}"
      f"  cond(R) = {np.linalg.cond(R):.1f}")

road = 1.3 * km + rng.normal(0, 1e-6, n)     # রাস্তার দূরত্ব ≈ 1.3 × সোজা দূরত্ব
X2 = np.column_stack([np.ones(n), km, road])
true = np.array([10, 1, 1.5])
y2 = X2 @ true                               # গোলমাল নেই, তাই সঠিক উত্তর জানা
print(f"\ncond(X2) = {np.linalg.cond(X2):.1e}  cond(X2ᵀX2) = {np.linalg.cond(X2.T @ X2):.1e}")
for name, w in zip(["normal", "lstsq", "QR"], three_ways(X2, y2)[:3]):
    err = np.linalg.norm(w - true) / np.linalg.norm(true)
    print(f"{name:6s} w = {np.array2string(w, precision=4)}  আপেক্ষিক ভুল {err:.1e}")

normal w = [9.0774 2.5542 0.9562]
lstsq  w = [9.0774 2.5542 0.9562]
QR     w = [9.0774 2.5542 0.9562]
সর্বোচ্চ পার্থক্য: 7.3e-14
cond(X) = 50.4  cond(XᵀX) = 2542.9  cond(R) = 50.4

cond(X2) = 3.8e+07  cond(X2ᵀX2) = 1.4e+15
normal w = [10.      0.6569  1.7639]  আপেক্ষিক ভুল 4.3e-02
lstsq  w = [10.   1.   1.5]  আপেক্ষিক ভুল 5.9e-10
QR     w = [10.   1.   1.5]  আপেক্ষিক ভুল 2.8e-10

Gram–Schmidt নিজে লিখে NumPy-র সঙ্গে মেলানো। নিচের function ঠিক stepper-এর নিয়ম মানে। একটি ছোট পরিবর্তনও আছে: classical সংস্করণে ছায়া মাপা হয় মূল কলাম a_j থেকে, আর modified সংস্করণে এখন পর্যন্ত ছেঁটে রাখা v থেকে। গাণিতিকভাবে দুটি এক; float-এ নয়। পরীক্ষাটি একটি কঠিন matrix-এ: 1, t, t², …, t⁹-এর কলাম (polynomial feature, যা পর্ব ২৩ (Linear model গভীরে — basis function, interaction আর regularization path)-এ আসবে), যাদের অনেকগুলো প্রায় একই দিকে।

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

def gram_schmidt(A, modified=False):
    """A-এর কলাম থেকে orthonormal কলামের Q আর upper-triangular R"""
    m, k = A.shape
    Q, R = np.zeros((m, k)), np.zeros((k, k))
    for j in range(k):
        v = A[:, j].copy()
        for i in range(j):
            # classical: মূল কলামের ছায়া; modified: এখন পর্যন্ত বাকি অংশের ছায়া
            R[i, j] = Q[:, i] @ (v if modified else A[:, j])
            v -= R[i, j] * Q[:, i]           # q_i দিকের অংশ বাদ
        R[j, j] = np.linalg.norm(v)          # যা রইল তার দৈর্ঘ্য
        Q[:, j] = v / R[j, j]
    return Q, R

A = np.array([[1.0, 3.0], [2.0, 3.0], [2.0, 0.0]])
Q, R = gram_schmidt(A)
print("3·Q =\n", 3 * Q); print("R =\n", R)
Qn, Rn = np.linalg.qr(A)                     # LAPACK-এর Householder QR
print("NumPy R =\n", Rn)
s = np.sign(np.diag(Rn))                     # কলামের চিহ্ন মিলিয়ে নেওয়া
print("চিহ্ন মেলালে এক?", np.allclose(Qn * s, Q), np.allclose(s[:, None] * Rn, R))
print("QᵀQ = I?", np.allclose(Q.T @ Q, np.eye(2)), "  QR = A?", np.allclose(Q @ R, A))

V = np.vander(np.linspace(0, 1, 30), 10, increasing=True)   # 1, t, t², …, t⁹
print(f"\ncond(V) = {np.linalg.cond(V):.1e};  ‖QᵀQ − I‖ (যত ছোট তত ভালো):")
for name, Qx in [("classical GS", gram_schmidt(V)[0]),
                 ("modified GS", gram_schmidt(V, True)[0]),
                 ("np.linalg.qr", np.linalg.qr(V)[0])]:
    print(f"  {name:13s} {np.linalg.norm(Qx.T @ Qx - np.eye(10)):.1e}")

3·Q =
 [[ 1.  2.]
 [ 2.  1.]
 [ 2. -2.]]
R =
 [[3. 3.]
 [0. 3.]]
NumPy R =
 [[-3. -3.]
 [ 0. -3.]]
চিহ্ন মেলালে এক? True True
QᵀQ = I? True   QR = A? True

cond(V) = 3.5e+06;  ‖QᵀQ − I‖ (যত ছোট তত ভালো):
  classical GS  5.0e-02
  modified GS   5.4e-11
  np.linalg.qr  1.4e-15

কারণটি জ্যামিতিক। প্রায় একই দিকের দুটি কলামের ছায়া বাদ দিলে যা থাকে তা খুব ছোট, আর classical সংস্করণে আগের q-গুলোর ছোট ছোট ভুল সেই ছোট অংশে বড় অনুপাতে ঢুকে যায়। Modified সংস্করণ প্রতিটি ধাপের পর বাকি অংশ থেকে আবার মাপে, তাই আগের ভুলও সঙ্গে সঙ্গে ছেঁটে যায়। Householder পদ্ধতি একেবারে আলাদা পথে যায়: q-গুলো ধাপে ধাপে বানানোর বদলে X-কে একের পর এক Householder আয়না দিয়ে প্রতিফলিত করে upper-triangular বানায়; প্রতিটি আয়না orthogonal, তাই ভুল বাড়ে না। LAPACK-এর geqrf routine এই পথে চলে, আর np.linalg.qr সেটিই ডাকে। Givens rotation আরেকটি বিকল্প, যা একবারে একটি করে ঘর শূন্য করে; sparse matrix-এ বা একটি একটি করে সারি যোগ হওয়া ডেটায় কাজে লাগে।

৯ · ML-এ কোথায় কোথায় — layer, PCA, leverage আর সঠিক solver

১০ · AI-এর যুগে — matrix-এর হিসাব আর AI-এর দাবি যাচাই

AI assistant linear algebra-র derivation আত্মবিশ্বাসের সঙ্গে লেখে, আর চিহ্ন ও সূচক নির্ভুল দেখায় বলে ভুল চোখে পড়ে না। কোডেও সে এমন অভ্যাস দেয় যা এই পর্বের জ্যামিতির সরাসরি লঙ্ঘন।

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

• প্রতিটি matrix-এর shape লিখে দেখেছি গুণগুলো মেলে কি না (m×k, k×k, k×m)
• একটি ২×২ বা ৩×২ উদাহরণ হাতে করে AI-এর সংখ্যার সঙ্গে মিলিয়েছি
• দাবিগুলো কোডে যাচাই করেছি: allclose(P @ P, P), X.T @ r ≈ 0, Q.T @ Q ≈ I
• inv-এর জায়গায় solve/lstsq/QR ব্যবহার করেছি, আর cond(X) ছেপে দেখেছি
• ব্যবহার করা function আর parameter আমার installed version-এর অফিসিয়াল docs-এ আছে কি না দেখেছি
• Model-এর মান মাপার আগে ডেটা ভাগ করে test set আলাদা রেখেছি; train-এর residual দেখে ভালো-মন্দ সিদ্ধান্ত নিইনি, আর একটি সরল baseline-এর সঙ্গে তুলনা করেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি