最小二乘：把最优参数直接解出来

小象实战讲义 · 机器学习算法精讲

上一节末尾留了一个悬念：求线性模型的参数，其实可以不走梯度下降——把问题写成向量形式，直接解出来。本节兑现这个承诺。我们把 $N$ 个训练样本看成一个线性方程组，装进三个矩阵；把损失函数展开成一条三项式；再借助向量求导的阶梯公式，一步走到著名的正规方程与最小二乘解析解。推导不短，但每一步只用转置、乘法与求导——跟上节奏，你会发现”直接解”比想象中平易近人。

💡 核心导读

• 残差量的是哪个方向：y 方向的差与垂直距离之间只差一个由夹角决定的常数，两种量法不改变对 $\theta$ 的估计。
• 从方程组到矩阵三件套：$N$ 个样本就是 $N$ 个方程；设计矩阵 $A$（$N \times (n+1)$，行=样本、列=特征）、参数 $\theta$、目标 $\mathbf{y}$，一行矩阵乘法 $A\theta = \mathbf{y}$ 写尽全部方程。
• 展开与合并：$J(\theta) = \Vert A\theta - \mathbf{y}\Vert^2$ 展开成三项，中间两项靠”标量的转置还是自己”合并成 $-2\theta^TA^T\mathbf{y}$。
• 求导三级阶梯：从乘积法则 $\frac{d(U^TV)}{dx}$ 出发，经特例 $\frac{d(x^Tx)}{dx} = 2x$，推广到 $\frac{d(x^TBx)}{dx} = (B+B^T)x$——一路通往正规方程 $A^TA\theta = A^T\mathbf{y}$ 与解析解。

先回答一个问题：残差量的是哪个方向

回顾一下损失函数。对于线性拟合，cost function 是：

$$J(\theta) = \sum_{i=1}^{N} (h_\theta(x_i) - y_i)^2$$

画图看：给定输入 $x_i$，预测值 $h_\theta(x_i)$ 是函数图像上的一个点，$y_i$ 是真实目标点，两者的误差是竖直方向上的一小段——每一次都这样量。

可能有同学会问：实际上的点到直线的距离并不是这么算的，还可以求 perpendicular distance——垂直距离。那我们为什么不算”真正的距离”？

讲师的回答分两层。第一，求垂直距离当然也完全没有问题；第二，竖直方向的差与垂直距离之间满足一个常数倍数关系，这个常数由这条线的夹角决定。关键在于：给定了一个 $h_\theta$ 的方程，$\theta$ 就定了，夹角也随之固定——所以无论算直线距离还是垂线距离，两种目标对 $\theta$ 的”排序”完全一致，不会改变关于 $\theta$ 的估计。所以不必执着于垂直距离；几何上”点到线”讲究垂直，而我们每次实际算的，就是 y 轴方向上的差。

记号约定：以讲义为准

动手推导前，先声明记号。本课程所有章节统一 notation：$i$ 从 1 到大 $N$，代表 $N$ 组数据；每一个 $x_i$ 属于 $n$ 维向量空间——这一点每一节课都强调。课件里偶尔借用网上或其他来源的图，那些图的记号与我们的不完全一样；遇到不一致时，以讲义里写的为准。

这套记号从第 4 章一路用到现在：大 $N$ 记样本数，小 $n$ 记特征维度，$x_0 \equiv 1$ 收编偏置 $\theta_0$。

从数据到线性方程组

现在从头开始推：如果不用梯度下降，我们怎么解出最优估计？

给定数据 $D$：$(x_1, y_1), (x_2, y_2), \dots, (x_N, y_N)$。给定一个线性函数 $h_\theta(x) = \theta_0 + \theta_1 x_1 + \dots + \theta_n x_n$。如果我们”奢望”模型完美命中每一个样本——$x_1$ 输进来得到 $y_1$、$x_2$ 输进来得到 $y_2$、一直到大 $N$——那我们实际上列出了一个线性方程组：

$$\begin{aligned} \theta_0 + \theta_1 x_1^{(1)} + \theta_2 x_2^{(1)} + \dots + \theta_n x_n^{(1)} &= y_1 \\ \theta_0 + \theta_1 x_1^{(2)} + \theta_2 x_2^{(2)} + \dots + \theta_n x_n^{(2)} &= y_2 \\ &\vdots \\ \theta_0 + \theta_1 x_1^{(N)} + \theta_2 x_2^{(N)} + \dots + \theta_n x_n^{(N)} &= y_N \end{aligned}$$

$N$ 个方程、$n+1$ 个未知数。我们一直在做的事——不断调整 $\theta$ 使输入对上输出——本质上就是解这个线性方程组。

装进矩阵：三件套

方程组写起来啰啰嗦嗦、麻麻烦烦。为了把它收拢，我们给 $\theta_0$ 也配上一个 $x_0$——而 $x_0 \equiv 1$，取值不变，纯粹为了形式统一、矩阵看着完整。于是定义三个矩阵：

$$A = \begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} & \dots & x_n^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} & \dots & x_n^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_1^{(N)} & x_2^{(N)} & \dots & x_n^{(N)} \end{bmatrix}, \qquad \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix}, \qquad \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{bmatrix}$$

$A$ 称为设计矩阵。读它的方式：每一行代表一个数据（第一行是第一个 example，第二行是第二个……第 $N$ 行是第 $N$ 个）；每一列代表一个特征部分——第一列全是 1（$x_0$ 的部分），第二列是 $x_1$ 的部分，依此类推。

做一遍维度检查：$A$ 是 $N \times (n+1)$（第一列 $x_0$，加上 $x_1$ 到 $x_n$ 共 $n+1$ 列）；$\theta$ 是 $(n+1) \times 1$ 的列向量；$\mathbf{y}$ 是 $N \times 1$。于是：

$$N \times (n+1) \;\to\; (n+1) \times (n+1) \;\to\; (n+1) \times 1 \;\Longrightarrow\; N \times 1$$

完全符合矩阵乘法，整条方程组被一行式子收编：

$$A\theta = \mathbf{y}$$

损失函数的矩阵形式

我们想解出 $\theta$，但心里始终记着：这个方程组未必有解，我们真正要的是损失最小。以前的损失函数是每个样本误差的平方和：

$$J(\theta) = \sum_{i=1}^{N} (h_\theta(x_i) - y_i)^2$$

换成矩阵记号，它有一个更紧凑的写法。$A\theta$ 乘出来是一个 $N \times 1$ 的列向量，减去 $\mathbf{y}$ 还是列向量——每一行都是一个误差；误差的平方和，恰好就是这个向量的模的平方：

$$J(\theta) = \Vert A\theta - \mathbf{y}\Vert^2$$

这里的模取的是 $L_2$ 范数。范数家族——$L_0$、$L_1$、$L_2$、$L_\infty$——第 5 章讲距离测算时已经学过。而 $L_2$ 模的平方有个方便的性质：向量乘上自己的转置。对列向量 $r$，$r^Tr$ 就是各分量的平方和。

展开：三条路走成一条

现在把 $J(\theta)$ 展开成能求导的形式。先把模平方写成转置相乘：

$$J(\theta) = (A\theta - \mathbf{y})^T (A\theta - \mathbf{y})$$

把转置放进括号里——记住转置要倒序：$(A\theta)^T = \theta^T A^T$，于是：

$$J(\theta) = (\theta^T A^T - \mathbf{y}^T)(A\theta - \mathbf{y})$$

像普通代数一样展开（但矩阵分左乘右乘，顺序不能搞混）：

$$J(\theta) = \theta^T A^T A \theta - \theta^T A^T \mathbf{y} - \mathbf{y}^T A \theta + \mathbf{y}^T \mathbf{y}$$

看中间两项。$\theta^T A^T \mathbf{y}$ 的维度流转：$\theta^T$ 是 $1 \times (n+1)$，$A^T$ 是 $(n+1) \times N$，$\mathbf{y}$ 是 $N \times 1$——$(n+1)$ 相消、$N$ 相消，最后是 $1 \times 1$，一个标量。标量的转置等于它自己：

$$\mathbf{y}^T A \theta = (\theta^T A^T \mathbf{y})^T = \theta^T A^T \mathbf{y}$$

所以中间两项是同一个数，合并：

$$J(\theta) = \theta^T A^T A \theta - 2\theta^T A^T \mathbf{y} + \mathbf{y}^T \mathbf{y}$$

二次项、一次项、常数项——一条标准的三项式。

求导置零：三级阶梯

想求 $J(\theta)$ 的最优解，回到最朴素的微积分：给定一个函数，最优解基本停留在导数为 0 的位置。梯度下降的走法是先随机放一个 $\theta$，看梯度方向，不断调整，一步步逼近最完美的结果；而我们现在直接求”导数为 0”——殊途同归，但一步到位。

对 $\theta$ 求导，三项里最难的正是二次项 $\theta^T A^T A \theta$。讲师为此搭了一座三级阶梯（公式的严格证明属于数学课的范围，本课直接使用，细节可参考数学基础部分）。

第一级：向量版的乘积法则。 若 $U$、$V$ 都是关于 $x$ 的函数：

$$\frac{d(U^TV)}{dx} = \frac{dU^T}{dx}V + \frac{dV^T}{dx}U$$

这就是”两个变量相乘的组合函数怎么求导”在向量世界的翻版。

第二级：特例 $x^Tx$。 把 $U$、$V$ 都取成 $x$：

$$\frac{d(x^Tx)}{dx} = x + x = 2x$$

第三级：中间垫一个方阵 $B$。 把 $V$ 换成 $Bx$（$B$ 是方阵）：

$$\frac{d(x^TBx)}{dx} = Bx + B^Tx = (B + B^T)x$$

为什么 $B$ 必须是方阵？因为矩阵与它的转置要能相加，必须同形——行列相等，只有方阵满足。

现在把阶梯搭到我们的二次项上。取 $x = \theta$、$B = A^TA$：它恰好是方阵，而且对称——$(A^TA)^T = A^TA$，转置翻过来还是它自己。于是 $B + B^T = 2A^TA$：

$$\frac{\partial (\theta^T A^T A \theta)}{\partial \theta} = 2A^TA\theta$$

一次项里 $A^T\mathbf{y}$ 与 $\theta$ 无关，常数项求导为零。三项合起来：

$$\frac{\partial J(\theta)}{\partial \theta} = 2A^TA\theta - 2A^T\mathbf{y}$$

正规方程与解析解

令导数为零：

$$2A^TA\theta - 2A^T\mathbf{y} = 0$$

系数 2 直接约掉，得到正规方程（Normal Equation）：

$$A^TA\,\theta = A^T\mathbf{y}$$

若 $A^TA$ 可逆，两边左乘它的逆，最优参数一步解出：

$$\hat{\theta} = (A^TA)^{-1}A^T\mathbf{y}$$

不需要学习速率、不需要迭代、不需要初值——只要 $A$ 与 $\mathbf{y}$ 本身，矩阵运算直接给出最好估计的 $\hat{\theta}$。

数值验证：推导的每一环都亲手过一遍

用一小段纯 Python（手写矩阵运算，零第三方库）验证本节推导的每个环节：构造 5 个样本（$y$ 大致沿 $2x+1$ 分布、带手工噪声），解正规方程，检查三种写法的 $J$ 是否相等、解处梯度是否为零、以及扰动后损失是否上升。

# 数值验证：正规方程推导的每一步
# 数据：5 个样本（x 从 0 到 4），y 大致沿 y = 2x + 1 分布，带少量手工噪声
x = [0, 1, 2, 3, 4]
y = [1.0, 2.9, 5.2, 6.8, 9.1]

# 设计矩阵 A：第一列全 1（对应 x0 = 1 的约定），第二列为 x
A = [[1.0, xi] for xi in x]

# 矩阵运算三件套（纯 Python 手写，不依赖任何库）
def transpose(M):
    return [list(col) for col in zip(*M)]

def mat_mult(M, N):
    return [[sum(M[i][k] * N[k][j] for k in range(len(N)))
            for j in range(len(N[0]))] for i in range(len(M))]

def mat_vec(M, v):
    return [sum(M[i][k] * v[k] for k in range(len(v))) for i in range(len(M))]

def dot(v1, v2):
    return sum(a * b for a, b in zip(v1, v2))

# ---- 第一步：正规方程的两个构件 A^T A 与 A^T y ----
AT = transpose(A)
ATA = mat_mult(AT, A)
ATy = [dot(row, y) for row in AT]

print("设计矩阵 A（5 行 2 列，行=样本，列=[x0, x1]）:")
for row in A:
    print(" ", row)
print("A^T A =", ATA)
print("A^T y =", ATy)

# ---- 第二步：解 2x2 正规方程（克莱姆法则）----
det = ATA[0][0] * ATA[1][1] - ATA[0][1] * ATA[1][0]
theta0 = (ATA[1][1] * ATy[0] - ATA[0][1] * ATy[1]) / det
theta1 = (ATA[0][0] * ATy[1] - ATA[1][0] * ATy[0]) / det
theta = [theta0, theta1]
print(f"\ndet(A^T A) = {det}")
print(f"最小二乘解: theta0 = {theta0}, theta1 = {theta1}")

# ---- 第三步：三种写法的 J(theta) 必须相等 ----
# 写法一：残差平方和（定义）
y_hat = [theta0 + theta1 * xi for xi in x]
resid = [yi - h for yi, h in zip(y, y_hat)]
J1 = sum(r * r for r in resid)

# 写法二：||A*theta - y||^2（范数形式）
r_vec = [row[0] * theta[0] + row[1] * theta[1] - yi for row, yi in zip(A, y)]
J2 = dot(r_vec, r_vec)

# 写法三：theta^T A^T A theta - 2 theta^T A^T y + y^T y（推导中的展开式）
J3 = dot(theta, mat_vec(ATA, theta)) - 2 * dot(theta, ATy) + dot(y, y)

print(f"\n三种写法的 J(theta_hat):")
print(f"  残差平方和           = {J1:.10f}")
print(f"  ||A*theta - y||^2    = {J2:.10f}")
print(f"  展开式（三项代数式） = {J3:.10f}")

# ---- 第四步：梯度在 theta_hat 处应为零向量 ----
grad = [2 * g for g in [ATA[0][0] * theta[0] + ATA[0][1] * theta[1] - ATy[0],
                        ATA[1][0] * theta[0] + ATA[1][1] * theta[1] - ATy[1]]]
print(f"\n梯度 2A^T A theta - 2A^T y 在解处的取值: [{grad[0]:.3e}, {grad[1]:.3e}]")

# ---- 第五步：扰动证明这是最低点（碗底） ----
theta_p = [theta0 + 0.3, theta1 - 0.2]
y_hat_p = [theta_p[0] + theta_p[1] * xi for xi in x]
J_p = sum((yi - h) ** 2 for yi, h in zip(y, y_hat_p))
print(f"\n扰动后 theta' = ({theta_p[0]:.2f}, {theta_p[1]:.2f}) 的 J(theta') = {J_p:.6f}（大于 {J1:.6f}）")

print("\n逐点预测对比:")
for xi, yi, h, r in zip(x, y, y_hat, resid):
    print(f"  x={xi}, y={yi}, y_hat={h:.2f}, 残差={r:+.2f}")

设计矩阵 A（5 行 2 列，行=样本，列=[x0, x1]）:
  [1.0, 0]
  [1.0, 1]
  [1.0, 2]
  [1.0, 3]
  [1.0, 4]
A^T A = [[5.0, 10.0], [10.0, 30]]
A^T y = [25.0, 70.1]

det(A^T A) = 50.0
最小二乘解: theta0 = 0.98, theta1 = 2.01

三种写法的 J(theta_hat):
  残差平方和           = 0.0990000000
  ||A*theta - y||^2    = 0.0990000000
  展开式（三项代数式） = 0.0990000000

梯度 2A^T A theta - 2A^T y 在解处的取值: [0.000e+00, 0.000e+00]

扰动后 theta' = (1.28, 1.81) 的 J(theta') = 0.549000（大于 0.099000）

逐点预测对比:
  x=0, y=1.0, y_hat=0.98, 残差=+0.02
  x=1, y=2.9, y_hat=2.99, 残差=-0.09
  x=2, y=5.2, y_hat=5.00, 残差=+0.20
  x=3, y=6.8, y_hat=7.01, 残差=-0.21
  x=4, y=9.1, y_hat=9.02, 残差=+0.08

五个环节全部兑现：$A^TA = \begin{bmatrix}5 & 10 \\ 10 & 30\end{bmatrix}$、$A^T\mathbf{y} = (25, 70.1)$ 都能手算复核（$\sum x = 10$、$\sum x^2 = 30$、$\sum y = 25$、$\sum xy = 70.1$）；解 $\hat{\theta} = (0.98, 2.01)$ 恰好贴着生成数据的真实关系 $2x+1$；三种写法的 $J$ 严格相等（$0.099$）；解处梯度精确为零（连浮点毛边都没有）；往旁边挪一步，损失从 $0.099$ 涨到 $0.549$——碗底名副其实。

方程解不出来怎么办

最后把话说完整。$A\theta = \mathbf{y}$ 这个方程组，有的时候是解不出来的——要看小 $n$ 与大 $N$ 的数量关系，要看 $A$ 这个矩阵是否满秩的各种情况。

但注意我们的态度：我们并不关心怎么把它”精确解”出来，因为我们本来就不是真的在解这个方程。把 $\hat{\theta}$ 代回模型、把每个 $x_i$ 输进去算出来的，是一组 $\hat{y}_i$（读作 y hat——hat 就是帽子，表示对 $y$ 的估计）。我们真正期待的，是所有 $y_i - \hat{y}_i$ 的平方和在一切可能下最小。损失最小，才是我们的目标；“解方程”只是这个目标在理想情形下的巧合。

这套方法的名字：least mean square——中文就是最小二乘法。本节从代数上把它推了出来；下一节将换两副眼镜重新看它——几何上它是投影，概率上它连着极大似然。

📝 动手练一练

1. 标量合并的维度论证：设 $A$ 为 $5 \times 3$、$\theta$ 为 $3 \times 1$、$\mathbf{y}$ 为 $5 \times 1$。(a) 写出 $\theta^TA^T\mathbf{y}$ 的完整维度流转，确认它是标量；(b) 用”标量的转置等于自己”证明 $\mathbf{y}^TA\theta = \theta^TA^T\mathbf{y}$；(c) 思考：为什么 $J(\theta)$ 的展开式里，只有中间两项可以合并，而首尾两项不行？

2. 共线性手算：设计矩阵 $A = \begin{bmatrix}1 & 0 & 0 \\ 1 & 1 & 2 \\ 1 & 2 & 4\end{bmatrix}$（第三个特征恒等于第二个特征的 2 倍）。计算 $A^TA$ 与它的行列式，判断正规方程是否有唯一解；并说说第 9 章的岭回归打算怎么救场。

本章小结

• 残差的方向：竖直方向的差与垂直距离差一个由夹角决定的常数倍——给定 $\theta$ 夹角固定，两种量法不改变对 $\theta$ 的估计，所以直接用 y 方向的差。
• 矩阵三件套：$N$ 个样本即 $N$ 个方程；$A$（$N \times (n+1)$，行=样本、列=特征）、$\theta$（$(n+1) \times 1$）、$\mathbf{y}$（$N \times 1$），$A\theta = \mathbf{y}$ 一行写尽。
• 损失函数矩阵化与展开：$J(\theta) = \Vert A\theta - \mathbf{y}\Vert^2 = \theta^TA^TA\theta - 2\theta^TA^T\mathbf{y} + \mathbf{y}^T\mathbf{y}$，中间两项靠标量转置合并。
• 求导三级阶梯：乘积法则 $\to$ 特例 $2x$ $\to$ 推广 $(B + B^T)x$；$B = A^TA$ 对称，导数 $2A^TA\theta$，置零得正规方程 $A^TA\theta = A^T\mathbf{y}$，解出 $\hat{\theta} = (A^TA)^{-1}A^T\mathbf{y}$。
• 态度：方程组可能无解（要看 $n$、$N$ 与满秩），但我们不追求”精确解方程”——要的是 $\sum_i (y_i - \hat{y}_i)^2$ 最小。这套方法叫 least mean square，最小二乘法。

📋 行动清单

• 默写正规方程 $A^TA\theta = A^T\mathbf{y}$，并口头报出 $A$、$\theta$、$\mathbf{y}$ 三者的维度
• 不看讲义重推一遍三级阶梯：$\frac{d(U^TV)}{dx} \to \frac{d(x^Tx)}{dx} \to \frac{d(x^TBx)}{dx}$，并说清 $B$ 必须为方阵、$A^TA$ 对称的原因
• 手算本节数据的 $A^TA$ 与 $A^T\mathbf{y}$（$\sum x = 10$、$\sum x^2 = 30$、$\sum xy = 70.1$），与代码输出对账
• 预习下一节：最小二乘的几何视角（投影）与概率视角（高斯噪声）

—— 小象教研组