numpy 编写梯度下降的例子 （二）

正规方程（Normal Equation）是线性回归在最小二乘准则下的解析解（闭式解，Closed-Form Solution）。与梯度下降通过迭代逼近极值点不同，正规方程利用微积分极值条件，一步直接求出最优参数。

---

一、正规方程的数学推导

1. 矩阵与符号定义

设数据集包含 \(n\) 个样本，\(d\) 个特征。定义包含常数偏置项的增广设计矩阵 \(X\)、真实标签向量 \(\mathbf{y}\) 与待求参数向量 \(\theta\)：

\[X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1d} & 1 \\ x_{21} & x_{22} & \cdots & x_{2d} & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{nd} & 1 \end{bmatrix}_{n \times (d+1)}, \quad \theta = \begin{bmatrix} w_1 \\ w_2 \\ \vdots \\ w_d \\ b \end{bmatrix}_{(d+1) \times 1}, \quad \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}_{n \times 1} \]

模型预测值为：

\[\hat{\mathbf{y}} = X\theta \]

2. 构造损失函数

采用无约束最小二乘误差（MSE 的常数系数并不影响极值点位置，推导中常取残差平方和的一半 \(\frac{1}{2}\sum e_i^2\) 以简化求导常数）：

\[J(\theta) = \frac{1}{2} \Vert{}\hat{\mathbf{y}} - \mathbf{y}\Vert{}_2^2 = \frac{1}{2} (X\theta - \mathbf{y})^T (X\theta - \mathbf{y}) \]

3. 展开二次型

利用转置性质 \((A - B)^T = A^T - B^T\) 展开：

\[J(\theta) = \frac{1}{2} \left( (X\theta)^T - \mathbf{y}^T \right) (X\theta - \mathbf{y}) = \frac{1}{2} \left( \theta^T X^T X \theta - \theta^T X^T \mathbf{y} - \mathbf{y}^T X \theta + \mathbf{y}^T \mathbf{y} \right) \]

因为 \(\theta^T X^T \mathbf{y}\) 和 \(\mathbf{y}^T X \theta\) 都是标量（维度为 \(1 \times 1\)），且标量的转置等于其自身：

\[(\theta^T X^T \mathbf{y})^T = \mathbf{y}^T (X^T)^T (\theta^T)^T = \mathbf{y}^T X \theta \]

因此中间两项相等，可直接合并：

\[J(\theta) = \frac{1}{2} \left( \theta^T X^T X \theta - 2\mathbf{y}^T X \theta + \mathbf{y}^T \mathbf{y} \right) \]

4. 矩阵求导与极值条件

对向量 \(\theta\) 求梯度。需要用到两条经典矩阵微分法则：

• 线性项求导： \(\nabla_\theta (\mathbf{a}^T \theta) = \mathbf{a}\)
• 二次型求导： 若矩阵 \(A\) 为对称矩阵，则 \(\nabla_\theta (\theta^T A \theta) = 2 A \theta\)。这里矩阵 \(A = X^T X\) 显然满足对称性：\((X^T X)^T = X^T (X^T)^T = X^T X\)。

对损失函数求偏导：

\[\nabla_\theta J(\theta) = \frac{1}{2} \left( 2 X^T X \theta - 2 X^T \mathbf{y} \right) = X^T X \theta - X^T \mathbf{y} \]

因为 \(J(\theta)\) 是凸函数，当梯度为零向量时取得全局极小值：

\[X^T X \theta - X^T \mathbf{y} = \mathbf{0} \]

整理得到正规方程的标准形式：

\[X^T X \theta = X^T \mathbf{y} \]

5. 求解最优参数

若矩阵 \(X^T X\) 满秩可逆，在方程两边左乘 \((X^T X)^{-1}\)：

\[\boxed{\theta^* = (X^T X)^{-1} X^T \mathbf{y}} \]

注： 若特征之间存在严重线性相关（多重共线性）或样本量小于特征数（\(n < d\)），\(X^T X\) 可能奇异不可逆。实际工程中会采用 Moore-Penrose 伪逆（np.linalg.pinv），或引入 \(L_2\) 正则化（岭回归，变为 \((X^T X + \lambda I)^{-1} X^T \mathbf{y}\)，保证矩阵严格正定可逆）。

---

二、几何视角：正交投影

正规方程在几何上代表正交投影：

• 预测值 \(\hat{\mathbf{y}} = X\theta\) 是设计矩阵 \(X\) 的列向量张成的超平面（记为 \(\operatorname{col}(X)\)）中的一个点。
• 真实值 \(\mathbf{y}\) 通常不在这个超平面上。
• 残差向量 \(\mathbf{e} = \mathbf{y} - X\theta\) 代表从点 \(\mathbf{y}\) 到超平面的连线。当该连线垂直（正交）于该平面时，残差模长最小。
• 正交条件意味着残差向量与 \(X\) 的每一列内积全为 0：

\[X^T \mathbf{e} = \mathbf{0} \implies X^T (\mathbf{y} - X\theta) = \mathbf{0} \implies X^T X \theta = X^T \mathbf{y} \]

---

三、正规方程 vs. 梯度下降

比较维度	正规方程 (Normal Equation)	梯度下降 (Gradient Descent)
求解机制	解析解：矩阵求逆一步到位	迭代法：沿负梯度逐步逼近
超参数调节	无需选择学习率 \(\eta\)；无需指定迭代轮数	需调优学习率 \(\eta\)、迭代次数、收敛容差
特征缩放需求	不需要特征缩放（对量纲差异不敏感）	高度依赖特征缩放（否则等高线扁平导致震荡）
计算时间复杂度	矩阵求逆复杂度约为 \(\mathcal{O}(d^3)\) 或 \(\mathcal{O}(d^{2.8})\)	每轮迭代复杂度约为 \(\mathcal{O}(k \cdot n \cdot d)\)
内存占用	需在内存中计算并存储 \(X^T X\)（\(\mathcal{O}(d^2)\)）	支持小批量（Mini-batch），内存开销极低
适用场景	特征维度小 (\(d \le 10,000\))、数据可全部放入内存	大规模特征 (\(d \gg 10^5\))、海量样本、流式数据
算法拓展性	仅适用于线性回归、岭回归等简单凸优化	普适于逻辑回归、深度神经网络等复杂非线性架构

---

四、NumPy 验证对比

用前文提到的两特征房价数据，直接对比两种求解方式的数值输出：

import numpy as np

# 样本数据 (面积 10m², 房间数, 价格万元)
raw_X = np.array([
    [5.0, 1.0],
    [7.0, 2.0],
    [9.0, 2.0],
    [11.0, 3.0],
    [13.0, 4.0]
])
y = np.array([102.0, 148.0, 179.0, 229.0, 275.0])
n = len(y)

# 构造增广矩阵
X = np.c_[raw_X, np.ones(n)]

# 1. 正规方程直接求解
# 优先使用 np.linalg.pinv (伪逆) 避免矩阵不可逆崩溃
theta_normal = np.linalg.pinv(X.T @ X) @ X.T @ y

# 2. 梯度下降求解 (学习率 0.005, 30000 轮)
theta_gd = np.zeros(3)
lr = 0.005
for _ in range(30000):
    grad = (2.0 / n) * (X.T @ (X @ theta_gd - y))
    theta_gd -= lr * grad

print("正规方程解析解 :", np.round(theta_normal, 4))
print("梯度下降近似解 :", np.round(theta_gd, 4))

运行结果：

正规方程解析解 : [17.5714  9.7143  4.0714]
梯度下降近似解 : [17.5714  9.7143  4.0714]

当特征维度 \(d\) 不大时，正规方程简单直接且没有调参负担；但在工业界和现代深度学习中，由于参数动辄数百万到数百亿，计算逆矩阵在算力与显存上完全不可行，梯度下降（以及 Adam、SGD 等变体）是无可替代的核心基石。

误差

导数确实只是当前点无限小邻域内的“瞬时变化率”，一旦离开这个点，导数本身就会改变。因此，直接拿固定步长去减，必然存在误差。

理解这个误差的来源以及算法如何解决它，需要从泰勒展开和迭代机制来看。

---

一、误差从何而来：一阶泰勒展开的近似截断

考虑一元函数 \(f(w)\)。当参数从 \(w\) 移动一个微小增量 \(\Delta w\) 时，根据泰勒展开（Taylor Expansion）：

\[f(w + \Delta w) = f(w) + f'(w) \Delta w + \frac{1}{2} f''(w) (\Delta w)^2 + \mathcal{O}((\Delta w)^3) \]

梯度下降选择的更新步长为：

\[\Delta w = -\eta f'(w) \]

将 \(\Delta w\) 代入泰勒展开式中：

\[f(w - \eta f'(w)) \approx \underbrace{f(w)}_{\text{原损失}} - \underbrace{\eta [f'(w)]^2}_{\text{希望下降的量}} + \underbrace{\frac{1}{2} \eta^2 f''(w) [f'(w)]^2}_{\text{产生的误差（二阶截断误差）}} \]

• 理想情况（局部线性）： 如果步长极小（\(\eta \to 0\)），高阶项 \(\eta^2\) 可以忽略不计，函数值确实下降了 \(\eta [f'(w)]^2\)。
• 现实情况（非线性曲率）： 只要步长不是无限小，就必定会受到二阶项（曲率 \(f''(w)\)）的干扰。这就是你感受到的“更新误差”。

---

二、既然有误差，为什么模型不会崩？

既然每一步都不精确，梯度下降能收敛主要依赖两个机制：

1. “盲人下山”的动态修正机制

梯度下降从来没有指望一步就算准全局最优解。

• 导数只在脚下的一小块区域有效。
• 我们沿着当下的斜率迈出很小的一步（参数更新）。
• 踩到新位置后，旧的瞬时变化率立刻作废；我们在新位置重新测量此时此刻的瞬时变化率，再决定下一步往哪迈。
  这种高频的“微调-重新采样-再微调”闭环，不断吞噬掉上一步产生的二阶误差。

2. 学习率 \(\eta\) 的调控作用

学习率 \(\eta\) 本质上是对“局部线性假设”成立范围的约束：

• 如果 \(\eta\) 适中偏小： 负项 \(-\eta [f'(w)]^2\) 远大于正向误差项 \(\frac{1}{2} \eta^2 f''(w) [f'(w)]^2\)，每一步总能保证函数值净下降。
• 如果 \(\eta\) 太大： 此时忽略二阶误差的代价太大，误差项超过了下降项，导致损失不降反升，甚至彻底震荡发散。

---

三、能消除这个误差吗？（牛顿法与二阶优化）

如果你希望消除这种“用瞬时变化率代替整体弯曲”的误差，数学上有现成的方案——牛顿法（Newton's Method）。

牛顿法直接把泰勒二阶项也纳入考量，让一阶导数与二阶导数共同决定步长：

\[w_{\text{new}} = w - \frac{f'(w)}{f''(w)} \]

在多维空间中，二阶导数对应的是 Hessian 矩阵（海森矩阵，包含所有二阶偏导）：

\[\theta_{\text{new}} = \theta - H^{-1} \nabla J(\theta) \]

牛顿法直接用抛物面去逼近曲面，更新时几乎没有线性截断误差，在极少数步数内即可收敛。

为什么深度学习通常不用牛顿法？

• 如果模型有 1 亿个参数，梯度向量只有 1 亿维（\(\mathcal{O}(d)\)）；
• 但海森矩阵的大小是 \(1\text{ 亿} \times 1\text{ 亿}\)（\(\mathcal{O}(d^2)\)），求逆矩阵的计算复杂度更是高达 \(\mathcal{O}(d^3)\)。

因此，现代机器学习做出的妥协是：宁愿承担每一步用瞬时变化率带来的微小近似误差，靠小步快跑（一阶梯度迭代）来换取极高的计算吞吐量。