线性单元与梯度下降深度教程

从一条拟合直线到化学里最通用的学习算法

面向对象：刚入学、方向是「人工智能 × 化学」的研究生 需要的预备知识：一元微积分、向量点积、一点点 Python 不需要：任何机器学习基础、任何第三方库 配套代码：code/ 目录（纯 Python，零依赖，可直接运行） 配套插图：figures/ 目录（全部由图中的代码真实训练后生成） 手绘示意图：images/ 目录（手绘风格插图，提示词见 images/prompts/）

---

0. 写在前面

这份教程只讲两个东西：线性单元（linear unit）和梯度下降（gradient descent）。

线性单元是机器学习里最小的一个模型：把输入乘上一组权重、加起来、再加一个截距， 就得到一个预测值。梯度下降是机器学习里最小的一个优化器：算出往哪个方向走能减少误差， 然后朝那个方向走一小步，重复。

这两个东西加在一起，只有十几行代码。但你会发现，后面所有的东西—— 逻辑回归、神经网络、支持向量机、Transformer 的训练——用的都是同一套骨架。 区别只在于"模型更复杂"和"损失函数换了一个"。

但对化学方向的研究生来说，花时间深挖线性单元有四个非常实际的理由。

第一，它是你每天都在用的东西，只是名字不一样。 Hammett 方程、Arrhenius 方程、Clausius-Clapeyron 方程、Langmuir 吸附等温线、 Beer-Lambert 定律、Debye-Hückel 极限公式、Eyring 方程——这些你在物化课上背过的东西， 全都是"一个截距加几个斜率"。它们和线性单元是同一个数学对象。 线性单元做的事情只有一件：把"画一条最合适的直线"从手工变成自动。

第二，它是你自己的科研里最容易上手、也最容易被低估计的工具。 当你花三个月攒了 200 个数据点，第一件事不是上神经网络， 而是先跑一个线性单元。它能告诉你两件事：这个数据集的信息有多少落在一个线性方向上； 以及，如果线性模型只有 0.6 的 $R^2$，你手上就多了一个证据—— 这个体系里确实存在非线性效应（强关联、协同效应、溶剂化重构……）。 这两个数字比模型本身的预测值更值钱。

第三，它是你理解"训练"这个词的起点。 今天的模型动辄有几十亿参数，"梯度下降"这四个字却从来没变过。 在只有三个参数的线性单元上把梯度下降看懂，比在几十亿参数上猜它为什么收敛要容易得多。

第四，它是你识别"漂亮但假的结论"的第一道防线。 本教程第 5.3 节会给你看一个真实的例子：同一份数据、同一个模型， 随机划分数据集得到 $R^2 = 0.75$，按化学体系划分得到 $R^2 = -0.12$。 差的那 0.87 不是模型的能力，是评估方式的水分。 这个坑，化学文献里每年都有人踩。

本文的写法有五条纪律：

1. 每个公式都配一个几何图像或一段可运行的代码。 数学不是装饰， 它是让你知道代码里那一行 -= 为什么长这样。
2. 代码不用 NumPy。 用显式 for 循环手写向量运算，虽然慢， 但每一个乘法你都看得见。等你彻底搞懂了，换成 NumPy 是一个下午的事；反过来则很难。
3. 所有化学案例的数字都是真跑出来的。 图里的拟合线、准确率、残差， 全部来自 code/ 里的脚本。你在自己电脑上跑一遍，得到的数字应该完全一样。
4. 有测试。 包括用有限差分校验梯度、用正规方程校验梯度下降的结果， 确保代码没有隐藏的 bug。
5. 有负面结论。 教程花了大量篇幅讲线性单元什么时候会失败、 为什么失败、以及失败里藏着什么化学信息。

---

1. 背景：从「是/否」到「多少」

1.1 感知器留下的那道坎

1958 年，Rosenblatt 提出了感知器。它是历史上第一个能从数据中自动调节权重的模型。

但感知器的输出是一个非常粗暴的东西：

$$\hat{y} = \varphi(z) = \begin{cases} 1, & z > 0 \\ 0, & z \le 0 \end{cases}, \qquad z = \mathbf{w}\cdot\mathbf{x} + b$$

它只会说"是"或者"不是"。这个设定在分类任务上很自然，但它有两个致命的后果。

后果一：梯度消失了。 阶跃函数在除了 $z = 0$ 之外的每一个点上导数都是 0。而梯度下降的全部依据是"导数"。 导数恒为零意味着：你稍微动一动权重，输出一点变化都没有， 于是你无法知道该往哪个方向调整。

后果二：它丢掉了化学里最值钱的那部分信息。 假设你研究的是取代基对反应速率的影响。如果模型只告诉你 "这个取代基会让反应变快"，你几乎学不到新东西。 你真正想知道的是快多少——因为那个数字背后是电子效应、位阻效应、 以及一个可以写进论文讨论的物理图像。

图 1　左边：同样一套权重，两种输出方式。阶跃函数把连续的信号压成 0/1，同时也把梯度压成了 0；恒等函数保留了全部信息。右边：化学研究真正关心的问题。带隙的分类标签（半导体 / 绝缘体）只有一个比特，而带隙的数值 $E_g$ 有连续的物理含义。

1.2 1960：Widrow 和 Hoff 的关键转身

答案在两年后出现。1960 年，斯坦福的 Bernard Widrow 和他的学生 Ted Hoff 提出了一种叫 Adaline（Adaptive Linear Neuron）的模型。

Hoff 这个名字你也许在别处见过——他就是后来 Intel 4004 微处理器的发明者之一。

Adaline 和感知器的区别只有一个地方，但它是决定性的： 把输出层的激活函数从阶跃函数换成恒等函数。

$$\hat{y} = \varphi(z) = z = \mathbf{w}\cdot\mathbf{x} + b$$

现在输出是连续的了。既然输出连续，就可以定义"差了多少"； 既然能定义"差了多少"，就可以定义"误差有多大"； 既然误差是一个可导的函数，梯度就回来了。

Widrow 和 Hoff 给出的学习规则叫 LMS（Least Mean Squares，最小均方）， 后来更常见的名字是 delta 规则或者 Widrow-Hoff 规则：

$$\mathbf{w} \leftarrow \mathbf{w} + \eta\,(y - \hat{y})\,\mathbf{x}$$

这个式子读起来非常直白：如果预测低了，就把权重往输入的方向推一点； 如果预测高了，就往反方向推一点。 推的幅度由误差 $(y-\hat{y})$ 和学习率 $\eta$ 决定。

这条规则不是拍脑袋来的。它是今天所有深度学习模型使用的更新规则的一个特例—— 在本文第 3.3 节，我们会把它从损失函数的梯度里严格推出来。

顺便说一句：LMS 至今还在你的手机里工作。电话的回声消除、调制解调器的信道均衡， 用的都是 1960 年这个算法。

1.3 化学比任何学科都更需要线性单元

这里有一个很有意思的时间线对照。

年份	机器学习	化学
1937	—	Hammett 提出线性自由能关系 $\log(k/k_0) = \rho\sigma$
1943	McCulloch & Pitts 阈值神经元	—
1949	Hebb 学习律	—
1951	—	Taft 提出极性取代基常数 $\sigma^*$
1958	Rosenblatt 感知器	—
1960	Widrow & Hoff：LMS / Adaline	—
1962	—	Hansch 分析：$\log(1/C) = a\log P + b\sigma + c$
1964	—	Free-Wilson 加和模型
1969	Minsky & Papert 证明 XOR 问题	—
1974	Werbos 反向传播	—
1986	Rumelhart 等推广反向传播	—
1990s	核方法、SVM	QSAR / QSPR 大规模应用
2010s	深度学习	材料基因组、高通量计算
2017	Transformer	—
2020s	大语言模型	AlphaFold、GNoME、MACE 等

看第三行和第五行：化学家比计算机科学家早 20 年就用上了线性模型。

而且不是被动地用。Hammett 和 Taft 做的事情，用今天的话说就是"特征工程 + 线性回归"： 他们手工设计了一组能描述电子效应的数字（$\sigma$、$\sigma^*$）， 然后拟合了一个斜率（$\rho$、$\rho^*$），并且给这个斜率赋予了物理意义—— $\rho$ 的大小反映反应中心对电子效应的敏感程度，$\rho$ 的符号反映电荷的积累方向。

这件事说明线性模型有一个今天的深度网络很难提供的东西： 它的每一个参数都能被直接翻译成一句化学。 本教程第 5.6 节会用带隙案例演示怎么读这些参数。

1.4 线性单元藏在今天的哪里

线性单元没有停留在教科书里。下面这些地方，最后一步都是它。

场景	线性单元在哪里
图神经网络（GNN）的最后一层	把学到的图表示线性映射成性质预测值
深度学习模型的微调	冻结主干、只训最后一层线性头（linear probe）
注意力机制	用线性单元生成 query / key / value 三组投影
大语言模型的对齐层	隐状态到词表 logits 的一次矩阵乘法
分子性质预测的基线模型	描述符 × 权重 + 截距
可解释性研究	用线性探针读出中间表征里"编码了什么"

所以它不是一个"过时的玩具"，而是贯穿整个领域的一根主轴。 你以后读到的每一篇 AI 化学论文，模型再复杂， 最后一步几乎都是"若干个特征的加权求和"。

1.5 一句提醒

学完这份教程，千万不要得出"线性模型就够了"的结论。

正确的结论是：线性单元是你的第一把尺子。 它告诉你数据集里有多少信息是线性的、 哪些样本在犯难、模型错在哪里、以及下一步该往哪个方向加复杂度。 跳过这把尺子直接上深度学习，你会在两个地方吃亏： 不知道该期待多高的性能，以及在模型表现不好的时候不知道该改什么。

---

2. 模型：线性单元是什么

2.1 数学定义

线性单元接受一个 $n$ 维输入向量 $\mathbf{x} = (x_1, x_2, \dots, x_n)$， 输出一个实数 $\hat{y}$：

逐分量写法（最适合对应到代码里的 for 循环）：

$$\hat{y} = w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + b$$

向量写法（最适合做理论推导）：

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

增广写法（最适合写矩阵公式）：把 $b$ 吸收进权重向量， 在输入末尾补一个恒为 1 的分量，记 $\tilde{\mathbf{x}} = (1, x_1, \dots, x_n)$， $\boldsymbol\theta = (b, w_1, \dots, w_n)$，于是

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

三种写法完全等价。在本文里，凡是需要"从零实现"的时候用第一种， 凡是需要"推导公式"的时候用第二种或第三种。 这不是数学上的偷懒， 而是因为不同的写法对应不同的直觉。

三个名字的含义：

• $\mathbf{w}$：权重（weight）。它回答"这个描述符每增加一个单位，预测值变化多少"。
• $b$：偏置 / 截距（bias / intercept）。它回答"所有描述符都是 0 时，预测值是多少"。
• $\hat{y}$：预测值，用来和真实值 $y$ 比较。

2.2 几何图像

在一维输入的情况下（只有一个描述符），线性单元就是一条直线。

图 2　左边：一维输入的线性单元就是一条直线，橙色的竖线是每个样本的残差 $r_i = y_i - \hat{y}_i$。拟合的目标就是让这些竖线"总体上最短"。右边：把所有可能的 $(w, b)$ 组合画成一张平面，每个点对应一个模型，颜色表示它在这份数据上的误差。最优解是那个星号——碗底。

二维输入时，它是一张平面；三维输入时是一个超平面； $n$ 维输入时是一个 $n$ 维空间里的超平面。

但比"它是几维的"更重要的是另一件事：

线性单元把所有可能的世界分成了两类：在这个平面上，和不在这个平面上。

这既是它的全部能力，也是它的全部局限。如果真实世界里 $y$ 和 $\mathbf{x}$ 的关系 恰好接近一个超平面，它会表现得非常好；如果真实关系是弯的，它就只能画一条 "平均而言最不坏"的直线。

2.3 化学里到处都是线性单元

下表把你在化学课上见过的东西排在一起。注意看最后一列—— 它们无一例外都是"斜率的物理含义"。

名称	线性形式	纵坐标	横坐标	斜率的含义
Hammett 方程	$y = \rho\sigma$	$\log(k/k_0)$	取代基常数 $\sigma$	反应对电子效应的敏感度
Taft 方程	$y = \rho^*\sigma^*$	$\log(k/k_0)$	极性取代基常数 $\sigma^*$	对诱导效应的敏感度
Arrhenius 方程	$\ln k = \ln A - \dfrac{E_a}{R}\cdot\dfrac{1}{T}$	$\ln k$	$1/T$（或 $1000/T$）	$-\dfrac{E_a}{R}$，活化能
Eyring 方程	$\ln\dfrac{k}{T} = \ln\dfrac{k_B}{h} + \dfrac{\Delta S^\ddagger}{R} - \dfrac{\Delta H^\ddagger}{R}\cdot\dfrac{1}{T}$	$\ln(k/T)$	$1/T$	$-\dfrac{\Delta H^\ddagger}{R}$，活化焓
Clausius-Clapeyron	$\ln P = C - \dfrac{\Delta H_{\text{vap}}}{R}\cdot\dfrac{1}{T}$	$\ln P$	$1/T$	$-\dfrac{\Delta H_{\text{vap}}}{R}$，汽化焓
Langmuir 等温线	$\dfrac{1}{q} = \dfrac{1}{q_m K}\cdot\dfrac{1}{C} + \dfrac{1}{q_m}$	$1/q$	$1/C$	吸附平衡常数与饱和容量
Freundlich 等温线	$\ln q = \ln K_F + \dfrac{1}{n}\ln C$	$\ln q$	$\ln C$	吸附强度的不均匀性
Beer-Lambert	$A = \varepsilon c l$	吸光度	浓度	摩尔吸光系数
Debye-Hückel 极限	$\log\gamma_\pm = -A z_+ z_- \sqrt{I}$	$\log\gamma_\pm$	$\sqrt{I}$	离子强度效应
QSAR（Hansch）	$\log(1/C) = a\log P + b\sigma + c$	活性	$\log P, \sigma$	疏水性 / 电子效应权重

这张表想说的是一句话：

你在化学里学过的"作图法"，在机器学习里叫"一元线性回归"， 而"多元线性回归"就是线性单元。

但这里有一个陷阱。既然这些都是线性模型，那是不是说化学家早就会机器学习了呢？

不是。差别在数量级和自动化上：

• 化学传统的线性模型，特征是手工设计的（$\sigma$ 是人想出来的）， 系数是手工拟合的（在坐标纸上画一条线，或者按计算器）。 参数通常只有 2–5 个。
• 机器学习的线性单元，特征可以是几百上千个 （分子指纹、DFT 计算的描述符、图神经网络的嵌入）， 系数由梯度下降自动拟合，而且配套了一整套评估、诊断、防止自欺的流程。

所以本文真正要教你的，不只是"怎么拟合一条直线"—— 那你在第一堂实验课上就会了。本文要教的是那套流程： 怎么定义损失、怎么判断模型有没有学好、怎么知道它在什么地方会失灵。

---

3. 学习：梯度下降是怎么工作的

上一节我们知道了模型长什么样。但模型里有 $n+1$ 个未知参数， 怎么把它们算出来？这一节回答这个问题。

思路非常朴素，分三步：

1. 定一个标准：什么叫做"拟合得好"？——损失函数。
2. 找一个方向：往哪个方向调参数，损失会下降？——梯度。
3. 走一小步：走多远？——学习率。

3.1 第一步：定一个标准

对于第 $i$ 个样本，定义残差（residual）：

$$r_i = y_i - \hat{y}_i = y_i - (\mathbf{w}\cdot\mathbf{x}_i + b)$$

残差有正有负。如果直接把它们加起来，正负会互相抵消—— 一个把正误差和负误差都犯得很离谱的模型，会得到"总误差为 0"这个荒谬的结果。

所以我们要把残差变成正数。有两种常见做法：取绝对值，或者取平方。 我们选平方，定义均方误差（Mean Squared Error, MSE）：

$$\boxed{\;J(\mathbf{w}, b) = \frac{1}{2m}\sum_{i=1}^{m}\bigl(y_i - \hat{y}_i\bigr)^2\;}$$

其中 $m$ 是样本数。$J$ 越小，模型越好。训练的目标就是

$$(\mathbf{w}^*, b^*) = \arg\min_{\mathbf{w}, b} J(\mathbf{w}, b)$$

两个细节：

• 为什么分母是 $2m$ 而不是 $m$？ 纯粹是为了求导时舒服。 $\frac{1}{2}$ 和平方项的导数 $2$ 抵消，结果里不会留下多余的数字。 这不影响最优解的位置——把一个函数整体除以 2，最小值点不变。 （详细的讨论见附录 A.1。）
• 为什么是"均"方？ 除以 $m$ 让损失和样本数无关。 否则同一份数据你复制一遍，损失就变成四倍，学习率也得跟着改。

3.2 为什么用平方：从最大似然到最小二乘

"取平方"看起来是一个随手的选择。其实它有一个很硬的理由。

假设真实关系和你的模型之间差了一个随机误差：

$$y_i = \mathbf{w}\cdot\mathbf{x}_i + b + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \sigma^2) \ \text{独立同分布}$$

也就是说，误差是均值为 0、方差为 $\sigma^2$ 的高斯噪声。 化学上的测量误差、DFT 的计算误差、样本的天然涨落，往往可以用这个假设近似。

在这个假设下，给定 $\mathbf{w}$ 和 $b$，观测到 $y_i$ 的概率密度是

$$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)$$

整个数据集同时出现的概率（似然函数）是各样本概率的乘积。取对数， 再取相反数得到负对数似然：

$$-\log L = \frac{m}{2}\log(2\pi\sigma^2) + \frac{1}{2\sigma^2}\sum_{i=1}^{m}(y_i - \hat{y}_i)^2$$

第一项和 $\mathbf{w}, b$ 无关。所以

$$\arg\min_{\mathbf{w}, b} \bigl(-\log L\bigr) = \arg\min_{\mathbf{w}, b}\sum_{i=1}^{m}(y_i - \hat{y}_i)^2 = \arg\min_{\mathbf{w}, b} J(\mathbf{w}, b)$$

结论：在高斯噪声假设下，"最小化平方误差"和"最大化似然"是同一件事。

这不是巧合，而是最小二乘法在两百多年里一直有效的真正原因。 反过来说，如果你的误差明显不服从高斯分布—— 比如你的目标变量是计数（应该用泊松）、是概率（应该用交叉熵）、 或者数据里有大量离群点（平方会被极端值主导）—— 那么换成别的损失函数，往往比换模型更有效。

3.3 第二步：梯度

现在我们知道要最小化 $J$ 了。$J$ 是一个关于 $n+1$ 个参数的函数。 怎么找到它的最小值？

站在参数空间里的某一点，我们想知道往哪个方向走能最快地增加 $J$—— 然后往反方向走就行了。这个"最陡上升方向"就是梯度：

$$\nabla J = \left(\frac{\partial J}{\partial w_1},\ \dots,\ \frac{\partial J}{\partial w_n},\ \frac{\partial J}{\partial b}\right)$$

下面把这个梯度算出来。关键在于链式法则，思路是： 先看 $J$ 对预测值 $\hat{y}_i$ 的导数，再看 $\hat{y}_i$ 对每个参数的导数。

第一步，$J$ 对某个预测值的偏导。

$$\frac{\partial J}{\partial \hat{y}_i} = \frac{\partial}{\partial \hat{y}_i}\left[\frac{1}{2m}\sum_{k=1}^m (y_k - \hat{y}_k)^2\right] = \frac{1}{2m}\cdot 2(y_i - \hat{y}_i)\cdot(-1) = -\frac{r_i}{m}$$

求和里只有第 $i$ 项含有 $\hat{y}_i$，其余项导数为零。

第二步，预测值对参数的偏导。 因为 $\hat{y}_i = w_1 x_{i1} + \cdots + w_n x_{in} + b$， 所以

$$\frac{\partial \hat{y}_i}{\partial w_j} = x_{ij}, \qquad \frac{\partial \hat{y}_i}{\partial b} = 1$$

第三步，合并。 对某个 $w_j$，用链式法则：

$$\frac{\partial J}{\partial w_j} = \sum_{i=1}^{m} \frac{\partial J}{\partial \hat{y}_i}\cdot\frac{\partial \hat{y}_i}{\partial w_j} = \sum_{i=1}^{m}\left(-\frac{r_i}{m}\right)x_{ij}$$

于是我们得到全部两条梯度公式：

$$\boxed{\;\frac{\partial J}{\partial w_j} = -\frac{1}{m}\sum_{i=1}^{m} r_i\, x_{ij} \qquad (j = 1, \dots, n)\;}$$

$$\boxed{\;\frac{\partial J}{\partial b} = -\frac{1}{m}\sum_{i=1}^{m} r_i\;}$$

这就是整个教程最核心的两个公式。它们说的话非常朴素：

第 $j$ 个权重的梯度，等于"残差"和"第 $j$ 个特征"的（带负号的）内积的平均。

注意这个式子的结构：它把模型错在哪里（$r_i$）和输入长什么样（$x_{ij}$） 乘在一起。如果某个特征对每个样本都是正数，那么残差为正的样本会把这个权重往下推， 残差为负的样本会把它往上推，最后停在一个"总体平衡"的位置。

3.4 损失曲面：一片碗

把 $J$ 看成参数的函数，我们可以把它画出来。 只有 $w$ 和 $b$ 两个参数时，$J(w,b)$ 是一张三维曲面。

图 3　左：线性单元的损失函数在参数空间里是一个朝上的抛物面。右：把左图投影成等高线图。红色的箭头是 $-∇J$，也就是梯度下降每一步要走的方向；它永远垂直于你脚下的那条等高线。

这张图里最重要的一件事是形状：它是一个碗，而不是一片有多个坑的山地。

为什么？因为 $J$ 是参数的二次函数，它的二阶导数矩阵（Hessian）是

$$\mathbf{H} = \frac{1}{m} \begin{pmatrix} \mathbf{X}^\top \mathbf{X} & \mathbf{X}^\top \mathbf{1} \\ \mathbf{1}^\top \mathbf{X} & m \end{pmatrix}$$

这个矩阵是半正定的：对任意向量 $\mathbf{v}$，

$$\mathbf{v}^\top \mathbf{H} \mathbf{v} = \frac{1}{m}\sum_{i=1}^{m}\bigl(\mathbf{v}_{\mathbf{w}}\cdot\mathbf{x}_i + v_b\bigr)^2 \ge 0$$

因为它是若干个平方项的和，所以永远非负。

半正定意味着 $J$ 是一个凸函数，而凸函数有一个非常好的性质：

任何局部最小值都是全局最小值，而且最小值点构成一个凸集。

它直接给了我们两个保证：

1. 不存在"局部最优陷阱"。 你不用担心梯度下降会卡在一个小坑里出不来—— 因为这片地形里根本没有小坑。线性单元不会像神经网络那样对初始化敏感。
2. 不需要调初始化。 从任何起点出发，最终都会滑到同一个碗底附近。

深度学习里所有的痛苦（梯度消失、局部极小、鞍点、初始化敏感）， 在线性单元这里都不存在。这就是为什么它是最好的教学起点。

那唯一的困难是什么？ 看第 3.8 节——碗的形状。如果这个碗被拉成了一个又长又扁的 雪茄，梯度下降就会在窄的方向上反复横跳，走得极其缓慢。

3.5 第三步：梯度下降

有了梯度，算法就出来了。这是批量梯度下降（Batch Gradient Descent）：

输入：数据 (X, y)，学习率 η，迭代轮数 T

1. 初始化 w = 0，b = 0
2. 重复 T 次：
     a. 用当前参数算出所有样本的预测值 ŷ_i = w·x_i + b
     b. 算出所有残差 r_i = y_i - ŷ_i
     c. 计算梯度：
          ∂J/∂w_j = -(1/m) Σ_i r_i · x_ij
          ∂J/∂b   = -(1/m) Σ_i r_i
     d. 更新参数：
          w_j ← w_j - η · ∂J/∂w_j
          b   ← b   - η · ∂J/∂b
3. 输出 w, b

把梯度的负号代进去，更新步骤可以写成更好记的形式：

$$\boxed{\;w_j \leftarrow w_j + \frac{\eta}{m}\sum_{i=1}^{m} r_i\,x_{ij} \qquad b \leftarrow b + \frac{\eta}{m}\sum_{i=1}^{m} r_i\;}$$

现在再回头看第 1.2 节里 Widrow 和 Hoff 那条规则：

$$\mathbf{w} \leftarrow \mathbf{w} + \eta (y - \hat{y})\mathbf{x}$$

你会发现它就是上面这个式子在"每次只用一条样本"（$m=1$）时的样子。 Widrow-Hoff 规则不是发明出来的，它是从最小二乘损失里推出来的。

图 4　左：负梯度场。在每个位置，箭头都指向损失下降最快的方向，并且垂直于该处的等高线。右：从同一个起点出发，三种学习率走出的三条路。浅蓝色（η = 0.05）步子太小，50 步只走了一小段；深蓝色（η = 0.45）几乎直奔碗底；红色（η = 1.75）在窄方向上反复横跳。

3.6 学习率：三种命运

学习率 $\eta$ 是梯度下降唯一真正重要的超参数。它决定了每一步走多远。

图 5　同一片损失曲面上的四条轨迹，只改了学习率。太大和太小的失败方式完全不同：太小时慢慢挪，太大时越过碗底、落到对面的碗壁上，然后越飞越远。

可以用一个一维例子把这件事讲清楚。假设损失是 $J(w) = \frac{1}{2}aw^2$ （碗的"陡峭程度"由 $a$ 决定），梯度是 $aw$，那么更新是

$$w_{t+1} = w_t - \eta\, a\, w_t = (1 - \eta a)\,w_t$$

每一步都要乘上因子 $(1-\eta a)$。于是 $w_t = (1-\eta a)^t w_0$，分三种情况：

条件	因子 $\lvert 1-\eta a\rvert$	行为
$0 < \eta < 1/a$	小于 1，且为正	单调收敛，没有过冲
$1/a < \eta < 2/a$	小于 1，但为负	振荡收敛，在碗底两侧来回跳
$\eta > 2/a$	大于 1	发散，每一步都跳得更远

所以稳定性的硬条件是

$$\boxed{\;\eta < \frac{2}{a}\;}$$

推广到多参数情况，$a$ 换成 Hessian 矩阵的最大特征值 $\lambda_{\max}$：

$$\eta < \frac{2}{\lambda_{\max}}$$

这个不等式解释了一件在实践中非常重要的事：

如果你的数据没有做标准化，最大的特征值会被数值最大的那一列主导， 于是学习率的上界被压得极低，训练慢得让人绝望。

图 6 是在真实化学数据上做的学习率扫描。

图 6　左：73 个化合物、2 个描述符、参数已标准化。四条损失曲线对应四个学习率。右：把每个学习率跑 200 轮后的损失画成柱状图。绿色是稳定收敛区，橙色是缓慢区，红色是发散区。

一个有意思的细节：标准化之后，稳定区间的上界是可以估计的。 标准化让 $X^\top X / m$ 的对角线元素全都等于 1，于是 $\lambda_{\max} \le \operatorname{tr}(X^\top X/m) = n+1$， 因此 $\eta < 2/(n+1)$ 是一个理论上的安全保证。 对这个 2 特征的数据集，$2/3 \approx 0.67$； 而实测的发散点出现在 $\eta \approx 1.65$——理论保证是保守的， 但它的方向是对的。实践建议：标准化之后，学习率从 0.01–0.1 起试。

3.7 三种梯度：批量、随机、小批量

上面算法里的"一次用全部 $m$ 个样本算梯度"，是最原始的做法。 在实践中还有两个变体，区别只在"每次用多少样本"。

名字	每次用多少样本	每个 epoch 更新几次	梯度质量	单次计算成本
批量梯度下降（BGD）	全部 $m$ 个	1 次	精确	最高
小批量梯度下降（mini-batch）	$B$ 个（通常 16–256）	$m/B$ 次	有噪声但可控	中等
随机梯度下降（SGD）	1 个	$m$ 次	噪声很大	最低

图 7　左：三种梯度下降在同一份化学数据上的损失曲线。右：真正的计算代价对比——"达到 J < 2.0 所需的总梯度计算次数"。

这里有一个反直觉但非常重要的结论。

在你以前可能读过的科普文章里，故事通常是"SGD 比批量梯度下降快得多"。 但在本教程这份 73 个样本、4 个特征的数据上，批量下降反而是最省的： 它只用了 27 次梯度计算就达到目标，而随机下降用了 292 次。

原因是：批量下降每一步都用上了全部信息，方向最准； 随机下降一步只看一个样本，方向很抖，需要用很多步把噪声平均掉。

SGD 的优势来自"数据量大"，而不是"算法更聪明"。 当你有 100 万个样本时，一次全量梯度要扫描 100 万次， 而随机下降一步只花 1 次——那时它才真的不可替代。

对线性单元 + 中等规模化学数据集（几百到几千个样本）来说， 批量下降或者直接用正规方程，通常就是最优选择。

3.8 收敛速度：条件数才是真正的敌人

梯度下降能收敛，不代表它收敛得快。这一节解释为什么有时候它会慢得让人绝望。

回到那个一维的例子上，收敛速度由因子 $(1-\eta a)$ 决定。 多维情况下，每个方向有各自的"陡峭程度"，对应 Hessian 的不同特征值。 取最优学习率 $\eta = 2/(\lambda_{\max}+\lambda_{\min})$ 时，收敛速率由

$$\boxed{\;\kappa = \frac{\lambda_{\max}}{\lambda_{\min}}\;}$$

控制。$\kappa$ 叫做条件数（condition number）。它与迭代次数的关系大致是

$$t \sim \kappa \log\frac{1}{\varepsilon}$$

也就是说：条件数翻多少倍，需要的迭代次数就大致翻多少倍。

$\kappa$	等高线的形状	梯度下降的行为
$1$	正圆	一步就走到圆心附近
$10$	略扁的椭圆	十几步收敛
$10^3$	很扁的椭圆	几千步
$10^6$	一根针	几百万步，实际等于跑不动

而化学数据的条件数通常极其糟糕，因为化学描述符的量纲差异极大：

图 8　同一份带隙数据，只用"电负性差"和"摩尔质量"两个描述符。左：不缩放时条件数是 228 559，等高线被拉成一根细长的雪茄，梯度下降在窄方向上反复横跳，400 步之后还在挣扎。右：标准化之后条件数降到 1.9，等高线接近正圆，几步就走到碗底。

最漂亮的例子来自化学自己的传统。看这张图：

图 9　同一组 Arrhenius 数据，三种横坐标表示。(a) 损失曲线；(b) 三种表示的条件数，从 5×10⁷ 一路降到 1；(c) 学到的活化能随迭代的演化——只有标准化之后它才稳定地收敛到文献值附近。

这三个数字值得盯着看一会儿：

横坐标	条件数 $\kappa$	5000 轮后的 $R^2$
直接用 $1/T$	$5.1\times 10^7$	0.0002
用 $1000/T$	$6.1\times 10^3$	0.889
标准化后的 $1000/T$	$1$	0.9997

化学教材里那句"用 $1000/T$ 作图"，翻译成机器学习的话就是"把特征缩放到量纲相当的尺度"。

$1/T$ 的数值大约是 0.003，平方之后是 $10^{-5}$，梯度小得可怜； 乘 1000 之后数值变成 3 左右，梯度回到正常量级； 再标准化，条件数降到 1，问题变得"圆"，几步就收敛。

一百年前做图的人靠直觉躲开的坑，今天叫条件数。

实践上怎么办？ 三个办法，按性价比排序：

1. 标准化（z-score）：每一列减去均值、除以标准差。 这是最常用、最有效的一步。代价是权重的物理量纲被破坏了， 解读时要还原（见第 5.6 节）。
2. 换一个量纲合适的特征：就像用 $1000/T$ 而不是 $1/T$。 这一步同时保留了物理含义，是化学里最优雅的做法。
3. 用更好的优化器：动量法、Adam 等（见第 3.10 节）。 它们能缓解病态问题，但不能解决它。

一条铁律：在你调任何超参数之前，先做标准化。

3.9 解析解：既然能一步算出来，为什么还要迭代

线性单元有一个特殊性：它的最小二乘解有闭式公式，不需要迭代。

把截距吸收进权重，记增广矩阵 $\tilde{\mathbf{X}}$（第一列全是 1）， 损失对 $\boldsymbol\theta$ 的梯度是 $\frac{1}{m}\tilde{\mathbf{X}}^\top(\tilde{\mathbf{X}}\boldsymbol\theta - \mathbf{y})$。 在最优点梯度为零，于是

$$\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}\,\boldsymbol\theta = \tilde{\mathbf{X}}^\top\mathbf{y}$$

如果 $\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}$ 可逆，就得到正规方程：

$$\boxed{\;\boldsymbol\theta^* = \left(\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}\right)^{-1}\tilde{\mathbf{X}}^\top\mathbf{y}\;}$$

在本文的带隙案例上，正规方程和梯度下降给出的参数最大差异是 $5.6\times10^{-9}$—— 它们是同一个答案，只是一个用解析法、一个用迭代法得到。 但耗时差了很多：正规方程 0.0001 秒，梯度下降 0.76 秒。

那为什么还要学梯度下降？ 四个理由：

1. 复杂度。正规方程要算矩阵求逆，复杂度是 $O(n^3)$。 描述符有几百上千个时（在化学里非常常见），这一步会很吃力。
2. 不可逆的情况。当特征数超过样本数，或者特征之间高度共线时， $\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}$ 是奇异的，正规方程直接失效。 岭回归（见 3.10 节）能救它，但那时你已经不在用纯正规方程了。
3. 没有解析解的模型。逻辑回归、神经网络、支持向量机都没有闭式解， 但它们的更新规则和这里一模一样。梯度下降是通往深度学习唯一的路。
4. 数据流。SGD 可以处理源源不断到来的数据，不需要把整个数据集放进内存。

一句话总结：在这个规模上正规方程更快；但梯度下降是通用语言。

3.10 走得更聪明：动量、自适应步长、正则化

标准梯度下降有三个可以改进的地方，每一个都对应一个常见的扩展。

改进一：动量（Momentum）——记住上一段路的方向。

在又长又扁的碗里，梯度下降会在窄方向上反复横跳。动量的做法是 把历史梯度做一个指数加权平均：

$$\mathbf{v}_t = \beta\,\mathbf{v}_{t-1} - \eta\,\nabla J(\boldsymbol\theta_{t-1}), \qquad \boldsymbol\theta_t = \boldsymbol\theta_{t-1} + \mathbf{v}_t$$

$\beta$ 通常取 0.9。它的效果是：在方向一致的方向上加速，在来回横跳的方向上抵消。 在本文的测试里，动量让一个病态问题达到同样损失所需的迭代次数显著减少。

改进二：自适应步长——每个参数有自己的学习率。

AdaGrad、RMSProp、Adam 都属于这一类。它们记录每个参数梯度平方的历史， 然后用它来给每个参数单独缩放步长：

$$\theta_j \leftarrow \theta_j - \frac{\eta}{\sqrt{G_j} + \epsilon}\,\frac{\partial J}{\partial \theta_j}, \qquad G_j = \sum_{t} \left(\frac{\partial J}{\partial \theta_j}\right)^2$$

直觉是：一直梯度很大的参数，说明它在陡方向上，步子该小一点； 一直梯度很小的参数，说明它在平方向上，步子该大一点。 这几乎自动地解决了条件数问题，所以 Adam 是今天深度学习里最常用的优化器。

改进三：正则化——给参数加一点"体重"。

在损失里加一项 $\frac{\lambda}{2}\|\mathbf{w}\|^2$，得到

$$J_{\text{ridge}}(\mathbf{w}, b) = \frac{1}{2m}\sum_{i=1}^{m}(y_i - \hat{y}_i)^2 + \frac{\lambda}{2}\sum_{j=1}^{n}w_j^2$$

对应的梯度只是多了一项 $\lambda w_j$，更新规则变成

$$w_j \leftarrow w_j + \frac{\eta}{m}\sum_{i=1}^{m} r_i x_{ij} - \eta\lambda w_j$$

这个式子的效果很直观：每一步除了往减少误差的方向走，还额外往 0 的方向拉一点点。 这叫做岭回归（Ridge Regression）或 $L_2$ 正则化。

它的用处是压住共线性造成的参数爆炸。但注意：正则化解决的是"参数不稳定"， 不是"描述符缺了物理"。 在带隙案例里，$\lambda$ 从 0 加到 100， 训练 $R^2$ 从 0.79 单调掉到 0.02，模型只是变得更糟—— 因为那个模型的瓶颈根本不是过拟合。

---

4. 纯 Python 实现

4.1 设计原则

code/ 目录里的所有算法代码不使用 NumPy，也不使用任何第三方库， 只用 Python 标准库里的 math 和 random。

这不是为了自虐，而是三个具体的原因：

1. 每一个数学符号都能在代码里找到对应的那一行。 NumPy 的 X @ w 一行就完成了矩阵乘法，但你没法从这一行看出 "哪个下标在求和、哪个在遍历样本"。
2. 你能亲手感受到 for 循环有多慢。 当你看到 73 个样本 × 8000 轮要跑 0.7 秒， 而一个向量化实现只要 0.001 秒时，你会真正理解"为什么工程实现要向量化"。
3. 没有隐藏的数值细节。 广播、视图、dtype 提升这些 NumPy 的行为， 在你还不熟悉的时候会制造很多"结果对不上"的困惑。

换成 NumPy 是一个下午的事，反过来则很难。

画图脚本 make_figures.py 确实需要 matplotlib（因此也间接用到了 NumPy）， 但它和算法完全分离——你删掉它，教程的全部结论依然成立。

4.2 最小实现：30 行看懂全部

下面是一个完整可用的线性单元，包括训练和预测：

def dot(a, b):
    """向量点积 w·x"""
    return sum(ai * bi for ai, bi in zip(a, b))


class MinimalLinearUnit:
    def fit(self, X, y, lr=0.05, n_epochs=1000):
        m, n = len(X), len(X[0])
        self.w = [0.0] * n      # 权重初始化为 0
        self.b = 0.0            # 截距初始化为 0

        for epoch in range(n_epochs):
            gw = [0.0] * n      # 权重的梯度累加器
            gb = 0.0            # 截距的梯度累加器

            # ---- 遍历所有样本，累加梯度 ----
            for xi, yi in zip(X, y):
                r = yi - (dot(self.w, xi) + self.b)      # 残差 r = y - ŷ
                for j in range(n):
                    gw[j] -= r * xi[j]                   # -(y - ŷ) · x_j
                gb -= r                                  # -(y - ŷ)

            # ---- 取平均，然后更新参数 ----
            for j in range(n):
                self.w[j] -= lr * gw[j] / m
            self.b -= lr * gb / m
        return self

    def predict(self, X):
        return [dot(self.w, x) + self.b for x in X]

这 30 行就是全部。请对照第 3.3 节的两条梯度公式，逐行确认它们一一对应。

4.3 完整实现

code/linear_unit.py 里的 LinearUnit 类在这个基础上加了：

功能	参数 / 方法	说明
三种梯度下降	mode="batch" / "minibatch" / "sgd"	一批批扫数据
批大小	batch_size	minibatch 模式下的 $B$
L2 正则化	l2=λ	岭回归
动量	momentum=β	通常 0.9
学习率策略	lr_schedule="constant"/"inverse"/"step"	迭代中衰减
收敛判据	tol	梯度范数小于它就提前停
训练历史	model.history	每轮的损失、梯度范数、验证损失
参数轨迹	record_params=True	记录每一步的 (w, b)，用于画图
评估	model.score(X, y)	返回 $R^2$

同一个文件里还提供了配套工具：

• StandardScaler / MinMaxScaler：特征标准化，第 3.8 节的那个"救命稻草"。
• normal_equation()：正规方程解析解，用手写的高斯消元实现，用来给梯度下降对答案。
• numerical_gradient() / check_gradient()：有限差分梯度校验。
• condition_number()：用幂迭代 + 反幂迭代估计 $\kappa$， 正是第 3.8 节里那些数字的来源。
• PolynomialFeatures：多项式特征，用于第 5.4 节。
• train_test_split() / kfold_indices()：数据划分，支持按组划分。
• mse / rmse / mae / r2_score：评估指标。

4.4 标准化

标准化的公式是

$$x'_j = \frac{x_j - \mu_j}{\sigma_j}$$

其中 $\mu_j$ 和 $\sigma_j$ 是第 $j$ 列在训练集上的均值和标准差。

class StandardScaler:
    def fit(self, X):
        n = len(X[0])
        self.mean_ = [mean(col(X, j)) for j in range(n)]
        self.std_ = [std(col(X, j)) for j in range(n)]
        for j in range(n):                     # 常数列不能除以 0
            if self.std_[j] < 1e-12:
                self.std_[j] = 1.0
        return self

    def transform(self, X):
        return [[(row[j] - self.mean_[j]) / self.std_[j]
                 for j in range(len(row))] for row in X]

两个必须记住的纪律：

1. 均值和标准差只能从训练集算，然后应用到测试集。 用全部数据算统计量再划分，是数据泄漏的一种，会让测试误差偏乐观。
2. 标准化会破坏权重的物理量纲。 要解读权重，必须还原回去（见第 5.6 节）。

4.5 怎么知道梯度写对了

这是本文想强调的一个工程习惯，也是很多人忽略的一步。

梯度公式推错了、符号抄反了，代码依然能跑——它只是收敛得慢，或者收敛到错误的地方。 这种 bug 非常难通过观察损失曲线发现。

解决办法是有限差分。导数可以近似成

$$\frac{\partial f}{\partial \theta_j} \approx \frac{f(\theta_j + h) - f(\theta_j - h)}{2h}$$

这叫中心差分，截断误差是 $O(h^2)$，取 $h = 10^{-6}$ 时精度足够。

def numerical_gradient(f, theta, h=1e-6):
    """中心差分：(f(x+h) - f(x-h)) / 2h"""
    grad = [0.0] * len(theta)
    for i in range(len(theta)):
        orig = theta[i]
        theta[i] = orig + h
        f_plus = f(theta)
        theta[i] = orig - h
        f_minus = f(theta)
        theta[i] = orig
        grad[i] = (f_plus - f_minus) / (2.0 * h)
    return grad

然后在测试里比较解析梯度和数值梯度：

analytic, numeric, max_rel_error = check_gradient(X, y, theta)
assert max_rel_error < 1e-6

在本文的测试里，这个相对误差是 $7.8 \times 10^{-11}$。 如果你自己写了一个梯度下降，请一定先做这一步。 它会为你省下几个小时的调试时间。

4.6 三种梯度下降的对比

cd code
python3 demo_playground.py

这个脚本会把梯度下降的每一步都打印出来。下面是它的真实输出片段 （用 Hammett 的 $\sigma$ 常数拟合相对速率对数）：

    迭代          w          b         J(w,b)          dJ/dw          dJ/db
------------------------------------------------------------------------------
     0     0.0000     0.0000       0.108869       -0.22163       -0.04769
     1     0.1108     0.0238       0.084948       -0.19562       -0.01958
     2     0.2086     0.0336       0.066791       -0.17309       -0.00603
     3     0.2952     0.0367       0.052654       -0.15338        0.00031
     4     0.3719     0.0365       0.041557       -0.13602        0.00311
   ...
    50     0.9745     0.0102       0.000360       -0.00057        0.00003
   ...
   499     0.9770     0.0101       0.000360       -0.00000       -0.00000

跑 500 步之后：  w = 0.977017,  b = 0.010115
解析解（正规方程）：w = 0.977017,  b = 0.010115
两者之差：4.44e-16 / 4.34e-17  —— 梯度下降确实走到了同一个地方

两件事值得注意：

• dJ/dw 那一列逐渐趋近于 0。 "梯度接近 0"就是"到达碗底"的数学写法。
• 越靠近碗底，步子越短。 因为步长正比于梯度，而梯度本身在变小。 这是梯度下降最后会"慢下来"的原因，也是为什么需要 500 步才能收敛到 $10^{-16}$ 精度。

4.7 怎么算条件数

condition_number() 的实现值得单独说一句，因为它看起来比实际难。

$\kappa = \lambda_{\max}/\lambda_{\min}$，其中 $\lambda$ 是 $\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}/m$ 的特征值。算全部特征值很麻烦，但我们只需要最大的和最小的两个：

• 最大特征值用幂迭代：反复做 $\mathbf{v} \leftarrow \mathbf{A}\mathbf{v}/\|\mathbf{A}\mathbf{v}\|$， $\mathbf{v}$ 会收敛到主特征向量，$\mathbf{v}^\top\mathbf{A}\mathbf{v}$ 就是 $\lambda_{\max}$。
• 最小特征值用反幂迭代：每次解一次 $\mathbf{A}\mathbf{z} = \mathbf{v}$， 这样就等价于对 $\mathbf{A}^{-1}$ 做幂迭代，收敛到 $\lambda_{\min}$。 解线性方程组正好复用正规方程里的那个高斯消元。

需要注意：condition_number(X, fit_intercept=True) 默认会把一列 1 拼到 X 前面， 因为梯度下降真正感受到的条件数来自增广矩阵，而那一列 1 往往就是最大特征值的来源。

4.8 测试

cd code
python3 tests_linear_unit.py

41 个测试，分七组。其中最重要的几条：

测试	断言	为什么重要
梯度 vs 有限差分	相对误差 < $10^{-6}$	确保推导没抄错
无噪声数据最优点的梯度	恰好为 0	确保损失函数定义正确
梯度下降 vs 正规方程	参数一致到 4 位小数	确保优化过程正确
常数列的标准化	不产生除零	化学数据里经常有一整列相同
标准化降低条件数	降低超过 1000 倍	第 3.8 节的论断
按组划分无泄漏	组没有跨越训练/测试	第 5.3 节的论断
Arrhenius 活化能	与文献值相差 < 2%	化学结论可复现
烷烃外推	$[n, n^{2/3}]$ 比 $[n]$ 好 5 倍以上	第 5.4 节的论断

最后两条是**"化学结论的回归测试"**：如果有人不小心改坏了代码， 这些测试会立刻发现。这是一个很值得借鉴的习惯—— 把你论文里的关键数字写成断言。

4.9 如何运行

cd code

# 1. 先确认实现无误（41 个测试）
python3 tests_linear_unit.py

# 2. 把梯度下降的每一步打印出来
python3 demo_playground.py

# 3. 三个化学案例
python3 demo_arrhenius.py         # Arrhenius 与 Clausius-Clapeyron
python3 demo_bandgap.py           # 带隙预测（十个小节，约 30 秒）
python3 demo_alkane.py            # 烷烃沸点与外推

# 4. 重新生成全部插图（需要 matplotlib）
python3 make_figures.py

不需要安装任何东西，只要 Python ≥ 3.8。

---

5. 化学与材料中的应用

5.1 应用地图

在动手之前，先想清楚一个问题：什么时候该用线性单元？

适合的场景	不适合的场景
目标量是连续的（能量、带隙、溶解度、$T_g$）	目标是分类标签，且类间边界明显非线性
描述符和目标之间有单调趋势	存在强协同效应、阈值效应、活性悬崖
样本量小（几十到几百），特征数不多	样本量巨大（$>10^5$）且特征维度很高
你需要一个可解释的基线	你需要最高的预测精度
你需要判断"非线性有多重要"	你需要处理图 / 三维结构 / 序列

下面三个案例会把这张表具体化。它们的数字全部来自 code/，你可以自己复现。

5.2 案例一：Arrhenius 方程——一个教科书级的线性单元

科学问题：$\mathrm{N_2O_5}$ 的气相分解反应

$$2\,\mathrm{N_2O_5} \rightarrow 4\,\mathrm{NO_2} + \mathrm{O_2}$$

的活化能是多少？

数据（教科书里常见的五个温度点）：

$T$ / K	298.15	308.15	318.15	328.15	338.15
$k$ / s⁻¹	3.46×10⁻⁵	1.35×10⁻⁴	4.98×10⁻⁴	1.50×10⁻³	4.87×10⁻³

模型：Arrhenius 方程 $k = A e^{-E_a/RT}$。两边取对数：

$$\ln k = \ln A - \frac{E_a}{R}\cdot\frac{1}{T}$$

对照线性单元 $\hat{y} = w x + b$：

• $y = \ln k$
• $x = 1000/T$（注意那个 1000）
• $w = -E_a/R$
• $b = \ln A$

这就是一个一元线性单元，没有任何区别。 代码不需要改一行。

图 10　左：把 ln k 对 1000/T 作图，五个数据点几乎完美地落在一条直线上（R² = 0.9997），橙色的短竖线是残差。右：回到原始坐标看，同一份模型给出的是一条指数曲线。

结果：

量	梯度下降的结果	文献值	相对误差
斜率 $w$	−12.411	—	—
活化能 $E_a$	103.19 kJ/mol	103.6 kJ/mol	0.40%
指前因子 $A$	4.18×10¹³ s⁻¹	—	—
$R^2$	0.999736	—	—

五个数据点、一个线性单元、几行代码，把一个化学反应的能量壁垒算到了 1% 以内。

同一套代码，换到水的汽化焓。 Clausius-Clapeyron 方程

$$\ln P = C - \frac{\Delta H_{\text{vap}}}{R}\cdot\frac{1}{T}$$

和 Arrhenius 方程形式完全一样。把水的饱和蒸气压数据（0–100 °C）代进去， 梯度下降给出

$$\Delta H_{\text{vap}} = 43.27\ \text{kJ/mol}$$

而文献值：25 °C 时是 43.99 kJ/mol，100 °C 时是 40.65 kJ/mol。

拟合值正好落在两个端点之间——因为线性单元给出的是整个温度区间的平均值， 而真实的 $\Delta H_{\text{vap}}$ 本身随温度变化。 这不是模型算错了，而是模型假设错了：它假设 $\Delta H_{\text{vap}}$ 是常数。

这个案例的两个结论：

1. 化学里大量的"作图法"就是线性单元。你已经会的一半化学，可以直接翻译成一行 model.fit(X, y)。
2. 但这个案例真正的技术含量在特征上，不在算法上。 用 $1/T$ 而不是 $1000/T$ 画图，看起来只是乘了一个常数， 但它让条件数从 $5\times10^7$ 降到 $6\times10^3$—— 这是梯度下降能不能跑起来的区别。细节见第 3.8 节。

5.3 案例二：用线性单元预测二元化合物的带隙

科学问题：给定一个化合物的化学式，能不能预测它的带隙 $E_g$？

这是一个经典的"成分—性质"问题，也是材料信息学里最常被拿来练手的题目之一。 它的好性质在于：带隙是连续的，有明确的物理意义， 而且实验数据丰富。

数据：73 个化合物，包括单质（Si、Ge）、III-V（GaAs、GaN……）、 II-VI（ZnS、CdTe……）、碱金属卤化物（NaCl、KI……）、氧化物（TiO₂、SiO₂……）、 3d 过渡金属一氧化物（MnO、FeO、CoO、NiO、CuO）和硫族化物（PbS、MoS₂……）。 带隙为室温附近的文献常见值。

描述符：全部由元素周期表算出，不含任何实验测量值。

先写一个化学式解析器——它是真正的第一步。 把 "Al2O3" 这样的字符串翻译成 {"Al": 2, "O": 3}， 才能算出加权平均：

def parse_formula(formula):
    """parse_formula("Ca(OH)2") -> {"Ca": 1, "O": 2, "H": 2}"""
    ...

然后得到四个描述符：

描述符	定义	化学含义
$\Delta\chi$	电负性极差	化学键的离子性
$\bar{\chi}$	化学计量加权平均电负性	电子的平均"被拉"程度
$\bar{r}$	加权平均共价半径 /pm	原子大小、轨道重叠
$M$	摩尔质量 /(g/mol)	尺寸的另一个代理

顺带说一句：写这个解析器只花了 40 行，但它体现了化学数据处理的第一步—— 把"字符串形式的化学式"翻译成"可以被加权的数字向量"。 你在真实项目里会反复做这件事。

结果：

图 11　左：预测值 vs 实验值。虚线是完美预测线。四个描述符给出 R² = 0.788、RMSE = 1.51 eV。五个被标注的化合物是误差最大的。右：按化合物类别拆开看平均绝对误差——模型在哪里失灵一目了然。

$$R^2 = 0.788, \qquad \text{RMSE} = 1.51\ \text{eV}$$

怎么读这个结果？

• $R^2 = 0.79$ 意味着这个四维超平面解释了约 79% 的方差。 对只有四个描述符、且完全不含结构信息的模型来说，这已经不差了。
• 但它离"能用"还差得远：RMSE 1.51 eV，比整个可见光光子能量范围（1.6–3.1 eV）还大。
• 结论：它是个好基线，不是一个好模型。

特征越多越好吗？ 看消融实验：

描述符组合	$R^2$	RMSE / eV	条件数 $\kappa$
$\Delta\chi$	0.7227	1.730	16.5
$\Delta\chi + \bar{\chi}$	0.7536	1.631	381
$\Delta\chi + \bar{\chi} + \bar{r}$	0.7851	1.523	4.1×10⁶
四个全用	0.7884	1.511	2.3×10⁷

性能在涨，但边际递减（从 0.785 到 0.788 只涨了 0.003）， 条件数从 16 涨到 2300 万。 特征越多，问题越"扁"，梯度下降越难走。 这就是第 3.8 节讲的问题在真实数据上的样子。

然后是这份数据里最重要的一个发现。

同一份数据、同一个模型、同样的训练代码，只改数据划分方式：

做法 A：随机划分。 把 73 个化合物随机分成训练集和测试集，跑 6 个随机种子：

$$R^2 = 0.746 \quad (\text{RMSE } 1.54\ \text{eV})$$

做法 B：按体系划分。 整个碱金属卤化物族要么全在训练集里，要么全在测试集里：

$$R^2 = -0.124 \quad (\text{RMSE } 2.88\ \text{eV})$$

$R^2$ 从 0.75 掉到 −0.12。 而这两次实验之间，模型、特征、超参数、 连随机种子都一样，只改了数据怎么划分。

为什么差别这么大？

• 随机划分时，测试集里的 KI 和训练集里的 KCl、KBr 是"邻居"： 同族、同结构、描述符几乎一样。模型只要学会插值就够了。
• 按体系划分时，模型必须去预测它从没见过的化学族。 这才是真实科研里的情形：你要预测的是一个新体系。

掉的不是模型的能力，是评估的水分。

这是化学机器学习里最常见、也最致命的一个错误： 把同族样本随机分到训练集和测试集两侧，然后报告一个漂亮却虚假的数字。 详细讨论见参考文献 [19]、[20]。

残差里藏着化学。 按体系拆开平均绝对误差：

图 12　左：残差随电负性差的变化，红色三角是 3d 过渡金属一氧化物。右：残差 vs 预测值的标准诊断图。低带隙一侧被系统性高估，高带隙一侧被系统性低估，说明真实关系是弯的。

体系	样本数	平均 $\lvert$残差$\rvert$ / eV
单质（Si、Ge）	2	0.12
碱金属卤化物 I-VII	20	0.53
III-V	14	0.84
II-VI	14	1.05
硫族化物	6	1.22
氧化物	11	1.81
3d 过渡金属一氧化物	5	2.32

预测得最离谱的八个化合物：

化学式	实验 / eV	预测 / eV	残差 / eV	体系
SiO₂	9.00	4.30	+4.70	oxide
BeO	10.60	6.76	+3.84	II-VI
CuO	1.20	4.56	−3.36	TMO
Al₂O₃	8.80	5.55	+3.25	oxide
InN	0.64	3.87	−3.23	III-V
CdO	2.20	5.21	−3.01	II-VI
LiF	13.60	10.59	+3.01	I-VII
BN	6.40	3.67	+2.73	III-V

三件事值得停下来想一想：

1. 碱金属卤化物拟合得最好（平均误差 0.53 eV），因为它们的带隙几乎只由离子性决定， 而 $\Delta\chi$ 正好描述离子性。模型学到的规律是真的。
2. 3d 过渡金属一氧化物最差（平均误差 2.32 eV）。这是一个物理故事： 这些氧化物的带隙由 d 电子的强关联能（Hubbard $U$）决定，不是由电负性决定。 想让模型学会它，你需要加一个描述 d 电子占据数或 $U$ 的特征， 或者换成能表达电子结构的模型。
3. SiO₂ 被低估了 4.7 eV。 模型看到的是"电负性差不大、原子又小"， 于是预测成一个中等带隙的材料；但真实的 SiO₂ 是强共价网络固体，带隙 9 eV。 缺的那一维叫"结构"。

不同适用范围下的表现：

数据集	$n$	$R^2$	RMSE / eV
全部化合物	73	0.7884	1.511
去掉 3d 过渡金属氧化物	68	0.8305	1.383
只留 I-VII / II-VI / III-V	48	0.9007	1.072
只留碱金属卤化物	20	0.9259	0.554

这不是在挑数据。这是在回答一个更诚实的问题： "我的模型在什么范围内是可用的？" 科研里报一个"适用范围"比报一个全局 $R^2$ 有价值得多。

5.4 案例三：正构烷烃沸点——外推是怎么失效的

科学问题：正构烷烃 $\mathrm{C}_n\mathrm{H}_{2n+2}$ 的沸点随碳数怎么变？能不能预测 C20 的沸点？

这个案例的价值在于它演示的是一种特别危险的失败： 模型在训练集上表现完美，一到外推区间就错得离谱。

数据：C1–C20 的常压沸点，全部是手册里的真实数值。 前 10 个（C1–C10）用来训练，后 10 个（C11–C20）完全不给模型看。

$n$	1	2	3	4	5	…	10
$T_b$ / °C	−161.5	−88.6	−42.1	−0.5	36.1	…	174.1

$n$	11	12	…	20
$T_b$ / °C	195.9	216.3	…	342.7

图 13　同一个线性单元、同一份训练数据，只是特征不同。浅蓝色区域是训练区间（C1–C10），右侧是外推区间。红色虚线标出 C20 的真实沸点 342.7 °C。

第一幕：特征只有碳数 $n$。

这是最自然的想法。线性单元给出

$$T_b = 35.5\,n - 159$$

	训练 RMSE	外推 RMSE	$n=20$ 预测
只用 $n$	17.3 °C	129.1 °C	551 °C

训练集看着相当好（$R^2 = 0.97$），但 C20 预测 551 °C，真实值 342.7 °C——差 209 °C。

为什么？因为沸点随碳数的增长是"越来越慢"的。相邻烷烃的沸点差：

72.9, 46.5, 41.6, 36.6, 32.6, 29.7, 27.3, 25.1, 23.3   (°C)

从 72.9 一路掉到 23.3，而且还在继续变小。一条直线只会保持同一个斜率， 所以往外推得越远，错得越离谱。

第二幕：再加一个 $n^2$。

既然曲线是弯的，那就加一个平方项。（严格说这仍然是线性单元—— 因为对参数 $w$ 而言它依然是线性的，我们只是把特征从 $[n]$ 换成了 $[n, n^2]$。）

	训练 RMSE	外推 RMSE	$n=20$ 预测
$[n, n^2]$	5.6 °C	132.9 °C	95 °C

训练误差从 17.3 掉到 5.6，看起来是进步。 但外推预测 95 °C——比第一幕还糟，几乎回到了室温附近。

原因：二次抛物线在 C1–C10 上确实是向上弯的，但它必然有一个顶点。 拟合出来的顶点大概在 $n \approx 13$ 附近，过了顶点，抛物线就开始往下掉。 模型在训练区间里学到的是"上升"，在训练区间外给出的却是"下降"—— 这是纯粹的数学外推，没有任何化学在里面。

这是本教程最重要的一课：更灵活的特征会让训练误差变小，同时让外推误差变大。 只看训练集，你永远发现不了这件事。

第三幕：加一个 $n^{2/3}$，因为化学这么要求。

现在换一个思路：不要问"什么函数能拟合得更好"，而要问"什么物理量在控制沸点"。

答案：正构烷烃的沸点主要由色散力（范德华力）决定， 而色散能正比于分子之间接触的表面积。表面积怎么随 $n$ 变？做一次量纲分析：

• 对任何形状相近的凝聚相物体，表面积 $\propto$ 体积$^{2/3}$；
• 正构烷烃每加一个 CH₂，体积增加一个固定量，所以体积 $\propto n$；
• 于是接触表面积 $\propto n^{2/3}$。

这不是严格定理（真实分子是柔软的链，会卷曲），但它给出了一个可以检验的预言： 指数应该在 0.67 附近。我们把它跑出来：

指数 $p$（特征为 $[n, n^p]$）	0.2	0.4	0.5	0.6	0.7	0.8	1.0
外推 RMSE / °C	22.9	11.5	5.3	1.5	8.1	15.2	129.1

最优指数是 0.60，量纲分析预言的 0.67 就在它旁边，两者相差不到 0.1。

注意 $p = 1.0$ 那一行：外推误差反而暴涨到 129。因为 $[n, n^1]$ 就是 $[n, n]$—— 两列完全一样的特征。模型有 3 个参数却只有 2 个独立方向，权重变得任意。 这就是多重共线性。

三幕的对照：

特征	训练 RMSE	外推 RMSE	$n=20$ 预测
$[n]$	17.3 °C	129.1 °C	551 °C
$[n, n^2]$	5.6 °C	132.9 °C	95 °C
$[n, n^{2/3}]$	1.8 °C	5.8 °C	334 °C
（真值）	—	—	342.7 °C

三行数字，三个结论：

1. 训练误差不能用来选模型。 第二幕的训练误差比第一幕小了三倍，外推误差却更大。 任何一个只看训练集的评估方式都会选错。
2. 加特征不等于加信息。 $n^2$ 和 $n^{2/3}$ 都是"一个额外的数字"， 但它们携带的信息量天差地别。
3. 唯一能让外推变好的办法，是让特征里带上物理。 这不是调参能解决的事，它需要你知道一点关于分子间作用力的事—— 也就是需要你是一个化学家。

5.5 从权重里读出化学

线性单元的一个独特优势是：它的参数可以直接翻译成化学语言。 但前提是你得把权重还原到原始量纲。

在标准化空间里训练出来的权重，反映的是"每变化一个标准差的影响"， 比较大小很有用，但不能直接说"每单位"。 还原的方法是

$$w_j^{\text{raw}} = \frac{w_j^{\text{scaled}}}{\sigma_j}, \qquad b^{\text{raw}} = b^{\text{scaled}} - \sum_j w_j^{\text{raw}}\mu_j$$

带隙案例里还原后的权重：

描述符	标准化空间权重	原始量纲权重	读法
$\Delta\chi$	+2.93	+3.52 eV / 单位	电负性差越大，带隙越大
$\bar{\chi}$	−0.96	−2.62 eV / 单位	平均电负性越大，带隙越小
$\bar{r}$	−0.44	−0.017 eV / pm	原子越大，轨道重叠越好，带隙越小
$M$	−0.38	−0.005 eV / (g/mol)	与半径高度共线，几乎是冗余的
$b$	+4.45	+8.25 eV	外推到描述符全为 0 的截距

怎么读这些数：

• $\Delta\chi$ 的权重是正的，和化学直觉一致：离子晶体的价带和导带来自不同的原子， 能量差自然大。
• $\bar{r}$ 的权重是负的：原子越大，轨道重叠越好，能带越宽，带隙越小。
• $M$ 的权重很小，而且和半径高度共线。它基本上是一个"尺寸"的重复表达， 删掉它模型几乎不变。

但这里有一个必须说清楚的警告。

权重是相关，不是因果。$\bar{\chi}$ 的权重是负的， 你完全可以编一个故事来解释它（"电负性大意味着电子被束缚得紧，能带窄，所以带隙小"）—— 但这个故事和 $\Delta\chi$ 的故事在某些数据上是等价的。 在相关的描述符之间，权重可以随意互换。

这就是多重共线性。它也是为什么本教程反复强调：

要同时看权重、条件数、残差三样东西，不能只看权重。

一个更稳妥的做法是：先算条件数（第 3.8 节）， 如果 $\kappa > 30$ 就要警惕；然后用残差图找物理； 最后才去解释权重，并且明确说明哪些权重是不稳定的。

5.6 三个案例的交叉结论

把三个案例放在一起，会发现它们讲的是同一个故事。

案例	数据规模	算法	决定成败的东西
Arrhenius	5 个点	一元线性单元	特征变换（$1/T \to 1000/T$）和标准化
带隙	73 个化合物	四元线性单元	描述符设计和数据划分方式
烷烃沸点	20 个点	一元 / 二元线性单元	特征的物理动机

三个案例里，梯度下降本身从来没有出过问题—— 它每次都老老实实地收敛到了全局最优。

出问题的地方永远在别处：

• 特征的定义（$1/T$ 还是 $1000/T$？$n$ 还是 $n^{2/3}$？）
• 数据怎么划分（随机还是按体系？）
• 怎么评估（只看训练集还是看外推？）

这是本文最想让你记住的一句话：

在化学机器学习里，把 90% 的时间花在数据和特征上，把 10% 花在算法上， 通常是对的。

5.7 线性单元在科研流水线里的位置

图 14　线性单元在「AI × 化学」流水线里的位置。它不是一个终点，而是一把尺子：它的残差会告诉你下一步该往哪个方向加复杂度。

这张图想表达的是一个工作流，而不是一个模型：

1. 拿到数据（实验 / DFT / 数据库）。
2. 设计描述符（这一步决定了 80% 的结果）。
3. 标准化（几乎永远是必需的）。
4. 跑线性单元（本文的全部内容）。
5. 残差诊断 + 交叉验证（按体系划分，不要随机划分）。
6. 读出化学 / 决定是否上更重的模型。

第 6 步有两个出口：

• 如果线性单元已经给出 $R^2 > 0.9$，而且残差里没有结构， 那你大概不需要深度学习。把时间花在新的实验数据上，回报更大。
• 如果残差里有明显的结构（像本教程带隙案例里的过渡金属氧化物那样）， 那你已经知道该补什么了——这个信息比任何超参数搜索都有价值。

永远先跑线性单元。 这一步只要十分钟，但它给你一个基线， 还告诉你有多少信息藏在非线性里。

---

6. 优缺点：一份诚实的清单

6.1 优点

1. 它是凸问题，没有局部最优。 梯度下降从任何起点出发都会滑到同一个碗底。不需要调初始化，不需要担心鞍点。 这一点在化学里特别重要：你的结果不依赖于运气。

2. 它可以在几秒钟内训练完。 73 个化合物、4 个描述符、8000 轮，在纯 Python 里跑 0.76 秒； 换成 NumPy 是毫秒级。这意味着你可以做上百次交叉验证、 做完整的超参数扫描、做自助法（bootstrap）估计误差棒。

3. 参数可以直接翻译成化学。 每个权重就是"这个描述符每增加一个单位，性质变化多少"。 这是深度网络给不了的东西。本教程第 5.5 节演示了怎么读。

4. 它告诉你非线性的下限。 如果线性模型只有 $R^2 = 0.6$，你就有了一个有力的论据： 这个体系确实存在非线性效应。反过来，如果线性模型已经有 0.95， 上深度学习的边际收益可能远小于你的预期。

5. 它不需要大量数据。$n$ 个特征大约需要 $10n$ 个样本就能得到稳定的估计。 在化学里，这往往意味着你不需要再跑三个月 DFT。

6. 它没有超参数焦虑。 真正需要调的只有一个：学习率。而如果你用了标准化， 0.01–0.1 通常都能用。相比之下，一个 GNN 有十几个超参数。

7. 数值上有解析解可以对答案。 正规方程（第 3.9 节）给你一个独立的验证手段。 本文的实现里，梯度下降的结果和正规方程相差 $5.6\times10^{-9}$—— 这种"两条路走同一个答案"的验证，在别的模型上是没有的。

8. 它是最小的那个单元，因此完全可懂。 从输入到输出，中间没有任何黑箱。你可以逐行追踪每一个数字的来源。

6.2 缺点

1. 它只能表达线性关系。 这是定义决定的，无法绕过。如果你的体系里有协同效应、阈值效应、饱和效应， 线性单元在结构上就无法表达。

2. 它不会自动发现特征之间的相互作用。A 和 B 一起出现时的效果是单独效果的乘积——这种"协同项"$\Delta\chi_1 \cdot \Delta\chi_2$ 必须由你显式地构造出来。这不是模型的缺陷，是它的使用纪律。

3. 外推是危险的，而且它的外推错误是"直线式的"。 第 5.4 节的烷烃案例里，线性外推把 C20 的沸点预测成 551 °C， 真实值是 342.7 °C。线性模型在训练区间外会以惊人的自信给出错误的答案。

4. 对离群点极其敏感。 平方损失会让一个偏离很远的样本主导整个拟合。化学数据里的离群点 往往来自实验误差、晶型不同、或者一个真正有趣的异常体系—— 但线性单元不会区分这两种情况。解决办法是改用 Huber 损失或先做诊断。

5. 相关特征之间会互相"抢"权重。 多重共线性会让权重的符号变得不可信。带隙案例里 $M$ 和 $\bar{r}$ 的相关系数很高，它们的权重就变得几乎无法单独解读。

6. 它默认测量误差都在 $y$ 上。 如果描述符本身有不可忽略的误差（DFT 计算的描述符、实验测定的键长）， 普通最小二乘会有偏。应该考虑正交回归（total least squares）。

7. 它假设误差同方差。 如果你的误差随预测值变大（化学里非常常见，比如低浓度端的测量误差更大）， 应该用加权最小二乘（WLS），权重取 $1/\sigma_i^2$。

8. 它不能处理缺失值。 一个样本缺少任何一个描述符，就会让整行无法使用。 需要先做插补，而插补本身就是一门学问。

6.3 与其他方法的正面对比

方法	能表达非线性	需要调参	外推能力	可解释性	典型适用规模
线性单元 + 梯度下降	否	1 个（学习率）	以直线方式犯错	极高	小到中等
正规方程	否	0 个	同线性单元	极高	$n < 10^3$
岭回归（$L_2$）	否	1 个（$\lambda$）	同线性单元，更稳	高	共线特征多时
LASSO（$L_1$）	否	1 个（$\lambda$）	同线性单元	高（自带特征选择）	稀疏特征
决策树	是	3–5 个	很差（无法外推）	中	小数据
随机森林 / 梯度提升	是	5–10 个	很差（无法外推）	中低	中等数据
高斯过程回归	是	2–3 个	中等（会回归均值）	中（有不确定性）	小数据（$<10^4$）
核岭回归 / SVR	是	2–4 个	中等	低	小到中等
神经网络	是（最强）	10+ 个	差（不可控）	极低	大数据
图神经网络	是（用结构信息）	10+ 个	差	低	大数据 + 结构

这张表里有一个非常值得注意的细节：树模型不能外推。

决策树和随机森林的预测值永远是训练集里出现过的那些值的组合。 你用 C1–C20 的烷烃数据训练一个随机森林，它会给出 C20 附近的值作为 C30 的预测—— 它连"外推"这个动作都做不了。

所以：如果你的应用需要外推（化学里非常常见：预测一个还没合成的化合物）， 线性模型和核方法往往比树模型更合适，即使树模型在内部测试上表现更好。

6.4 十二个最常见的坑

这一节是本教程最有实用价值的部分。每一条都对应一个真实发生过的错误。

1. 不做标准化就直接梯度下降。 症状：损失爆炸成 inf，或者跑了几万轮还停在原地。 诊断：算一下条件数（第 4.7 节）。超过 1000 就一定要处理。 解法：StandardScaler，或者换一个量纲合适的特征。

2. 随机划分数据集，把同族样本分到两边。 症状：测试 $R^2$ 很漂亮，一预测新体系就崩。 诊断：检查训练集和测试集里有没有同一个化学族的成员。 解法：按骨架、按元素组成、按空间群、按论文来源分组划分。 这是化学机器学习里最常见的一个致命错误。 本教程第 5.3 节给了实测数字： 随机划分 $R^2 = 0.75$，按体系划分 $R^2 = -0.12$。

3. 用测试集来选特征。 症状：换了十组特征，选了在测试集上最好的那组；报告的测试误差偏乐观。 解法：用训练集内部的交叉验证来选特征，测试集只在最后用一次。

4. 拿模型的输出直接外推。 症状：烷烃案例里把 C20 预测成 551 °C。 解法：永远报告适用范围。"本模型适用于 $E_g < 5$ eV 的离子性化合物"， 这样一句话在论文里比 $R^2$ 更值钱。

5. 忽略多重共线性。 症状：权重的符号和化学直觉相反，或者加一个样本权重就翻转。 诊断：算条件数；算特征之间的相关系数矩阵。 解法：删掉冗余特征、改用岭回归、或者用主成分分析。 不要相信共线特征之间的权重分配。

6. 单位不统一。 症状：模型给的系数看起来很奇怪。 检查：能量是 kJ 还是 kcal（差 4.184 倍）？ 长度是 Å 还是 pm（差 100 倍）？温度是 K 还是 °C？ 在化学里这是一个高发错误。 本教程的所有数据都标了单位，请你也这样做。

7. 只报 $R^2$，不报带单位的误差。$R^2$ 无量纲，看起来很"科学"，但它无法告诉你这个模型能不能用。 化学家需要的是"预测误差是 1.5 eV"这种能和自己实验精度比较的数字。 本文的所有案例都同时报了 $R^2$ 和 RMSE。

8. 让离群点主导拟合。 症状：模型在 95% 的样本上很好，在 5% 上错得离谱， 而这 5% 正是你最关心的新体系。 诊断：画残差图（第 5.3 节的图 12 就是标准做法）。 解法：改用 Huber 损失、或者先做异常值分析。 但注意：不要无脑删掉离群点。 有时候它们才是发现新物理的入口。

9. 数据泄漏的隐蔽形式：用了"未来信息"。 举例：预测一个反应的产率，特征里放了产物的 NMR 数据； 预测材料的稳定性，特征里放了它已经被合成出来这个事实。 这类泄漏最难发现，因为它们看起来像是正常的描述符。 自查方法：想象你在真正做预测的那一刻，这个特征能不能拿到？

10. 目标变量该取对数却没取。 化学里很多量是正数、跨越几个数量级、而且相对误差比绝对误差更稳定： 速率常数、溶解度、分配系数、电导率。 这类量通常应该先取 $\log$ 再回归（本教程的 Arrhenius 案例就是）。 诊断：画残差 vs 预测值，如果呈现"喇叭形"，说明该做变换了。

11. 样本太少就上复杂模型。 经验法则：$n$ 个特征大约需要 $10n$ 个样本才能得到稳定估计。 你有 30 个样本、200 个描述符的时候，线性单元必须配合正则化， 而深度网络在这个规模上毫无意义。

12. 不做残差诊断就下结论。 残差图是化学机器学习里信息密度最高的一张图。 它会告诉你：模型在哪里失灵、失灵的样本有没有共同点、 真实关系是不是弯的、还缺哪一维物理。 本教程带隙案例里关于"过渡金属氧化物的 d 电子强关联"这个结论， 完全是从残差图里读出来的。

---

7. 接下来学什么

如果你跟着这份教程把代码跑通了，你现在站在一条很清楚的路口上。 下面五个方向，按"离你最近"排序。

第一步：把线性单元换一种输出。 把最后的恒等函数换成 sigmoid，你就得到逻辑回归—— 它做的是分类，但更新规则和本文一模一样（这是因为交叉熵损失配 sigmoid 时， 梯度恰好又是"残差 × 输入"）。学会逻辑回归，你就同时学会了 神经网络输出层、以及所有二分类问题的基线。

第二步：学会正则化和特征选择。 把 $L_2$ 换成 $L_1$（LASSO），你会发现一部分权重被精确地压成了 0—— 模型自动做了特征选择。这在化学里极其有用， 因为"200 个描述符里只有 8 个真正重要"是一个常见的情形。

第三步：处理非线性。 三条路线：

• 特征工程：本教程第 5.4 节的做法。用化学知识构造非线性特征 （$n^{2/3}$、$1/T$、交乘项）。可解释性最好，但需要你的化学直觉。
• 核方法：核岭回归、支持向量回归。把数据隐式映射到高维空间， 在那里找线性关系。数学优美，小数据上很强。
• 神经网络：让模型自己学非线性。表达力最强，但需要数据和调参， 而且可解释性会大幅下降。

第四步：处理结构信息。 分子不是一组数字，它有图结构。图神经网络（GNN） 把分子当作原子+化学键的图， 用消息传递来学习表示。这是今天分子性质预测的主流方法。 好消息是：GNN 的最后一层，仍然是一个线性单元。

第五步：学会问对问题。 真正的科研不是"给定数据预测性质"，而是"在有限的计算/实验资源下， 下一个该做哪个？"这属于主动学习（active learning）和 贝叶斯优化的范畴。它们和线性单元的关系是： 你需要一个能给出不确定性的模型，而线性单元的贝叶斯版本 （贝叶斯线性回归）恰好是最简单的起点。

一个具体的建议： 不要急着学 Transformer。 先把本文的十二个坑逐一在自己手头的数据上检验一遍。 你会发现问题里 80% 的困惑，在第 6.4 节里都有。

---

8. 习题

前八题用 code/ 里的代码就能做，后四题需要你自己写一点代码。

第一部分：读懂代码

题 1. 把 linear_unit.py 里 LinearUnit.gradient() 的四行核心代码 抄下来，逐行标注它对应第 3.3 节的哪一个数学符号。

题 2. 运行 python3 demo_playground.py，把第一部分表格里 第 0 步和第 1 步的 dJ/dw 手算一遍，验证代码算对了。 （提示：$\sigma$ 的 13 个值你都拿得到。）

题 3. 把 LinearUnit.__init__ 里的 lr 从 0.5 改成 1.0 和 2.5， 分别跑 demo_playground.py 的第一部分，观察损失曲线。 用 $\eta < 2/\lambda_{\max}$ 解释你看到的现象。

题 4. 删掉 StandardScaler，在 demo_bandgap.py 第 4 节里 找一个能让不缩放的数据收敛的学习率。它比标准化的那一个大还是小？差几个数量级？

第二部分：改代码

题 5. 给 LinearUnit 加一个 Huber 损失选项：

$$\ell(r) = \begin{cases} \frac{1}{2}r^2, & |r| \le \delta \\ \delta(|r| - \frac{1}{2}\delta), & |r| > \delta \end{cases}$$

推导出它的梯度，然后用有限差分验证。 最后在带隙数据里人为把一个样本的带隙改成 50 eV，比较 MSE 和 Huber 的结果。

题 6. 实现加权最小二乘：给每个样本一个权重 $s_i > 0$， 损失变成 $J = \frac{1}{2m}\sum_i s_i r_i^2$。 写出梯度，实现它，然后验证 $s_i$ 全部相等时结果和原来一致。

题 7. 实现 RidgeRegression 的正规方程解： $\boldsymbol\theta^* = (\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}} + \lambda \mathbf{I}')^{-1}\tilde{\mathbf{X}}^\top\mathbf{y}$， 其中 $\mathbf{I}'$ 是单位矩阵但截距对应的那一项为 0（截距不做正则化）。 验证它和梯度下降 + l2 给出的结果一致。

题 8. 给 condition_number() 换成 Lanczos 迭代，比较两者的结果和收敛速度。

第三部分：做点化学

题 9. 在 chem_data.py 里加入 20 个新的化合物 （可以从你手边的物理化学手册里找），重跑 demo_bandgap.py。 新数据和原数据上的误差分布一样吗？如果不一样，为什么？

题 10. 用 build_water_vapor_dataset() 的数据， 在 0–50 °C 和 50–100 °C 上分别拟合 $\Delta H_{\text{vap}}$。 两个值差多少？这说明了线性模型的什么假设在失效？

题 11. 烷烃案例里，把目标从 $T_b$ 换成 $\log T_b$（$T_b$ 用 K）， 特征只用 $\log n$。这等价于假设 $T_b = a n^b$。 拟合出来的 $b$ 是多少？它和外推效果最好的指数 0.6 有什么关系？

题 12. 设计一个新特征，让带隙模型的 $R^2$ 提高 0.02 以上。 提示：想想模型失败的那三类化合物有什么共同点 （过渡金属氧化物、SiO₂、Al₂O₃）。可以用的信息有： 是否有未满的 d 壳层、周期表中第几周期、阳离子的氧化态。 注意：不要用带隙本身来构造特征——那是数据泄漏。

---

9. 参考文献与延伸阅读

原始文献

1. Widrow, B., & Hoff, M. E. Adaptive Switching Circuits. IRE WESCON Convention Record, Part 4, 96–104, 1960. —— Adaline 与 LMS 规则的原始论文。本文第 1.2 节的历史即基于此文。
2. Rosenblatt, F. The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain. Psychological Review, 65(6), 386–408, 1958. —— 感知器的原始论文，用来对照"阶跃"与"恒等"的区别。
3. Hammett, L. P. The Effect of Structure upon the Reactions of Organic Compounds. Benzene Derivatives. Journal of the American Chemical Society, 59(1), 96–103, 1937. —— 线性自由能关系的起点， 比感知器早 21 年。
4. Hansch, C., Maloney, P. P., Fujita, T., & Muir, R. M. Correlation of Biological Activity of Phenoxyacetic Acids with Hammett Substituent Constants and Partition Coefficients. Nature, 194, 178–180, 1962. —— Hansch 分析，QSAR 的奠基工作。
5. Rumelhart, D. E., Hinton, G. E., & Williams, R. J. Learning Representations by Back-Propagating Errors. Nature, 323, 533–536, 1986. —— 反向传播。你会发现它和本文的梯度推导是同一个套路。

教材

6. Bishop, C. M. Pattern Recognition and Machine Learning. Springer, 2006. —— 第 3 章讲线性回归，从最大似然一路讲到贝叶斯线性回归。数学严谨。
7. Hastie, T., Tibshirani, R., & Friedman, J. The Elements of Statistical Learning. Springer, 2nd ed., 2009. —— 第 3 章讲最小二乘与岭回归，第 7 章讲模型选择。 免费电子版可在作者主页下载。
8. Goodfellow, I., Bengio, Y., & Courville, A. Deep Learning. MIT Press, 2016. —— 第 4 章讲数值计算，第 8 章讲优化。 有中文版，第 4 章讲条件数的那一节强烈推荐。
9. Nocedal, J., & Wright, S. J. Numerical Optimization. Springer, 2nd ed., 2006. —— 如果你想把"学习率、条件数、收敛速度"彻底搞懂，这一本是标准参考。
10. Leach, A. R., & Gillet, V. J. An Introduction to Chemoinformatics. Springer, 2007. —— QSAR / QSPR 的经典教材， 讲描述符、模型选择、以及怎么把结果写进论文。

化学信息学与材料信息学

11. Butler, K. T., Davies, D. W., Cartwright, H., Isayev, O., & Walsh, A. Machine Learning for Molecular and Materials Science. Nature, 559, 547–555, 2018. —— 入门必读的综述。
12. Schmidt, J., Marques, M. R. G., Botti, S., & Marques, M. A. L. Recent Advances and Applications of Machine Learning in Solid-State Materials Science. npj Computational Materials, 5, 83, 2019. —— 材料 ML 的综述，里面有大量"描述符—性质"的线性基线案例。
13. Ramakrishnan, R., Dral, P. O., Rupp, M., & von Lilienfeld, O. A. Quantum Chemistry Structures and Properties of 134 Kilo Molecules (QM9). Scientific Data, 1, 140022, 2014. —— 最常用的分子数据集之一， 适合拿来练手。

关于方法学的警告文献（强烈推荐）

14. Wallach, I., & Heifets, A. Most Ligand-Based Classification Benchmarks Reward Overfitting and Loss of Chemical Validity, with Lessons for the Coming Era of Deep Learning. Journal of Chemical Information and Modeling, 58(5), 916–932, 2018.
15. Kapoor, S., & Narayanan, A. Leakage and the Reproducibility Crisis in Machine-Learning-Based Science. Patterns, 4(9), 100804, 2023. —— 关于数据泄漏的必读文献，第 5.3 节和第 6.4 节的核心参考。
16. Chuang, K. V., & Keiser, M. J. Comment on "Predicting the Solubility of Organic Compounds by Graph Convolutional Networks with Multitask Learning". Journal of Cheminformatics, 10, 61, 2018. —— 一个关于"改进的划分方式会让模型看起来变差"的经典案例。

数据与工具

17. RDKit: https://www.rdkit.org/ —— 分子描述符与指纹的事实标准。
18. pymatgen: https://pymatgen.org/ —— 材料结构处理。
19. matminer: https://hackingmaterials.lbl.gov/matminer/ —— 材料描述符库。
20. Materials Project: https://materialsproject.org/ —— 材料性质数据库。
21. CRC Handbook of Chemistry and Physics —— 本文烷烃沸点与水的蒸气压数据来源。

---

附录 A：三个数学补充

A.1 那个 $\frac{1}{2}$ 从哪来

损失函数

$$J = \frac{1}{2m}\sum_i (y_i - \hat{y}_i)^2$$

里的 $\frac{1}{2}$ 纯粹是记号上的方便。求导时：

$$\frac{\partial}{\partial \hat{y}_i}\left[\frac{1}{2m}(y_i - \hat{y}_i)^2\right] = \frac{1}{2m}\cdot 2(y_i - \hat{y}_i)\cdot(-1) = -\frac{y_i - \hat{y}_i}{m}$$

如果分母写成 $m$ 而不是 $2m$，结果会多出一个因子 2， 梯度就是 $-\frac{2(y_i-\hat{y}_i)}{m}$，于是更新时学习的有效步长变成两倍。

这不影响最优解的位置——把目标函数整体乘一个正常数， 最小值点不变。它只影响梯度的大小，也就是影响"学习率取多少合适"。

在文献里你会看到两种约定都在用。只要你前后一致，就没有问题。 本文统一用 $\frac{1}{2m}$。

A.2 凸性的完整证明

定义：$f: \mathbb{R}^n \to \mathbb{R}$ 是凸函数，当且仅当对任意 $\mathbf{u}, \mathbf{v}$ 和 $t \in [0, 1]$，

$$f(t\mathbf{u} + (1-t)\mathbf{v}) \le t f(\mathbf{u}) + (1-t) f(\mathbf{v})$$

定理：如果 $f$ 二阶连续可微且 Hessian 矩阵 $\mathbf{H}(\mathbf{x}) \succeq 0$ （半正定）处处成立，那么 $f$ 是凸函数。

证明线性单元的损失满足这个条件。

把参数写成 $\boldsymbol\theta = (\mathbf{w}, b)$，输入写成增广向量 $\tilde{\mathbf{x}}_i = (\mathbf{x}_i, 1)$。损失是

$$J(\boldsymbol\theta) = \frac{1}{2m}\sum_{i=1}^m \bigl(y_i - \boldsymbol\theta^\top \tilde{\mathbf{x}}_i\bigr)^2$$

这是一个"仿射函数的平方和"。对 $\boldsymbol\theta$ 求梯度：

$$\nabla J = -\frac{1}{m}\sum_i \bigl(y_i - \boldsymbol\theta^\top\tilde{\mathbf{x}}_i\bigr)\tilde{\mathbf{x}}_i$$

再求一次导数（注意 $\tilde{\mathbf{x}}_i$ 与 $\boldsymbol\theta$ 无关）：

$$\mathbf{H} = \frac{\partial^2 J}{\partial \boldsymbol\theta\partial\boldsymbol\theta^\top} = \frac{1}{m}\sum_{i=1}^m \tilde{\mathbf{x}}_i\tilde{\mathbf{x}}_i^\top$$

现在验证它是半正定的。对任意 $\mathbf{v} \neq \mathbf{0}$：

$$\mathbf{v}^\top \mathbf{H}\mathbf{v} = \frac{1}{m}\sum_{i=1}^m \mathbf{v}^\top\tilde{\mathbf{x}}_i\tilde{\mathbf{x}}_i^\top\mathbf{v} = \frac{1}{m}\sum_{i=1}^m \bigl(\mathbf{v}^\top\tilde{\mathbf{x}}_i\bigr)^2 \ge 0$$

因为它是若干个实数的平方和。$\blacksquare$

这个证明告诉我们两件事：

1. 凸性不是假设，是结构决定的。 只要损失是"残差的平方和"，它就自动成立。
2. $\mathbf{H}$ 是奇异的当且仅当存在非零 $\mathbf{v}$ 使得 $\mathbf{v}^\top\tilde{\mathbf{x}}_i = 0$ 对所有 $i$ 成立。 也就是说，存在一个参数方向的改变完全不影响任何预测值—— 这正是"特征共线"和"样本数少于特征数"的情形。 这时最优解不唯一（存在一条直线上的最优解）， 但梯度下降依然会收敛到其中一个（具体是哪一个取决于初始化和轨迹）。

A.3 条件数与收敛速度

这一节把第 3.8 节的结论严格化。

考虑二次型目标

$$f(\boldsymbol\theta) = \frac{1}{2}\boldsymbol\theta^\top \mathbf{A}\boldsymbol\theta - \mathbf{b}^\top\boldsymbol\theta, \qquad \mathbf{A} \succ 0$$

它的最优解满足 $\mathbf{A}\boldsymbol\theta^* = \mathbf{b}$。梯度下降是

$$\boldsymbol\theta_{t+1} = \boldsymbol\theta_t - \eta(\mathbf{A}\boldsymbol\theta_t - \mathbf{b})$$

定义误差 $\mathbf{e}_t = \boldsymbol\theta_t - \boldsymbol\theta^*$，代入得

$$\mathbf{e}_{t+1} = (\mathbf{I} - \eta\mathbf{A})\mathbf{e}_t$$

$\mathbf{A}$ 是对称正定的，可以做特征分解 $\mathbf{A} = \mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^\top$， 特征值 $0 < \lambda_1 \le \lambda_2 \le \cdots \le \lambda_n$。 把误差投影到特征基上，每个分量独立地衰减：

$$e_t^{(j)} = (1 - \eta\lambda_j)^t\, e_0^{(j)}$$

所以收敛速度由最慢的那个分量决定，即 $\max_j |1 - \eta\lambda_j|$。 这个最大值由 $\lambda_1$ 和 $\lambda_n$ 端点的两个分量主导。

最优学习率 就是让两端的衰减因子相等的那个值：

$$\eta_{\text{opt}} = \frac{2}{\lambda_1 + \lambda_n}$$

代入端点：

$$1 - \eta_{\text{opt}}\lambda_1 = \frac{\lambda_n - \lambda_1}{\lambda_n + \lambda_1} = \frac{\kappa - 1}{\kappa + 1}$$

$$\left|1 - \eta_{\text{opt}}\lambda_n\right| = \frac{\lambda_n - \lambda_1}{\lambda_n + \lambda_1} = \frac{\kappa - 1}{\kappa + 1}$$

（中间的特征值给出的因子绝对值更小，所以不影响最坏情况。）于是

$$\boxed{\;\|\mathbf{e}_t\| \le \left(\frac{\kappa-1}{\kappa+1}\right)^t \|\mathbf{e}_0\|, \qquad \kappa = \frac{\lambda_n}{\lambda_1}\;}$$

要达到精度 $\varepsilon$，需要

$$t \ge \frac{\ln(1/\varepsilon)}{\ln\left(\frac{\kappa+1}{\kappa-1}\right)} \;\approx\; \frac{\kappa}{2}\ln\frac{1}{\varepsilon} \qquad (\kappa \gg 1)$$

这就是"条件数翻多少倍，迭代次数就翻多少倍"的严格版本。

代入本文的数字：带隙数据不缩放时 $\kappa \approx 2.3\times10^7$， 标准化后 $\kappa \approx 41$。两者相差 $5.5\times10^5$ 倍—— 这就是为什么一个跑不动、一个几步就收敛。

最后一点有用的推论。 上面假设 $\mathbf{A}$ 的特征值是 $\lambda$， 而在线性回归里 $\mathbf{H} = \frac{1}{m}\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}}$， 所以

$$\lambda_{\max}(\mathbf{H}) \le \operatorname{tr}(\mathbf{H}) = \frac{1}{m}\sum_{i=1}^m\|\tilde{\mathbf{x}}_i\|^2$$

标准化之后，每一列的均方都是 1，所以迹大约等于 $n+1$， 于是 $\lambda_{\max} \lesssim n+1$，稳定性条件 $\eta < 2/\lambda_{\max}$ 就给出 $\eta < 2/(n+1)$。

对本文的 2 特征数据集，这是 $\eta < 0.67$； 而实测的发散点出现在 $\eta \approx 1.65$。 理论上界是保守的（因为它用的是迹估计，不是真实的 $\lambda_{\max}$）， 但它的量级是正确的，而且它解释了为什么标准化之后学习率可以放心地取 0.01–0.1。

---

附录 B：文件结构

线性单元与梯度下降教程/
├── README.md                     快速上手与结论速览
├── 线性单元与梯度下降教程.md       Markdown 源文件（含 LaTeX 公式）
├── 线性单元与梯度下降教程.html     网页版（公式离线渲染，双击即可打开）
├── assets/
│   ├── tutorial.css              网页版样式
│   └── katex/                    KaTeX 公式渲染器（离线，勿删）
├── figures/
│   ├── fig01_perceptron_vs_linear_unit.png    图 1   感知器 vs 线性单元
│   ├── fig02_geometry.png                     图 2   残差与损失曲面
│   ├── fig03_loss_surface.png                 图 3   损失曲面
│   ├── fig04_gradient_path.png                图 4   负梯度场与轨迹
│   ├── fig05_learning_rate.png                图 5   学习率的四种命运
│   ├── fig06_lr_real_data.png                 图 6   真实数据上的学习率
│   ├── fig07_feature_scaling.png              图 8   缩放前后的等高线
│   ├── fig08_batch_vs_sgd.png                 图 7   批量 vs 随机
│   ├── fig09_arrhenius.png                    图 10  Arrhenius 拟合
│   ├── fig10_arrhenius_scaling.png            图 9   条件数实验
│   ├── fig11_bandgap_pred_vs_true.png         图 11  带隙预测
│   ├── fig12_bandgap_residuals.png            图 12  残差诊断
│   ├── fig13_alkane_extrapolation.png         图 13  烷烃外推
│   ├── fig14_pipeline.png                     图 14  科研流水线
│   └── （每个图都有同名 .svg 矢量版本）
└── code/
    ├── linear_unit.py           核心库（约 700 行，零依赖）
    │   ├── LinearUnit           线性单元 + 批量/小批量/随机梯度下降
    │   ├── StandardScaler       标准化
    │   ├── MinMaxScaler         [0, 1] 归一化
    │   ├── PolynomialFeatures   多项式特征
    │   ├── normal_equation      正规方程（手写高斯消元）
    │   ├── condition_number     幂迭代估计条件数
    │   ├── check_gradient       有限差分梯度校验
    │   └── train_test_split     支持按组划分
    ├── chem_data.py             化学数据集（元素表 + 化学式解析器 + 四组真实数据）
    ├── demo_playground.py       把梯度下降每一步打印出来
    ├── demo_arrhenius.py        案例一：Arrhenius 与 Clausius-Clapeyron
    ├── demo_bandgap.py          案例二：带隙预测（十节）
    ├── demo_alkane.py           案例三：烷烃沸点与外推（四幕）
    ├── tests_linear_unit.py     41 个自动化测试
    └── make_figures.py          生成全部插图（需要 matplotlib）

注：figures/ 里的文件编号是"生成顺序"，正文里的图号是"出现顺序"， 两者不一致。正文中的每一张图都标注了对应的文件名，查找时以上表为准。