10.1 線性迴歸完整重現#

本節目標#

• 用第 6 章的「正規方程」 \(A^T A \mathbf{x} = A^T \mathbf{b}\) 求出最小平方解，完整重現線性迴歸。
• 推導迴歸的目標函數梯度 \(\nabla f = 2A^T(A\mathbf{x} - \mathbf{b})\) ，理解「梯度就是矩陣運算」。
• 手寫梯度下降迴圈，畫出收斂曲線，並實驗學習率太大時的發散現象。
• 對照「解析解」與「迭代解」的結果，體會兩條路通往同一個答案。

迴歸是機器學習的起點，而它的核心完全是線性代數：**找到一個向量 \(\mathbf{x}\) ，讓 \(A\mathbf{x}\) 盡可能接近 \(\mathbf{b}\) **。這正是第 3 章講過的「投影」問題。

10.1.1 問題設定：資料從哪裡來#

我們用一個有雜訊的線性模型產生資料，模擬真實世界的觀測：

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

n, d = 100, 3            # 100 筆資料、3 個特徵
x_true = np.array([2.0, -1.5, 3.0])

# 設計矩陣 A：每列是一筆資料的特徵向量，第一行是常數項
X_raw = rng.normal(size=(n, d))
A = np.column_stack([np.ones(n), X_raw])     # 形狀 (100, 4)
b = A @ x_true + rng.normal(scale=0.5, size=n)  # 加上雜訊

print("A 的形狀:", A.shape)   # (100, 4)
print("b 的形狀:", b.shape)   # (100,)

注意一個重要的面向轉換：第 3 章我們問「 \(A\mathbf{x} = \mathbf{b}\) 有沒有解？」，這裡 \(A\) 是 \(100 \times 4\) ，方程式比未知數多，通常無精確解。於是退而求其次，問「哪個 \(\mathbf{x}\) 讓誤差最小？」——這就是最小平方問題。

10.1.2 解析解：正規方程#

定義目標函數（殘差平方和）：

對 \(\mathbf{x}\) 微分（這是矩陣微積分，結果先記住，第 10.2 節會大量使用同樣的手法）：

令梯度為零，得到正規方程 \(A^T A \mathbf{x} = A^T \mathbf{b}\) ，也就是第 6 章的結論： \(\mathbf{b}\) 投影到 \(A\) 的行空間，投影座標就是答案。

# 正規方程：A^T A x = A^T b
x_normal = np.linalg.solve(A.T @ A, A.T @ b)
print("正規方程解:", np.round(x_normal, 4))
# 預期接近 [2, -1.5, 3, 0(截距部分按設計而定)]

也可用 np.linalg.lstsq，它在 \(A^T A\) 奇異（共線性強）時仍穩定：

x_lstsq, *_ = np.linalg.lstsq(A, b, rcond=None)
print("lstsq 解:  ", np.round(x_lstsq, 4))

10.1.3 迭代解：手寫梯度下降#

解析解一次算完，但現實中（例如第 10.2 節的神經網路）目標函數非線性、寫不出解析解。此時用梯度下降：沿著負梯度方向一小步一小步走。

梯度我們已經推好了： \(\nabla f = 2A^T(A\mathbf{x} - \mathbf{b})\) 。整個演算法只需要矩陣乘法！

def grad(x):
    return 2 * A.T @ (A @ x - b)

def gradient_descent(lr, iters=60):
    x = np.zeros(A.shape[1])
    history = []
    for i in range(iters):
        history.append(0.5 * np.sum((A @ x - b) ** 2))
        x = x - lr * grad(x)
    return x, np.array(history)

x_gd, hist = gradient_descent(lr=0.005)
print("梯度下降解:", np.round(x_gd, 4))
print("兩解差距:  ", np.linalg.norm(x_gd - x_normal))

畫出收斂曲線——損失應該單調下降：

plt.figure(figsize=(6, 4))
plt.plot(hist, lw=2)
plt.xlabel("iteration"); plt.ylabel("loss")
plt.title("learning rate = 0.005")
plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

10.1.4 學習率的實驗：太小慢、太大爆#

學習率 \(\eta\) 決定每步走多遠。對二次目標函數可以證明： \(\eta < 1/\lambda_{\max}(A^TA)\) 才收斂（特徵值出現了！）。

lrs = [0.0005, 0.005, 0.02, 0.03]
plt.figure(figsize=(7, 4.5))
for lr in lrs:
    _, h = gradient_descent(lr, iters=60)
    plt.plot(h, label=f"lr={lr}", lw=2)
plt.yscale("log"); plt.xlabel("iteration"); plt.ylabel("loss")
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

# 檢查發散邊界
lam_max = np.linalg.eigvalsh(A.T @ A)[-1]
print("收斂門檻 1/λ_max =", 1 / lam_max)   # 約 0.03x，lr=0.03 已在邊緣

你會看到：lr=0.0005 下降太慢；lr=0.02 收斂漂亮；lr=0.03 直接指數爆炸（log 尺度下是直線上升）。這就是特徵值在告訴你步長的上限——線性代數不是理論裝飾，它預測了演算法的生死。

10.1.5 兩種解的對照與延伸#

一個值得動手的小實驗：把特徵縮放改成「標準化每個特徵」，觀察 \(A^T A\) 的條件數（第 6 章）下降、可容忍的學習率變大。條件數與收斂速度的關係是最佳化理論的經典結果。

與《用 Python 學機率統計》的銜接：那本書用統計觀點談迴歸（ \(t\) 檢定、 \(R^2\) 、信賴區間），本書補上幾何與代數觀點——迴歸就是投影，係數就是投影座標。兩個視角合起來，你才真正擁有這個工具。

小結與練習#

1. 正規方程 \(A^T A \mathbf{x} = A^T \mathbf{b}\) 是投影的代數形式；梯度下降是同一答案的迭代逼近。
2. 迴歸的梯度 \(2A^T(A\mathbf{x}-\mathbf{b})\) 全由矩陣乘法構成，這個模式會一路延伸到神經網路。
3. 學習率上限由 \(A^T A\) 的最大特徵值決定。

練習：

• (a) 對特徵做標準化後重跑梯度下降，比較 60 步內的損失下降幅度。
• (b) 把 scale=0.5 的雜訊改成 5.0，觀察兩種解是否仍一致。
• (c) 挑戰：實作「動量法」v = 0.9*v - lr*grad; x = x + v，比較收斂速度。