ক্যাপস্টোন — লিনিয়ার রিগ্রেশন সম্পূর্ণ গাণিতিকভাবে ডেরাইভ ও কোড করা

১ · মডেল ও লস ফাংশন

লিনিয়ার রিগ্রেশনের মডেল সহজ — প্রতিটি ইনপুট ভেক্টর $\mathbf{x} \in \mathbb{R}^d$-এর জন্য একটি সংখ্যা পূর্বাভাস করা: $$\hat y = \mathbf{w}^T\mathbf{x} + b$$ এটিই পাঠ ০১-এ দেখা একক-নিউরনের গণনা — এখানে ওজন ভেক্টর $\mathbf{w} \in \mathbb{R}^d$ ও বায়াস $b$ হলো শেখার মতো প্যারামিটার। $n$টি ডেটা পয়েন্ট $\{(\mathbf{x}_i, y_i)\}_{i=1}^n$ থাকলে, আমাদের লস ফাংশনMean Squared Error, MSEপ্রতিটি পূর্বাভাস ও প্রকৃত মানের পার্থক্যের বর্গের গড়। হলো: $$\mathcal{L}(\mathbf{w}, b) = \frac{1}{n}\sum_{i=1}^n (\hat y_i - y_i)^2 = \frac{1}{n}\sum_{i=1}^n (\mathbf{w}^T\mathbf{x}_i + b - y_i)^2$$ লক্ষ্য — এমন $\mathbf{w}, b$ খুঁজে বের করা যা $\mathcal{L}$-কে সর্বনিম্ন করে। এই একটি সমীকরণকে দুটি সম্পূর্ণ ভিন্ন উপায়ে সমাধান করা যায় — একটি ক্যালকুলাসে ভিত্তি করে সরাসরি সূত্র বসিয়ে (নরমাল ইকুয়েশন), আরেকটি ধাপে ধাপে গ্রেডিয়েন্ট অনুসরণ করে (গ্রেডিয়েন্ট ডিসেন্ট) — এবং দুটোই একই উত্তরে পৌঁছায়, কারণ $\mathcal{L}$ একটি কনভেক্স ফাংশন (পাঠ ১৬)।

২ · গ্রেডিয়েন্ট ডেরাইভ করা (পাঠ ১৪)

পাঠ ১৪-এর নিয়ম অনুযায়ী, $\mathcal{L}$-এর প্রতিটি প্যারামিটারের সাপেক্ষে পার্শিয়াল ডেরিভেটিভ বের করি। প্রথমে $b$-এর সাপেক্ষে (চেইন রুল, পাঠ ১৩ ব্যবহার করে): $$\frac{\partial \mathcal{L}}{\partial b} = \frac{1}{n}\sum_{i=1}^n 2(\hat y_i - y_i)\cdot\frac{\partial \hat y_i}{\partial b} = \frac{2}{n}\sum_{i=1}^n (\hat y_i - y_i)$$ কারণ $\partial \hat y_i/\partial b = 1$। এখন $\mathbf{w}$-এর সাপেক্ষে — এখানে $\partial \hat y_i/\partial \mathbf{w} = \mathbf{x}_i$ (পাঠ ১৪-এ দেখা $\nabla_{\mathbf{x}}(\mathbf{w}^T\mathbf{x}) = \mathbf{w}$-এর প্রতিসম সংস্করণ): $$\nabla_{\mathbf{w}} \mathcal{L} = \frac{1}{n}\sum_{i=1}^n 2(\hat y_i - y_i)\,\mathbf{x}_i = \frac{2}{n}\sum_{i=1}^n (\hat y_i - y_i)\,\mathbf{x}_i$$ সব ডেটা পয়েন্ট একসাথে ম্যাট্রিক্স আকারে লিখলে ($\mathbf{X}$ হলো $n \times d$ ম্যাট্রিক্স, প্রতিটি সারি একটি $\mathbf{x}_i^T$; $\hat{\mathbf{y}} = \mathbf{Xw} + b\mathbf{1}$, পাঠ ০৭): $$\nabla_{\mathbf{w}} \mathcal{L} = \frac{2}{n}\mathbf{X}^T(\hat{\mathbf{y}} - \mathbf{y}), \qquad \frac{\partial \mathcal{L}}{\partial b} = \frac{2}{n}\sum_i(\hat y_i - y_i)$$

৩ · ক্লোজড-ফর্ম সমাধান — নরমাল ইকুয়েশন (পাঠ ০৯)

যেহেতু $\mathcal{L}$ কনভেক্স (এটি একটি কোয়াড্রেটিক ফর্ম, পাঠ ১২), এর গ্লোবাল মিনিমাম ঠিক সেই বিন্দুতে যেখানে গ্রেডিয়েন্ট শূন্য (পাঠ ১৬)। গ্রেডিয়েন্ট শূন্য বসাই: $$\frac{2}{n}\tilde{\mathbf{X}}^T(\tilde{\mathbf{X}}\tilde{\mathbf{w}} - \mathbf{y}) = \mathbf{0} \implies \tilde{\mathbf{X}}^T\tilde{\mathbf{X}}\tilde{\mathbf{w}} = \tilde{\mathbf{X}}^T\mathbf{y}$$ এই সমীকরণকেই বলে নরমাল ইকুয়েশন (normal equations)। যদি $\tilde{\mathbf{X}}^T\tilde{\mathbf{X}}$ ইনভার্টিবল হয় (অর্থাৎ পূর্ণ র‍্যাঙ্ক, পাঠ ০৯ — এর মানে ফিচারগুলো একে অপরের সাথে সম্পূর্ণ রৈখিকভাবে নির্ভরশীল নয়, পাঠ ০৬), তাহলে উভয় পাশে $(\tilde{\mathbf{X}}^T\tilde{\mathbf{X}})^{-1}$ দিয়ে গুণ করে সরাসরি সমাধান পাওয়া যায়: $$\boxed{\tilde{\mathbf{w}} = (\tilde{\mathbf{X}}^T\tilde{\mathbf{X}})^{-1}\tilde{\mathbf{X}}^T\mathbf{y}}$$ এটিই ক্লোজড-ফর্ম সমাধান — কোনো ইটারেশন ছাড়াই, একটি ম্যাট্রিক্স ইনভার্স ও কয়েকটি গুণ দিয়ে সরাসরি সর্বোত্তম $\mathbf{w}, b$ বের করা যায়।

৪ · কেন এটি কাজ করে — গাউসিয়ান নয়েজের অধীনে MLE (পাঠ ৩০)

"স্কয়ারড এরর মিনিমাইজ করা" কি স্রেফ একটি সুবিধাজনক পছন্দ, নাকি এর পেছনে গভীর যুক্তি আছে? ধরি বাস্তব ডেটা তৈরি হয় এভাবে — সত্যিকারের রৈখিক সম্পর্ক প্লাস কিছুটা এলোমেলো নয়েজ: $$y = \mathbf{w}^T\mathbf{x} + b + \epsilon, \qquad \epsilon \sim \mathcal{N}(0, \sigma^2)$$ তাহলে (পাঠ ২৬-এর গাউসিয়ান PDF ব্যবহার করে) $y$, $\mathbf{x}$ দেওয়া থাকলে একটি গাউসিয়ান অনুসরণ করে: $y \mid \mathbf{x} \sim \mathcal{N}(\mathbf{w}^T\mathbf{x}+b, \sigma^2)$, অর্থাৎ $$P(y_i \mid \mathbf{x}_i; \mathbf{w}, b) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y_i - \hat y_i)^2}{2\sigma^2}\right)$$ $n$টি ডেটা পয়েন্ট iid ধরলে (পাঠ ৩০-এর মতো), likelihood হলো এদের গুণফল, ও log-likelihood যোগফল: $$\ell(\mathbf{w}, b) = \sum_{i=1}^n \log P(y_i \mid \mathbf{x}_i) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \hat y_i)^2$$ লক্ষ করুন — এই রাশির মধ্যে $\mathbf{w}, b$-এর উপর নির্ভরশীল একমাত্র অংশ হলো $-\frac{1}{2\sigma^2}\sum_i(y_i-\hat y_i)^2$ (প্রথম পদটি $\mathbf{w}, b$-নিরপেক্ষ ধ্রুবক)। তাই: $$\arg\max_{\mathbf{w},b} \ell(\mathbf{w},b) = \arg\max_{\mathbf{w},b}\left[-\sum_i (y_i-\hat y_i)^2\right] = \arg\min_{\mathbf{w},b}\sum_i(y_i-\hat y_i)^2$$ অর্থাৎ — গাউসিয়ান নয়েজের অধীনে ম্যাক্সিমাম লাইকলিহুড এস্টিমেশন ঠিক স্কয়ারড এরর মিনিমাইজ করার সমতুল্য। "কেন বর্গ (square), অন্য কোনো পাওয়ার নয়?" — এই প্রশ্নের সঠিক উত্তর হলো এটাই: নয়েজ গাউসিয়ান ধরে নিলে, বর্গ-এরর ঠিক সেই রাশি যা likelihood থেকে স্বাভাবিকভাবে বেরিয়ে আসে।

৫ · কোড: ক্লোজড-ফর্ম বনাম গ্রেডিয়েন্ট ডিসেন্ট — একই ডেটাসেটে

এখন এই দুটি সম্পূর্ণ ভিন্ন সমাধান পদ্ধতি — একটি সরাসরি সূত্র (ধারা ৩), আরেকটি পাঠ ১৭-এর ইটারেটিভ গ্রেডিয়েন্ট ডিসেন্ট (ধারা ২-এর গ্রেডিয়েন্ট ব্যবহার করে) — একই সিন্থেটিক ডেটাসেটে চালিয়ে দেখি তারা একই $\mathbf{w}, b$-এ পৌঁছায় কি না।

import numpy as np

rng = np.random.default_rng(42)

# --- সিন্থেটিক ডেটাসেট তৈরি: y = true_w . x + true_b + গাউসিয়ান নয়েজ ---
n, d = 200, 2
X = rng.normal(0, 1, size=(n, d))          # ফিচার ম্যাট্রিক্স (n x d)
true_w = np.array([3.0, -2.0])
true_b = 5.0
noise = rng.normal(0, 0.5, size=n)          # epsilon ~ N(0, 0.25)
y = X @ true_w + true_b + noise

# --- (A) ক্লোজড-ফর্ম সমাধান: নরমাল ইকুয়েশন ---
X_tilde = np.hstack([X, np.ones((n, 1))])   # বায়াসের জন্য কলাম যোগ
w_tilde_closed = np.linalg.inv(X_tilde.T @ X_tilde) @ X_tilde.T @ y
w_closed, b_closed = w_tilde_closed[:-1], w_tilde_closed[-1]

print("ক্লোজড-ফর্ম     -> w =", np.round(w_closed, 4), " b =", round(b_closed, 4))

# --- (B) গ্রেডিয়েন্ট ডিসেন্ট (পাঠ ১৭) ---
w_gd = np.zeros(d)
b_gd = 0.0
lr = 0.1
epochs = 1000

for epoch in range(epochs):
    y_hat = X @ w_gd + b_gd
    error = y_hat - y
    grad_w = (2 / n) * (X.T @ error)
    grad_b = (2 / n) * np.sum(error)
    w_gd -= lr * grad_w
    b_gd -= lr * grad_b

print("গ্রেডিয়েন্ট ডিসেন্ট -> w =", np.round(w_gd, 4), " b =", round(b_gd, 4))
print("প্রকৃত (true)      -> w =", true_w, " b =", true_b)

৬ · বোনাস — Ridge regularization (পাঠ ৩২-এর সাথে সংযোগ)

পাঠ ৩২-এ আমরা দেখেছিলাম প্রতিটি ওজনের উপর গাউসিয়ান প্রায়র বসালে MAP এস্টিমেশন $L2$ regularization-এ পরিণত হয়। এখানে সেই একই ধারণা প্রয়োগ করলে লস ফাংশন হয় $\mathcal{L}_{ridge} = \frac{1}{n}\|\tilde{\mathbf{X}}\tilde{\mathbf{w}} - \mathbf{y}\|^2 + \lambda\|\mathbf{w}\|_2^2$ (সাধারণত বায়াস $b$-কে regularize করা হয় না)। একই পদ্ধতিতে গ্রেডিয়েন্ট শূন্য বসিয়ে সমাধান করলে: $$\tilde{\mathbf{w}} = (\tilde{\mathbf{X}}^T\tilde{\mathbf{X}} + \lambda \mathbf{I}')^{-1}\tilde{\mathbf{X}}^T\mathbf{y}$$ যেখানে $\mathbf{I}'$ হলো identity ম্যাট্রিক্স, শুধু বায়াসের অবস্থানে $0$ (যাতে বায়াস regularize না হয়)। লক্ষণীয় — $\lambda \mathbf{I}'$ যোগ করা $\tilde{\mathbf{X}}^T\tilde{\mathbf{X}}$-কে সবসময় ইনভার্টিবল করে তোলে, এমনকি যখন এটি নিজে থেকে ইনভার্টিবল না-ও হয় (পাঠ ০৯-এর সেই সমস্যাটির একটি সরাসরি সমাধান)।

# উপরের একই X, y, X_tilde ব্যবহার করে Ridge regression
lam = 1.0
reg_matrix = lam * np.eye(d + 1)
reg_matrix[-1, -1] = 0   # বায়াসকে regularize করা হয় না

w_tilde_ridge = np.linalg.inv(X_tilde.T @ X_tilde + reg_matrix) @ X_tilde.T @ y
print("Ridge (λ=1.0) -> w =", np.round(w_tilde_ridge[:-1], 4), " b =", round(w_tilde_ridge[-1], 4))