线性回归进阶：从设计矩阵到概率解释与基函数

在基础篇中，我们已经从一条直线出发，学习了 MSE、最小二乘法、梯度下降、Python 实现与回归指标。但线性回归还有几层很重要的含义：

• 从代数看，它是在解一个矩阵方程。
• 从几何看，它是在把目标向量投影到特征张成的空间中。
• 从概率看，它是在高斯噪声假设下做极大似然估计。

这三种视角会回答一些基础篇没有展开的问题：多元回归的闭式解是怎么来的？矩阵不可逆时怎么办？为什么最小二乘偏偏使用平方误差？一条“线性”模型能不能拟合曲线？预测值周围的不确定性又从哪里来？

0. 阅读索引：哪些内容已经讲过

这是一篇承接文章。为了不把同一个知识点重复书写，建议按下面的索引查阅：

想学习的内容	对应文章
一元线性回归、MSE、手算斜率和截距、回归指标	线性回归：从直觉到最小二乘法
梯度、学习率、Batch/SGD/Mini-batch	梯度下降法：机器学习如何一步步找到最优解
欠拟合、过拟合、Ridge、Lasso、偏差与方差	正则化与过拟合：让模型不只记住训练集
设计矩阵、正规方程、秩、伪逆、基函数、极大似然、预测分布	本文

本文默认已经理解下面两个式子：

$$
\hat{y}=\mathbf{w}^{T}\mathbf{x}+b
$$

$$
\mathrm{MSE}=\frac{1}{N}\sum_{n=1}^{N}(\hat{y}_n-y_n)^2
$$

1. 先把偏置并入权重

多元线性回归通常写成：

$$
\hat{y}=w_1x_1+w_2x_2+\cdots+w_dx_d+b
$$

每次都单独写 $b$ 会让矩阵推导比较麻烦。可以给输入补一个恒为 1 的特征：

$$
\tilde{\mathbf{x}}=
\begin{bmatrix}
1\x_1\x_2\\vdots\x_d
\end{bmatrix},\qquad
\boldsymbol{\theta}=
\begin{bmatrix}
b\w_1\w_2\\vdots\w_d
\end{bmatrix}
$$

于是模型统一写成：

$$
\hat{y}=\boldsymbol{\theta}^{T}\tilde{\mathbf{x}}
$$

这里的“增广”没有改变模型，只是把偏置 $b$ 变成了一个普通权重。以后看到设计矩阵第一列全是 1，就知道它负责学习截距。

1.1 一个两特征的计算例子

假设用“学习时长 $x_1$”和“睡眠时长 $x_2$”预测考试分数：

$$
\hat{y}=5x_1+2x_2+35
$$

一个学生学习 4 小时、睡眠 7 小时，其增广输入与参数为：

$$
\tilde{\mathbf{x}}=
\begin{bmatrix}1\4\7\end{bmatrix},\qquad
\boldsymbol{\theta}=
\begin{bmatrix}35\5\2\end{bmatrix}
$$

因此：

$$
\hat{y}=\boldsymbol{\theta}^{T}\tilde{\mathbf{x}}
=35\times1+5\times4+2\times7=69
$$

增广形式的好处是：无论有多少特征，预测都可以写成一次向量内积。

2. 从一个样本走向设计矩阵

假设有 $N$ 个样本，每个样本有 $d$ 个特征。把每个增广输入放在一行，就得到设计矩阵（Design Matrix）：

$$
\mathbf{X}=
\begin{bmatrix}
1 & x_{11} & x_{12} & \cdots & x_{1d}\
1 & x_{21} & x_{22} & \cdots & x_{2d}\
\vdots & \vdots & \vdots & \ddots & \vdots\
1 & x_{N1} & x_{N2} & \cdots & x_{Nd}
\end{bmatrix}
$$

目标值组成列向量：

$$
\mathbf{y}=
\begin{bmatrix}y_1\y_2\\vdots\y_N\end{bmatrix}
$$

所有样本的预测就可以一次写成：

$$
\hat{\mathbf{y}}=\mathbf{X}\boldsymbol{\theta}
$$

2.1 用具体数字检查维度

给定三个样本，每个样本只有一个原始特征：

$$
x=[0,1,2],\qquad y=[1,3,5]
$$

加入偏置列后：

$$
\mathbf{X}=
\begin{bmatrix}
1&0\
1&1\
1&2
\end{bmatrix},\qquad
\boldsymbol{\theta}=
\begin{bmatrix}b\w\end{bmatrix}
$$

如果 $b=1,w=2$：

$$
\mathbf{X}\boldsymbol{\theta}=
\begin{bmatrix}
1&0\1&1\1&2
\end{bmatrix}
\begin{bmatrix}1\2\end{bmatrix}
=
\begin{bmatrix}1\3\5\end{bmatrix}
$$

这正好等于真实目标 $\mathbf{y}$。矩阵乘法并不是新的模型，它只是把三个样本的三次预测合并成了一次运算。

3. 正规方程是怎样得到的

为了让公式更整洁，把平方误差和写为：

$$
J(\boldsymbol{\theta})=
\frac{1}{2}|\mathbf{X}\boldsymbol{\theta}-\mathbf{y}|_2^2
$$

前面的 $\frac{1}{2}$ 不会改变最优点，只是求导后能与平方产生的 2 抵消。对参数向量求梯度：

$$
\nabla_{\boldsymbol{\theta}}J
=\mathbf{X}^{T}(\mathbf{X}\boldsymbol{\theta}-\mathbf{y})
$$

最优点的梯度为 0：

$$
\mathbf{X}^{T}(\mathbf{X}\boldsymbol{\theta}-\mathbf{y})=0
$$

整理得到正规方程（Normal Equation）：

$$
\mathbf{X}^{T}\mathbf{X}\boldsymbol{\theta}
=\mathbf{X}^{T}\mathbf{y}
$$

如果 $\mathbf{X}^{T}\mathbf{X}$ 可逆，则闭式解为：

$$
\boldsymbol{\theta}^{*}
=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}
$$

3.1 手算一次矩阵闭式解

继续使用：

$$
\mathbf{X}=
\begin{bmatrix}1&0\1&1\1&2\end{bmatrix},\qquad
\mathbf{y}=
\begin{bmatrix}1\3\5\end{bmatrix}
$$

先计算：

$$
\mathbf{X}^{T}\mathbf{X}
=
\begin{bmatrix}3&3\3&5\end{bmatrix},\qquad
\mathbf{X}^{T}\mathbf{y}
=
\begin{bmatrix}9\13\end{bmatrix}
$$

逆矩阵为：

$$
(\mathbf{X}^{T}\mathbf{X})^{-1}
=\frac{1}{6}
\begin{bmatrix}5&-3\-3&3\end{bmatrix}
$$

所以：

$$
\boldsymbol{\theta}^{*}
=\frac{1}{6}
\begin{bmatrix}5&-3\-3&3\end{bmatrix}
\begin{bmatrix}9\13\end{bmatrix}
=\begin{bmatrix}1\2\end{bmatrix}
$$

最终模型就是 $\hat{y}=2x+1$。

3.2 几何上发生了什么

$\mathbf{X}\boldsymbol{\theta}$ 一定落在 $\mathbf{X}$ 的列空间中。最小二乘是在这个空间里寻找离 $\mathbf{y}$ 最近的向量。最优残差：

$$
\mathbf{r}=\mathbf{y}-\mathbf{X}\boldsymbol{\theta}^{*}
$$

与设计矩阵的每一列正交，因此：

$$
\mathbf{X}^{T}\mathbf{r}=0
$$

把 $\mathbf{r}$ 代回去，就再次得到正规方程。这说明“令梯度为 0”和“做正交投影”其实是同一件事的代数与几何表达。

4. 为什么有时不能直接求逆

闭式解暗含一个条件：$\mathbf{X}^{T}\mathbf{X}$ 必须可逆，也就是设计矩阵的列需要线性无关。

下面几种情况会破坏这个条件：

• 某个特征可以由其他特征精确线性组合得到。
• 特征数量大于样本数量。
• 重复添加了相同的特征。
• 某个特征整列都是常数，并且已经有偏置列。

例如同时把房屋面积的“平方米”和“平方英尺”作为两个特征。二者满足固定换算关系，因此两列提供的是同一份信息。此时参数解可能不唯一：一个权重增加，另一个权重可以相应减少，而预测保持不变。

4.1 秩、伪逆与 SVD

矩阵中独立信息的数量称为秩（Rank）。如果 $\mathbf{X}$ 的列不满秩，$(\mathbf{X}^{T}\mathbf{X})^{-1}$ 就不存在。

常见处理方式有三种：

1. 删除重复或高度冗余的特征。
2. 使用基于 SVD 的最小二乘或 Moore-Penrose 伪逆 $\mathbf{X}^{+}$。
3. 使用 Ridge，在矩阵上加入 $\lambda\mathbf{I}$；详细推导见正则化与过拟合。

使用伪逆时，可以写成：

$$
\boldsymbol{\theta}^{*}=\mathbf{X}^{+}\mathbf{y}
$$

实践中通常不要显式计算逆矩阵。np.linalg.lstsq 会使用更稳定的数值方法直接求解最小二乘问题：

import numpy as np

X = np.array([
    [1.0, 0.0],
    [1.0, 1.0],
    [1.0, 2.0],
])
y = np.array([1.0, 3.0, 5.0])

theta, residuals, rank, singular_values = np.linalg.lstsq(
    X,
    y,
    rcond=None,
)

print("theta =", theta)  # [1. 2.]
print("rank =", rank)    # 2，说明两列线性无关

5. “线性”是对参数线性，不一定对输入线性

这是线性回归中最容易产生误解的一点。

普通直线模型：

$$
y(x)=w_0+w_1x
$$

无法描述明显弯曲的关系。但把输入先经过固定的非线性变换 $\phi_j(x)$，再对变换后的特征做线性组合：

$$
y(x,\mathbf{w})=
w_0+\sum_{j=1}^{M-1}w_j\phi_j(x)
=\mathbf{w}^{T}\boldsymbol{\phi}(x)
$$

只要模型对待学习参数 $\mathbf{w}$ 仍然是线性的，它依旧属于线性模型。

5.1 多项式基函数

取：

$$
\boldsymbol{\phi}(x)=[1,x,x^2,\ldots,x^M]^T
$$

模型就变成多项式回归：

$$
y(x)=w_0+w_1x+w_2x^2+\cdots+w_Mx^M
$$

例如 $x=2$，使用二次特征：

$$
\boldsymbol{\phi}(2)=[1,2,4]^T
$$

若权重为 $\mathbf{w}=[3,1.5,-0.25]^T$，则：

$$
y(2)=3\times1+1.5\times2-0.25\times4=5
$$

曲线对 $x$ 是非线性的，但对 $w_0,w_1,w_2$ 仍然是线性的。

5.2 高斯基函数

高斯基函数常写为：

$$
\phi_j(x)=\exp\left(-\frac{(x-\mu_j)^2}{2s^2}\right)
$$

• $\mu_j$ 决定“小山包”的中心位置。
• $s$ 决定它覆盖的宽度。

假设 $\mu_1=0,\mu_2=2,s=1$。当 $x=2$ 时：

$$
\phi_1(2)=e^{-2}\approx0.135,\qquad
\phi_2(2)=e^0=1
$$

这表示位于 $x=2$ 的基函数被完全激活，远处的基函数只产生较小响应。多个不同中心的高斯基函数可以拼出局部变化明显的曲线。

5.3 Sigmoid 基函数

Sigmoid 基函数可以写为：

$$
\phi_j(x)=\sigma\left(\frac{x-\mu_j}{s}\right)
=\frac{1}{1+\exp(-(x-\mu_j)/s)}
$$

$\mu_j$ 控制转折位置，$s$ 控制转折快慢。多个阶梯形响应的线性组合可以逼近较复杂的平滑函数。

5.4 基函数与可学习隐藏层的区别

这里的 $\phi_j$ 通常是预先选好的，训练只学习外层权重 $w_j$。MLP 中隐藏层的参数也参与学习，所以表达能力更强，但优化问题不再是简单的凸二次函数。神经网络的实现可继续阅读：MLP 多层感知机：从原理到 PyTorch 实现。

6. 为什么最小二乘等价于极大似然

前面只是把平方误差当作一个合理的距离。现在给它一个概率解释。

假设真实目标由模型预测加上随机噪声得到：

$$
t_n=\mathbf{w}^{T}\boldsymbol{\phi}(x_n)+\epsilon_n
$$

并假设噪声独立且服从零均值高斯分布：

$$
\epsilon_n\sim\mathcal{N}(0,\beta^{-1})
$$

$\beta$ 称为精度（Precision），它是方差的倒数：

$$
\beta=\frac{1}{\sigma^2}
$$

因此单个目标的条件分布为：

$$
p(t_n\mid x_n,\mathbf{w},\beta)
=\mathcal{N}\left(
t_n\mid \mathbf{w}^{T}\boldsymbol{\phi}(x_n),
\beta^{-1}
\right)
$$

假设样本噪声彼此独立，整个数据集的似然是：

$$
p(\mathbf{t}\mid\mathbf{X},\mathbf{w},\beta)
=\prod_{n=1}^{N}p(t_n\mid x_n,\mathbf{w},\beta)
$$

取负对数并去掉与 $\mathbf{w}$ 无关的常数：

$$
-\log p(\mathbf{t}\mid\mathbf{X},\mathbf{w},\beta)
=\frac{\beta}{2}
\sum_{n=1}^{N}
\left[t_n-\mathbf{w}^{T}\boldsymbol{\phi}(x_n)\right]^2
+\text{constant}
$$

因为 $\beta>0$，最大化似然就等价于最小化平方误差和：

$$
\arg\max_{\mathbf{w}}p(\mathbf{t}\mid\mathbf{X},\mathbf{w},\beta)
=
\arg\min_{\mathbf{w}}
\sum_{n=1}^{N}(t_n-y(x_n,\mathbf{w}))^2
$$

所以平方误差并不是凭空出现的：它对应“观测误差是独立同方差高斯噪声”这一建模假设。

6.1 用残差估计噪声大小

在极大似然解 $\mathbf{w}_{ML}$ 下，噪声方差估计为：

$$
\sigma_{ML}^2
=\frac{1}{N}\sum_{n=1}^{N}
\left[t_n-y(x_n,\mathbf{w}_{ML})\right]^2
$$

精度就是：

$$
\beta_{ML}=\frac{1}{\sigma_{ML}^2}
$$

例如三个真实值为 $[2.2,3.9,6.1]$，模型预测为 $[2,4,6]$，残差为 $[0.2,-0.1,0.1]$。于是：

$$
\sigma_{ML}^2
=\frac{0.2^2+(-0.1)^2+0.1^2}{3}
=0.02
$$

$$
\beta_{ML}=\frac{1}{0.02}=50
$$

残差越分散，估计出的噪声方差越大、精度越小。

7. 从一个预测值到预测分布

普通回归代码往往只返回一个点预测：

$$
\hat{t}^{\mathrm{new}}
=(\mathbf{w}^{\mathrm{ML}})^{T}
\boldsymbol{\phi}(x^{\mathrm{new}})
$$

这里用上标 $\mathrm{ML}$ 表示极大似然估计得到的参数，用 $\mathrm{new}$ 表示新的输入与目标。

在前面的高斯噪声假设下，还可以写出新目标的分布：

$$
p(t^{\mathrm{new}}\mid x^{\mathrm{new}},
\mathbf{w}^{\mathrm{ML}},\beta^{\mathrm{ML}})
=\mathcal{N}\left(
t^{\mathrm{new}}\mid
(\mathbf{w}^{\mathrm{ML}})^{T}
\boldsymbol{\phi}(x^{\mathrm{new}}),
(\beta^{\mathrm{ML}})^{-1}
\right)
$$

这句话包含两层信息：

• 分布中心是回归模型给出的预测值。
• 分布宽度由训练残差估计出的观测噪声决定。

若只考虑高斯观测噪声，一个近似的 $95%$ 预测范围是：

$$
\hat{t}^{\mathrm{new}}\pm1.96\sigma^{\mathrm{ML}}
$$

需要注意：这只是“代入极大似然参数”后的噪声范围，没有包含参数 $\mathbf{w}$ 本身的不确定性。样本很少或新输入远离训练数据时，这个范围可能过于自信。要同时建模参数不确定性，需要进一步学习贝叶斯线性回归。

8. 为什么平方损失下最优预测是条件均值

假设给定输入 $x$ 后，目标 $t$ 仍然存在随机性。我们需要选择一个预测 $y(x)$，使期望平方损失最小：

$$
\mathbb{E}[L]
=\int\left(y(x)-t\right)^2p(t\mid x),dt
$$

对 $y(x)$ 求导并令其为 0：

$$
\frac{\partial\mathbb{E}[L]}{\partial y(x)}
=2\int(y(x)-t)p(t\mid x),dt=0
$$

由于 $\int p(t\mid x)dt=1$，可以得到：

$$
y^*(x)=\int t,p(t\mid x)dt
=\mathbb{E}[t\mid x]
$$

因此，在平方损失下，最优点预测是条件均值。这也解释了线性回归为什么在建模“给定输入后的平均目标值”，而不是保证每个样本都落在回归线上。

顺带对比：绝对误差损失下的最优点预测是条件中位数。损失函数不同，“最优预测”的统计含义也不同。

9. 手动模拟完整项目流程

下面把全文连起来：我们用二次基函数拟合一条曲线，再从残差估计噪声，并为新输入给出点预测与噪声范围。

9.1 第一步：准备数据

import numpy as np

x = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
t = np.array([4.2, 1.1, 0.1, 0.9, 4.1])

数据大致满足 $t=x^2$，但带有少量噪声。

9.2 第二步：选择基函数

使用：

$$
\boldsymbol{\phi}(x)=[1,x,x^2]^T
$$

因此设计矩阵为：

Phi = np.column_stack([
    np.ones_like(x),
    x,
    x ** 2,
])

print(Phi)

[[ 1. -2.  4.]
 [ 1. -1.  1.]
 [ 1.  0.  0.]
 [ 1.  1.  1.]
 [ 1.  2.  4.]]

9.3 第三步：求最小二乘解

w_ml, _, rank, singular_values = np.linalg.lstsq(
    Phi,
    t,
    rcond=None,
)

print("w_ml =", w_ml)
print("rank =", rank)

w_ml 的三个元素分别对应常数项、一次项和二次项。rank=3 表示三个基函数列彼此独立。

本例会得到近似结果：

w_ml = [ 0.0229 -0.0400  1.0286]
rank = 3

也就是拟合函数约为：

$$
\hat{t}=0.0229-0.04x+1.0286x^2
$$

9.4 第四步：计算残差与噪声

t_hat = Phi @ w_ml
residual = t - t_hat

sigma2_ml = np.mean(residual ** 2)
sigma_ml = np.sqrt(sigma2_ml)
beta_ml = 1.0 / sigma2_ml

print("residual =", residual)
print("sigma_ml =", sigma_ml)
print("beta_ml =", beta_ml)

这里的 Phi @ w_ml 就是 $\boldsymbol{\Phi}\mathbf{w}_{ML}$，而 sigma2_ml 是残差平方的平均值。

输出约为：

residual = [-0.0171  0.0086  0.0771 -0.1114  0.0429]
sigma_ml = 0.0641
beta_ml = 243.0556

9.5 第五步：预测新输入

假设新输入为 $x^{\mathrm{new}}=1.5$：

x_new = 1.5
phi_new = np.array([1.0, x_new, x_new ** 2])

mean = phi_new @ w_ml
lower = mean - 1.96 * sigma_ml
upper = mean + 1.96 * sigma_ml

print("point prediction =", mean)
print("noise interval =", (lower, upper))

计算结果约为：

point prediction = 2.2771
noise interval = (2.1514, 2.4029)

完整计算链为：

原始输入 x
  -> 基函数 phi(x)
  -> 设计矩阵 Phi
  -> 最小二乘求 w_ml
  -> 残差估计 sigma 和 beta
  -> phi(x_new)^T w_ml 得到分布中心
  -> mean +/- 1.96 * sigma 得到观测噪声范围

这段流程把“模型怎么拟合”与“预测有多不确定”连接了起来。

10. 正规方程、伪逆还是梯度下降

方法	适合场景	优点	注意点
正规方程	特征数量较少且矩阵条件良好	一步得到解析解	不建议显式计算矩阵逆
lstsq / SVD / 伪逆	可能秩不足或追求数值稳定	能处理更一般的最小二乘问题	仍受大规模矩阵分解成本限制
梯度下降	样本和特征规模很大、数据可分批读取	易扩展到大数据和复杂模型	需要学习率、迭代次数与特征缩放
Ridge	多重共线性明显、希望参数更稳定	改善条件数并控制过拟合	需要选择正则化系数 $\lambda$

如果只是写普通 NumPy 代码，优先使用 np.linalg.lstsq，而不是自己写 inv(X.T @ X)。如果使用 scikit-learn，则让 LinearRegression、Ridge 等成熟实现负责数值求解。

11. 把知识链串起来

现在可以把线性回归的学习路径分成四层：

1. 预测形式：$\hat{y}=\mathbf{w}^{T}\mathbf{x}+b$。
2. 优化形式：最小化平方误差，可用解析解或梯度下降。
3. 表示形式：使用基函数把原始输入映射到新的特征空间。
4. 概率形式：高斯噪声下，最小二乘就是极大似然，残差还能估计预测噪声。

继续往后学习时，可以按问题选择文章：

• 想研究迭代优化，阅读梯度下降法。
• 想理解模型复杂度、Ridge、Lasso、MAP 与偏差-方差权衡，阅读正则化与过拟合。
• 想从线性分数走向分类概率，阅读逻辑回归。
• 想让基函数也由数据学习，阅读MLP 多层感知机。

12. 总结

线性回归看起来只是一条直线，但它同时连接了线性代数、优化、概率与统计决策：

• 增广向量把权重和偏置统一起来。
• 设计矩阵把所有样本的预测写成 $\mathbf{X}\boldsymbol{\theta}$。
• 正规方程来自平方误差梯度为 0，也对应几何上的正交投影。
• 秩不足时不能盲目求逆，应考虑 SVD、伪逆、删去冗余特征或 Ridge。
• 线性模型只要求对参数线性，借助基函数仍然可以拟合曲线。
• 高斯噪声假设把极大似然与最小二乘连接起来。
• 平方损失下的最优预测是条件均值，残差方差描述了观测噪声。

掌握这些内容后，线性回归就不再只是第一个入门算法，而是一套能反复迁移到逻辑回归、神经网络、贝叶斯模型和正则化方法中的基础语言。

参考资料

1. Christopher M. Bishop, Pattern Recognition and Machine Learning, Chapter 3.
2. 周志华，《机器学习》，第 3 章。
3. NumPy Documentation, numpy.linalg.lstsq.