图宾根机器学习数值方法笔记-全-

图宾根机器学习数值方法笔记（全）

001：引言 🎯

在本节课中，我们将要学习数值算法在机器学习中的核心作用。我们将探讨为什么理解这些“引擎”的内部工作原理对于构建高效、可靠的机器学习系统至关重要，并概述本课程将涵盖的四大主题：线性代数、模拟、积分和优化方法。

概述：为什么需要关心数值算法？

机器学习算法是处理数据并生成可用于预测或行动的模型的计算机程序。输入是数据，输出是模型。而在这两者之间发生的，是计算机执行的计算过程。

与基于规则的经典人工智能不同，当代机器学习算法使用数值算法作为这些计算的基本单元。经典算法（如最近邻搜索）总是能给出正确答案。而数值算法是可能出错的方法，它们用于估计那些没有封闭形式解（或解在合理时间内不可计算）的数学量。

数值算法与原子操作的关键区别在于：原子操作（如指数函数 exp(x)）是“不会出错”的基础函数。而数值算法（如线性代数求解器）是“可能以微妙方式出错”的复杂方法，需要使用者确保其正确运行。

由于这些算法是机器学习的核心，作为机器学习从业者，你必须理解并思考它们。

数值算法的特点与挑战

以下是数值算法的一些典型特征：

  • 历史悠久：许多核心算法（如龙格-库塔法）发明于一百多年前（1905年）。
  • 接口复杂：调用时需要提供大量参数，必须清楚自己在做什么。
  • 结构精密：算法内部包含许多“魔法数字”，这些数字经过精心推导，任意改动都可能导致算法失效。
  • 黑盒化：通常以编译库（如Fortran代码）的形式提供，源代码难以查看和修改。

机器学习领域比这些数值方法年轻得多。当机器学习工程师遇到需要数值计算的任务时，他们发现早已存在由数学家、物理学家等领域专家编写的成熟方法。直接使用这些“古董”代码很方便，但也带来了问题：如果我们的任务与算法预设的任务不完全匹配，我们将束手无策。

本课程核心内容预览

接下来，我们将通过预览本课程涵盖的四大类数值算法，来理解为什么需要以新的视角看待它们。

1. 线性代数方法 🧮

上一节我们介绍了数值算法的普遍特性，本节中我们来看看线性代数在机器学习中的基础作用。

线性代数出现在几乎所有机器学习算法的底层。以高斯过程回归为例，其核心计算涉及求解形如 (K + σ²I)⁻¹y 的线性系统，其中 K 是核矩阵。直接计算矩阵逆 np.linalg.inv() 是错误做法。正确的做法是使用乔列斯基分解等数值稳定的方法。

# 正确做法：使用乔列斯基分解求解
from scipy.linalg import cho_factor, cho_solve
L, lower = cho_factor(K + sigma_sq * np.eye(N))
alpha = cho_solve((L, lower), y)

然而，scipy.linalg.cho_factor 底层调用的是BLAS/LAPACK库中的Fortran代码，对我们而言是一个黑盒。这导致了一个严重问题：对于 n 个数据点，乔列斯基分解的计算复杂度是 O(n³)。对于百万级的数据集，计算将变得完全不可行。因此，人们常认为核方法无法扩展到大数据集。

但在本课程中，Marvin 和 Jonathan 将讲授，如果了解数据的结构（如稀疏性、低秩性），或者接受近似解，我们可以使用更高效的迭代求解器（如共轭梯度法），将复杂度从 O(n³) 显著降低，甚至达到接近线性的复杂度。理解算法内部原理，才能突破“核方法太慢”的固有认知。

2. 模拟方法 🎮

上一节我们探讨了线性代数中的效率问题，本节中我们转向模拟方法，即求解动力系统轨迹的算法。

模拟方法解决的是微分方程，形式如 D[x](t) = f(x(t), t)，其中 D 是微分算子。这包括常微分方程（ODE）、偏微分方程（PDE）等。在机器学习中，模拟方法出现在两个领域：1）需要预测未来的智能体（如机器人、强化学习）；2）将物理知识编码为微分方程的科学机器学习。

一个典型例子是预测COVID-19病例。单纯用深度神经网络或高斯过程进行曲线拟合外推效果很差。我们需要利用机制性知识，例如SIR传染病模型：

dS/dt = -β * I * S / N
dI/dt = β * I * S / N - γ * I
dR/dt = γ * I

这里 β（接触率）和 γ（恢复率）是未知参数。传统ODE求解器（如龙格-库塔法）是一个初始值问题，要求方程完全已知。它无法直接利用观测数据来实时推断像 β(t) 这样的时变参数。

尽管现代框架（如JAX）重新实现了这些ODE求解器，并支持自动微分，但其核心算法（及内部的“魔法数字”）依然来自上世纪的Fortran代码。为了同时拟合数据和求解方程，我们不得不采用“模拟推断”的方法，在优化参数的外循环中，反复调用整个ODE求解器（内循环），计算效率低下。

Younès 和 Nathanaël 将讲授一种新思路：将模拟过程构建为一个马尔可夫链，其状态是概率分布。我们可以向这个链中插入不同的信息算子，例如“ODE成立”的算子，或“轨迹经过观测数据点”的算子。这样，我们可以在单次前向计算中，同时利用物理机制和观测数据，非参数地推断时变参数 β(t) 并预测未来。这种基于概率推断的框架比古老的、僵化的黑盒算法更灵活、更强大。

3. 积分方法 📊

上一节我们看到了如何将模拟重构为概率推断，本节中我们简要看看积分方法。

积分是概率推断中的基本操作，例如计算边缘分布 p(y) = ∫ p(y, z) dz。标准的做法是马尔可夫链蒙特卡洛（MCMC）方法。可以说，MCMC是“最普适”也是“在多数情况下效率最低”的算法。

在本课程中，我们将探讨在低维或具有特殊结构的问题中，存在更高效、性质更好的积分算法（如高斯积分、贝叶斯求积法）。理解这些算法有助于我们认识经典数值方法与基于概率的新计算范式之间的联系。

4. 优化方法 ⚙️

上一节我们提到了积分，本节我们来到机器学习最核心的数值操作：优化。

现代机器学习的基础模板是经验风险最小化：min_θ Σ_i L(θ; y_i)。对于可微程序（如深度神经网络），我们通过自动微分计算梯度并使用优化算法（如SGD、Adam）求解。

有趣的是，优化领域已经发生了变化。机器学习界发现，像BFGS这样的经典优化器（无需学习率、自动停止）在大规模深度学习中不再流行。取而代之的是SGD及其变体，但我们需要手动调整学习率、监控训练过程。

这是因为在实践中，没人直接优化整个损失函数。由于数据量 N 极大，我们使用随机批量来估计梯度和损失：
∇̂L(θ) ≈ (1/M) Σ_{j in batch} ∇L(θ; y_j)
此时，我们得到的梯度 g 是一个随机变量，其均值是我们关心的真实梯度，但还存在方差 Σ。

然而，当前的主流框架（如PyTorch）只返回梯度的均值（期望），而忽略了其方差（不确定性）信息。计算方差在技术上并不昂贵（只需对已有梯度张量进行平方和求和），但当前的优化器设计哲学仍停留在“处理确定值”的经典数值分析时代，没有将梯度视为一个随机变量来设计更新规则。

Frank, Lucas, Agustinus 和 Julia 将讲授，如果我们正视计算中的这种随机性，并利用梯度的方差等统计量，我们可以更深入地理解深度网络的训练动态，构建更稳定、自动化的训练方案，甚至无需额外成本地为网络输出提供不确定性估计（贝叶斯深度网络）。我们需要将优化重新构想为：利用随机变量探针，寻找一个泛化良好的点，而不仅仅是寻找一个函数的极小值。

核心理念：作为学习机器的数值算法

贯穿以上所有例子的一个共同主题是：数值算法本身就是一种学习机器。

我最初将数值算法定义为“估计一个可能出错的未知量的方法”。它们估计一个隐变量 z（如最小值的位置、积分的值、模拟的解），通过评估一些易处理的观测值 c（如积分点处的函数值、矩阵向量积、批量梯度）来进行估计。

这与贝叶斯推断完全同构：在给定观测数据 c 的条件下，估计隐变量 z 的后验分布。唯一的区别在于，数据的来源不是硬盘，而是处理器（CPU/GPU）。数值算法将计算机作为数据源。

因此，我们可以且应该用机器学习的语言（概率、推断、不确定性）来理解和设计数值算法。这类新算法被称为概率数值算法。它们将计算任务表述为概率分布，将CPU/GPU作为数据源来改进对解的估计，并返回一个包含不确定性的后验分布。

课程总结与特点

本节课中我们一起学习了数值算法对于机器学习的基础性作用，并预览了本课程将如何以全新的视角对待它们。

本课程不是传统的数值分析课：

  • 视角独特：我们将数值算法视为学习机器，使用机器学习的语言进行思考。
  • 注重实践：你将编写大量代码（Python/Julia），通过实践理解算法。
  • 前沿探索：这是全球首个专注于“机器学习数值方法”的课程之一，你将接触到该领域最新的思想。
  • 专家教学：课程由该领域深入研究的博士生们主讲，你可以向他们提出最深入的问题。

通过拥抱“计算即学习”这一观点，我们可以构建出运行更快、更可靠、更易于使用且能提供丰富信息的机器学习解决方案。我们不会抛弃经典的数值算法（那是超过100年的宝贵积累），但我们将摆脱经典数值分析的固有视角，转而采用机器学习的视角，将这些算法本身视为需要理解和改进的学习机器。

希望你能加入我们，共同参与这次激动人心的探索。

002：数值线性代数

在本节课中，我们将学习数值线性代数在机器学习中的核心作用。我们将首先探讨为何需要关注数值线性代数，然后介绍一些对机器学习至关重要的基本任务。最后，我们将通过一个实际例子，学习如何正确实现高斯过程回归，并了解算法实现细节对性能的巨大影响。

为何关注数值线性代数？

数值线性代数是机器学习的基石，是工具箱中最基础的工具集之一。从幻灯片上展示的众多不同领域可以看出，这个概念非常重要。

以下是几个快速示例，用以说明数值线性代数在机器学习中的普遍性：

  • 基础概率论：多元正态分布的概率密度函数中，已经包含了协方差矩阵的逆和行列式。这两个概念是本节课的重点。
  • 一般线性模型：在线性回归中，若使用平方损失并求解相应的优化问题，会出现矩阵的逆。这个矩阵通常具有特定结构（如对称正定），并且是外积形式。
  • 主成分分析：其核心是计算截断特征分解，即将数据协方差矩阵分解为三个具有特定结构的矩阵。这里也涉及矩阵向量积这一基本操作。
  • 高斯过程：这是一种用于回归的贝叶斯非参数方法。其核心公式中同样涉及矩阵的逆。这里的矩阵（核矩阵）具有生成性信息（由一个函数生成每个元素）、对称性和正定性。
  • 微分方程：线性偏微分方程的解可以表示为无限维函数空间中的线性系统。通过离散化，可以将其转化为有限维线性系统，该矩阵通常很大且高度结构化。
  • 优化：大多数局部优化方法（如牛顿法）的更新规则涉及矩阵（如海森矩阵或其逆）与向量（如梯度）的乘积。在深度学习中，由于使用小批量数据，这些矩阵和向量通常是带噪声的估计。
  • 深度学习：前馈神经网络的定义中包含矩阵向量积。反向传播本质上是高效的向量-雅可比积计算。在贝叶斯深度学习中，拉普拉斯近似需要计算并反转神经网络参数的海森矩阵。

机器学习中的核心线性代数任务

从上述例子可以看出，有四个基本线性代数操作在机器学习中无处不在，至关重要：

  1. 矩阵向量乘法：在深度学习和优化中需要大量使用，要求计算时间和内存高效。
  2. 线性系统求解：出现在高斯密度、高斯过程、线性回归、优化和牛顿方法中。
  3. 矩阵分解：用于解决PCA等问题，将矩阵分解为具有特定结构的矩阵（如下三角和上三角矩阵的乘积、奇异值分解、QR分解）。这对于高效进行高斯过程推断和卡尔曼滤波至关重要。
  4. 对数行列式：在高斯密度和高斯过程回归中需要。对于对称正定矩阵，对数行列式可以转化为矩阵对数的迹。这在GP超参数优化中尤其重要。

算法实现的重要性

在计算机上实现数值运算的方式至关重要。即使对于看似简单的任务，错误的实现也可能导致严重问题。

考虑一个任务：给定矩阵 A ∈ ℝ^(n×k)（其中 n >> k）和向量 x ∈ ℝ^k，计算 (A A^T) x。

一个直接的Python实现可能是 A @ A.T @ x。然而，由于运算符的从左到右结合性，这会先计算 A A^T（一个 n×n 矩阵），然后再与 x 相乘。当 n 很大时（例如 10^5），这个中间矩阵将占用约 80 GB 内存，这是不可行的。

正确的实现是利用结合律：A @ (A.T @ x)。这样先计算 A.T @ x（得到一个 k 维向量），再与 A 相乘，内存消耗仅为 O(nk)，计算也快得多。

这个例子说明，实现数学运算的算法对性能、稳定性和内存消耗有巨大影响。

从图形角度看，A A^T 是一个低秩矩阵（秩最多为 k）。对于低秩矩阵，存在比先形成完整矩阵再相乘更高效的方法。

这引出了一个更重要的观点：如果你知道矩阵的结构，就应该利用它。算法本身无法自动识别结构，需要你告知算法或使用专门的方法。

利用矩阵结构

低秩结构只是一个例子。一旦计算出了完整的 A A^T，再想事后识别其低秩结构是昂贵的。因此，如果你知道矩阵是低秩的，就应该使用专门处理低秩矩阵的方法。

矩阵结构多种多样，以下是一个不完整的列表：

  • 一般矩阵：无特殊结构。矩阵向量积的时间和空间复杂度至少为 O(mn)。
  • 低秩+对角矩阵：例如，带有 L2 正则化的回归问题中会出现。复杂度可降至 O(nk)，其中 k 是低秩部分的秩。
  • 稀疏矩阵：复杂度取决于非零元素的数量。表示方式也会影响空间需求。
  • 核矩阵：具有生成性信息（由核函数生成）、对称性和正定性。虽然矩阵向量积的时间没有加速，但可以节省大量内存（存储数据点和核函数，而非整个矩阵，空间复杂度从 O(n²) 降至 O(nd)，其中 d 是输入维度）。
  • 计算图相关的矩阵：如神经网络中的雅可比矩阵。在特定假设下，可以加速计算。

高斯过程回归：一个具体案例

现在，我们聚焦于一个具体算法：高斯过程回归。本节内容也适用于其他方法，如核岭回归。

问题回顾

目标是学习一个未知函数 f*: ℝ^d → ℝ。作为一种贝叶斯回归方法，我们在未知函数上放置一个先验（高斯过程），其均值函数为 μ，核函数为 k。我们观测到带噪声的函数值数据 y。应用贝叶斯规则后，我们得到后验高斯过程，其后验均值和协方差函数如下：

后验均值: μ* = K* (K + σ² I)^(-1) y
后验协方差: k* = K** - K* (K + σ² I)^(-1) K*^T

其中，K 是训练数据点之间的核矩阵，K* 是训练点与测试点之间的核矩阵，**K**** 是测试点之间的核矩阵。

为了进行超参数优化，我们还需要计算对数边际似然（用于模型选择）：
log p(y | θ) = -1/2 [ y^T (K + σ² I)^(-1) y + log det(K + σ² I) + n log(2π) ]

这里，θ 代表核函数的超参数（如输出尺度、长度尺度）。

计算挑战与所需操作

进行高斯过程回归需要：

  1. 求解线性系统：(K + σ² I) w = y （一次求解，用于后验均值）。
  2. 求解线性系统：(K + σ² I) W = K*^T （M 次求解，M 是测试点数量，用于后验协方差）。
  3. 计算 log det(K + σ² I) （用于超参数优化）。

利用结构：核矩阵的性质

我们的核矩阵 K + σ² I 是对称正定的，并且具有生成性信息。对于大规模数据集，我们需要利用这些结构来设计高效算法。

矩阵分解：LU 与 Cholesky

求解线性系统的基础是矩阵分解。我们首先介绍 LU 分解。

前向与后向替换

假设我们要解 L y = b，其中 L 是下三角矩阵。我们可以使用前向替换算法（递归/分治法）：

  1. 将 L 和 b 分解为块形式。
  2. 递归求解较小的下三角系统。
  3. 回代求解最后一个变量。
    其时间复杂度为 O(n²)。类似地，对于上三角矩阵 U x = y，可以使用后向替换。

LU 分解

对于一般可逆方阵 A，我们可以将其分解为 A = L U，其中 L 是单位下三角矩阵（对角线为1），U 是上三角矩阵。分解过程也采用递归：

  1. 选取一个主元（通常是左上角元素）。
  2. 计算该列下方的乘子向量 l。
  3. 更新右下角的子矩阵：A' = A' - l u^T。
  4. 递归地对子矩阵进行 LU 分解。

如果主元为零，需要进行部分选主元（行交换），以确保数值稳定性。通常选择绝对值最大的元素作为主元。加上行交换后，分解可表示为 P A = L U，其中 P 是置换矩阵。

LU 分解的计算成本约为 (2/3)n³ 次浮点运算。

利用 LU 分解求解和计算行列式

一旦得到 A = L U，求解 A x = b 就转化为：

  1. 解 L y = P b （前向替换，O(n²)）。
  2. 解 U x = y （后向替换，O(n²)）。

对于多个具有相同矩阵 A 但不同右侧向量 b 的线性系统，我们只需支付一次 O(n³) 的分解成本，之后每个系统只需 O(n²) 的替换成本。这在高斯过程回归中非常有用。

此外，利用 LU 分解可以高效计算行列式：
det(A) = det(L) det(U) = (∏ᵢ Lᵢᵢ) (∏ᵢ Uᵢᵢ) = ∏ᵢ Uᵢᵢ
因为 L 是单位下三角矩阵，其行列式为 1。计算成本为 O(n)。

Cholesky 分解：针对对称正定矩阵

对于对称正定矩阵（如 K + σ² I），有更高效的 Cholesky 分解：A = L L^T，其中 L 是下三角矩阵。

构造方式与 LU 分解类似，但利用了对称性。由于矩阵正定，所有主对角元都为正，因此不需要选主元。Cholesky 分解的计算成本约为 (1/3)n³，比 LU 分解快一倍，并且更节省内存（只需存储 L）。

应用于高斯过程回归

使用 Cholesky 分解 K + σ² I = L L^T，后验均值和协方差的计算变为：

  • 后验均值：μ* = K* (L^T \ (L \ y))。这里需要一次前向替换和一次后向替换。
  • 后验协方差：类似地，利用分解进行计算。
  • 对数行列式：log det(K + σ² I) = 2 ∑ᵢ log Lᵢᵢ。

扩展到大规模数据集

当数据集很大（例如 n=100,000）时，核矩阵无法放入内存，直接 Cholesky 分解（O(n³)）也不可行。有两种主要思路来解决这个问题。

方法一：利用核结构（精确方法）

选择具有特殊结构的核函数，使得线性代数运算大大简化。

例如，使用参数化先验（有限维特征映射）。假设我们有 p 个特征函数（p << n），权重先验为 w ~ N(0, Σ)。那么对应的核函数为 k(x, x') = φ(x)^T Σ φ(x')。此时的核矩阵为：
K = Φ Σ Φ^T
其中 Φ 是 n×p 的特征矩阵。那么 K + σ² I 具有“低秩 + 对角”的结构。

为了高效求逆，可以使用 矩阵求逆引理（Woodbury 恒等式）：
(Φ Σ Φ^T + σ² I)^(-1) = (1/σ²) I - (1/σ²) Φ [ Σ^(-1) + (1/σ²) Φ^T Φ ]^(-1) Φ^T (1/σ²)

关键点在于，括号内的矩阵是 p×p 的，而 p 很小。因此，求逆成本从 O(n³) 降为 O(np² + p³)。我们甚至不需要显式形成整个核矩阵。类似地，矩阵行列式引理可以高效计算对数行列式。

方法二：近似求逆（迭代/近似方法）

另一种思路是进行近似。观察 Cholesky 分解算法，它迭代地处理数据点，逐行生成下三角因子 L。我们可以提前终止算法，得到一个低秩近似。

然而，按原始顺序处理数据点，近似的收敛速度可能很慢。解决方案是选择处理数据点的顺序。这类似于在 Cholesky 分解中引入完全选主元（同时进行行和列置换），称为带主元的 Cholesky 分解。

通过智能地选择顺序（例如，基于空间覆盖或不确定性），我们可以用更少的迭代步数获得更好的低秩近似。这可以看作一个在线学习或主动学习问题：在每一步，根据当前已获得的信息，选择下一个最具信息量的数据点进行处理。

这种低秩近似不仅可以用于加速求解，还可以用于从后验分布中高效采样（因为 L 本身就是一个矩阵平方根）。

总结

本节课我们一起学习了以下内容：

  1. 数值线性代数是机器学习的基石，出现在几乎所有机器学习算法中。
  2. 作为机器学习从业者，我们需要关心数值实现，因为利用结构可以极大提升算法速度，并且实现细节需要非常小心。
  3. 我们重点讨论了四个核心任务：
    • 高效的矩阵向量乘法。
    • 线性系统求解。
    • 矩阵分解（如 LU、Cholesky）。
    • 对数行列式和矩阵平方根的计算（用于超参数优化和采样）。
  4. 我们通过多个例子看到，利用矩阵结构对于设计快速的机器学习算法至关重要。
  5. 我们以高斯过程回归为例，深入探讨了如何应用这些线性代数工具，并讨论了将其扩展到大规模数据集的两种思路：利用核的精确结构和使用近似方法。

在接下来的课程中，我们将继续探索结构如何影响高斯过程回归及其他机器学习模型。

感谢观看，下节课见。

003：高斯过程的扩展 🚀

在本节课中，我们将学习如何将高斯过程扩展到大型数据集。我们将看到，通过使用迭代方法求解线性系统，可以显著降低计算复杂度，从立方时间降至平方时间，甚至在线性时间内完成推理。

高斯过程回顾 📚

上一节我们介绍了使用楚列斯基分解求解线性系统的教科书方法，并将其应用于高斯过程回归。本节中，我们将探讨解决此类核矩阵线性系统的现代方法，并将其应用于扩展高斯过程。

回归的主要目标是学习一个未知函数 f*。我们期望模型能够：

  • 泛化到未见过的数据。
  • 简单可解释，便于理解模型所学内容。
  • 提供不确定性估计，以了解预测的可信度。
  • 计算快速高效。

高斯过程能同时满足以上所有要求。然而，正如上周所学，其训练和推理的计算成本非常高。高斯过程后验的均值和协方差计算都涉及对大小为 n x n 的格拉姆矩阵求逆，其中 n 是数据点数量。

对于模型选择（如优化核超参数），我们需要计算对数边际似然及其梯度，这同样涉及矩阵求逆和行列式计算。

楚列斯基分解可以精确计算这些量。然而，对于大型数据集（例如 n=10,000），形成核矩阵需要约 80 GB 内存，而分解本身在标准 CPU 上可能需要超过 80 小时。显然，楚列斯基分解在时间和空间上对于我们所关心的大型数据集都是不可行的。

迭代方法作为替代方案 🔄

为了应对上述挑战，我们将讨论迭代方法，目标是将复杂度降低到数据点数量的平方级别。

我们之前看到的楚列斯基分解的迭代形式，在提前停止时，可以得到核矩阵的低秩近似。但对于高斯过程回归，我们需要的是逆矩阵的近似，而非核矩阵本身的近似。

关键问题是：我们能否修改迭代楚列斯基算法，使其在构建矩阵近似 L_i L_i^T 的同时，也构建一个逆矩阵的近似 C_i，使得 C_i b ≈ A^{-1} b？

答案是肯定的。通过追踪额外的量，我们可以递归地构建一个同样低秩的逆矩阵近似 C_i。其每轮迭代的复杂度仍然是 O(n^2)，与仅近似矩阵本身的开销渐近相同。

让我们看看将其应用于高斯过程推理的效果。我们只需用逆近似 C 替换精确的核矩阵逆 K^{-1}。

以下是使用部分楚列斯基（迭代3次）的示例，其中我们按顺序选择了前三个数据点（蓝色点）。可以看到，近似效果并不理想，它本质上类似于只在这三个数据点上进行高斯过程回归。

选择哪些数据点（即“动作”或“枢轴”）至关重要。如果随机选择数据点，即使同样迭代10次，后验均值和不确定性的近似效果也会好得多。

因此，可以将选择数据点的过程视为一种学习行为：一个更智能的策略能让我们更快地收敛到真实的后验。

共轭梯度法：更智能的“动作”选择 ⚡

上一节我们看到了选择不同“动作”对近似速度的影响。本节中，我们来看看一种名为共轭梯度法的方法，它最初就是为求解对称正定线性系统（如高斯过程回归中的系统）而设计的。

我们的目标是求解 A x = b，其中 A 是核矩阵。我们希望用尽可能少的矩阵-向量乘法（每次 O(n^2)）来快速近似后验。

这个问题可以重新表述为一个二次优化问题：最小化函数 f(x) = (1/2) x^T A x - b^T x。其梯度为 ∇f(x) = A x - b，设梯度为零即等价于求解 A x = b。我们定义残差 r = b - A x，它也是当前点的负梯度方向。

如果使用最速下降法（沿负梯度方向），收敛路径可能会呈“之”字形。共轭梯度法的核心思想是：选择一组共轭方向（在由 A 诱导的内积下正交），使得算法能在最多 n 步内收敛（精确算术下）。

共轭梯度法可以解释为一种学习逆矩阵的算法。它同样会构建一个低秩的逆矩阵近似 C_i，其更新形式为外积 d_i d_i^T 的求和，其中 d_i 是共轭方向。

与部分楚列斯基法对比，两者结构相似：

  • 都选择一个“动作”向量（楚列斯基选单位向量，共轭梯度选残差）。
  • 都对该动作进行修正以获得搜索方向。
  • 都通过低秩更新（秩1）来估计逆矩阵。

对于高斯过程推理，共轭梯度法（使用残差作为动作）通常收敛得更快。因为残差向量 r 会全局地探测所有数据点，其大小反映了当前估计与数据的差距，这正是一个高效的探索策略。

收敛速度与预条件处理 🏎️

上一节我们直观地看到共轭梯度法收敛更快。本节中，我们具体分析其收敛速度，并介绍加速技巧——预条件处理。

共轭梯度法的收敛速度取决于矩阵 A 的条件数 κ(A)，即最大特征值与最小特征值之比。条件数越接近1，收敛越快。

对于高斯过程中的矩阵 K + σ^2 I，其特征值为 λ_i + σ^2。其条件数为 (λ_max + σ^2) / (λ_min + σ^2)。

  • 当观测噪声 σ^2 很小时，条件数接近 λ_max / λ_min。对于大多数核函数，这个值可能非常大（病态），导致收敛缓慢。
  • 当 σ^2 很大时，条件数接近1，但这也意味着我们假设观测噪声很大，从数据中学到的东西变少。

这就形成了一个权衡：更小的噪声有助于学习，但会使数值问题更难解。为了解决病态问题，我们引入预条件处理。

其思想是：我们有一个易于求逆的矩阵 P 来近似 A。我们求解等价的线性系统 P^{-1} A x = P^{-1} b。如果 P 选得好，新系统的条件数 κ(P^{-1} A) 会远小于原系统的条件数 κ(A)，从而极大加速共轭梯度法的收敛。

几何上，预条件子 P 的作用是“扭曲”空间，使二次型的等高线更接近圆形，从而优化路径更直接。

预条件子可以视为关于问题的先验知识。对于核矩阵，一个有效的预条件子是使用部分楚列斯基分解得到的低秩加对角矩阵 P = L L^T + σ^2 I。由于其结构特殊，我们可以利用矩阵求逆引理高效计算 P^{-1}。

实践表明，即使是使用随机选取数据点构建的简单部分楚列斯基预条件子，也能显著加速共轭梯度法。例如，对于一个 n=100,000 的数据集，使用预条件的共轭梯度法可以在约7分钟内将残差降到很低，而精确的楚列斯基分解则需要83小时。

超参数优化：对数行列式与迹估计 📈

到目前为止，我们主要讨论了如何求解线性系统以获得后验均值。但对于超参数优化，我们还需要计算对数边际似然及其梯度。

对数边际似然的梯度包含两项：一项涉及线性系统求解（我们已经能用共轭梯度法处理），另一项涉及迹的计算：tr(K^{-1} ∂K/∂θ)。

直接计算迹需要 O(n^3) 复杂度。我们希望用少量 O(n^2) 的矩阵-向量乘法来近似它。迹可以写成 tr(A) = sum_i e_i^T A e_i，其中 e_i 是单位向量。但这仍需要 n 次矩阵-向量乘法。

我们采用随机化迹估计。其核心思想是：如果 z 是一个零均值、协方差为 I 的随机向量，那么 tr(A) = E[ n * (z^T A z) ]。因此，我们可以通过蒙特卡洛方法，用少量随机向量 z_i 的平均来估计迹：tr(A) ≈ (n/L) * sum_{i=1}^L z_i^T A z_i。

然而，简单的蒙特卡洛估计器方差较大，可能需要数万个随机向量才能达到可接受的精度，这同样代价高昂。

解决方案依然是利用预条件子。我们可以将对数行列式重写为：
log det(K) = log det(P) + tr( log(K) - log(P) )
其中 P 是我们的预条件子。第一项 log det(P) 可以精确计算（例如，对于楚列斯基分解的 P，就是其对角线乘积的对数）。我们只需要估计第二项的迹。由于 P 接近 K，log(K) - log(P) 是一个小量，因此其迹的估计方差会大大降低，从而可以用少得多的随机向量（例如几十个）达到高精度。

这种预条件化的随机迹估计策略，与预条件共轭梯度法加速线性求解的思路一脉相承，都是通过注入关于问题的先验知识（预条件子）来提升效率。现代高斯过程库（如GPyTorch）正是使用这种组合方法来实现大规模数据集的超参数优化。

线性时间方法：变分近似与诱导点 🏃

上一节我们实现了平方时间复杂度的推理。本节最后，我们简要探讨能否达到线性时间复杂度 O(n)。

其核心思想是：数据集中往往存在大量冗余或相似的点。我们可以尝试将数据集压缩成一个更小的诱导点集 Z（大小为 m << n），用它来总结原始数据。

这通过变分近似实现：我们直接用一个依赖于诱导点 Z 的、不同的高斯过程来近似真实后验。这个近似后验的形式与标准高斯过程后验类似，但所有涉及 n x n 矩阵的操作都被替换为涉及 m x m 矩阵或 m x n 矩阵的操作。

我们需要优化三个量来使变分分布接近真实后验：

  1. 诱导点的位置 Z。
  2. 一个类似于观测值的向量 μ。
  3. 一个校正不确定性的项 Σ。

通过最小化变分分布与真实后验之间的KL散度，我们可以优化这些参数。其主要计算成本为 O(m^2 n)，当 m 固定时，即为 O(n)。

这种方法（如随机变分高斯过程）速度极快，并且适合使用随机梯度下降进行优化。然而，它有一个潜在问题：不确定性量化可能不准确。特别是在远离诱导点的区域，近似后验的不确定性可能与真实高斯过程后验相差甚远。

与迭代方法不同，变分方法需要预先设定诱导点数量 m，这决定了近似的“容量”。而迭代方法（如共轭梯度）允许我们随时根据计算预算停止，并得到一个在相应计算量下的最优近似，其不确定性估计通常更可靠。

总结 🎯

本节课我们一起学习了如何将高斯过程扩展到大型数据集：

  1. 近似是必需的：精确计算在内存和时间上对于大数据集都不可行。
  2. 迭代方法是核心工具：如部分楚列斯基法和共轭梯度法，它们允许我们以平方时间复杂度 O(n^2) 近似后验并进行超参数优化。
  3. 算法即学习器：这些迭代方法中选择“动作”（如探测哪些数据点）的过程，可以解释为一种主动学习算法。更智能的动作选择能带来更快收敛。
  4. 预条件处理即先验知识：通过向求解器注入关于问题的近似信息（预条件子），可以显著加速迭代方法，无论是线性求解还是迹估计。
  5. 存在更快的线性时间方法：如变分高斯过程，通过诱导点总结数据，实现 O(n) 复杂度推理。但其不确定性估计可能需要仔细评估。

最终，我们留下一个开放问题：能否设计一种方法，在获得线性时间复杂度的同时，还能提供可靠的不确定性量化？ 这将是未来研究的一个方向。

下节课，我们将更深入地探讨这些数值方法本身如何被理解为“学习机器”，从而统一优化与学习的视角。

004：计算感知高斯过程

在本节课中，我们将学习如何为高斯过程回归构建一种计算感知的近似方法。这种方法的核心在于，无论投入多少计算资源，我们都能精确地量化近似带来的不确定性。我们将从一个新的视角——概率数值计算——来重新审视迭代求解算法，并将其视为一个“学习代理”，从而实现对近似误差的跟踪和量化。

快速回顾：高斯过程与迭代近似

上一讲我们深入探讨了用于高斯过程的迭代近似算法。我们了解到，精确计算高斯过程后验涉及对核矩阵 K 进行求逆，其计算复杂度为 O(n³)，在大规模数据场景下代价高昂。

为了应对这个问题，我们引入了迭代算法，例如部分 Cholesky 分解和共轭梯度法。这些算法通过每次迭代执行一次矩阵-向量乘法（复杂度 O(n²)），来构建对精度矩阵 K⁻¹ 的低秩近似。

一个关键发现是，行动（即算法在每次迭代中选择关注的数据点或方向）的选择至关重要。它类似于机器学习中的主动学习：在信息量更大的位置进行“观测”，能更快地降低我们对真实函数认知的不确定性。

然而，这些快速近似方法（如基于诱导点的稀疏近似）存在一个显著问题：它们提供的不确定性估计可能是失准的。在数据稀疏的区域，模型可能过于自信；而在数据密集的区域，近似后验均值可能偏离真实值，但模型给出的不确定性却很低。这违背了使用高斯过程的初衷之一——获得可靠的不确定性量化。

那么，我们能否在享受计算效率的同时，修复这种不确定性失准的问题呢？这就是本节课要回答的核心问题。

核心思想：将权重估计视为学习问题

让我们重新审视高斯过程后验均值的表达式：

μ(x) = k(x, X) · v

其中，v* = K⁻¹y 被称为表征权重。后验均值可以理解为一系列以数据点为中心的核函数 k(x, x_j) 的线性组合，而 v* 就是这些组合的权重。

当我们使用迭代求解器时，我们得到的是表征权重 v* 的一个近似值 v_i。这导致我们的后验均值预测 μ_i(x) = k(x, X) · v_i 也是近似的。

核心思路：如果我们能以概率化的方式，量化近似权重 v_i 与真实权重 v* 之间的误差，并将此误差融入到我们的预测不确定性中，那么即使在近似计算下，我们也能获得校准良好的不确定性估计。

构建概率学习算法

为了实现上述思路，我们为真实表征权重 v* 建立一个高斯先验信念。这个信念包含一个均值估计 v_i（我们当前的猜测）和一个协方差矩阵 Σ_i（对我们猜测的不确定性的度量）。Σ_i 越大，表示我们认为 v_i 可能离 v* 越远。

我们如何更新这个信念呢？我们无法直接观测 v，但可以间接地通过计算残差*来获取信息。在第 i 步，我们选择一个“行动”向量 s_i。我们可以计算一个标量观测值：

r_i = s_iᵀ K (v - v_{i-1})*

这个观测值本质上衡量了在当前行动方向 s_i 上，我们的猜测与真实解之间的距离。这正是一个经典的线性高斯推断问题：我们有一个关于 v* 的高斯先验，并通过一个线性变换 s_iᵀ K 观测到了带有噪声（这里噪声为0）的数据 r_i。

利用高斯推断的公式，我们可以更新对 v* 的信念。更新后的后验均值 v_i 会朝着真实解方向移动，而后验协方差 Σ_i 会减小，反映出我们获得了新的信息。

这个更新过程在数学形式上，与此前见过的部分 Cholesky 和共轭梯度算法惊人地相似。行动 s_i 的选择（例如，选择单位向量或残差向量）直接决定了具体是哪种算法。

解决先验难题与不确定性传播

现在出现了一个循环论证：为了初始化对 v* 的信念，我们需要一个先验协方差 Σ_0。合理的先验应该是 Σ_0 = K⁻¹，但这恰恰是我们试图避免计算的矩阵！

关键的突破点在于，我们最终的目标不是 v* 本身，而是高斯过程的后验预测。当我们对 v* 具有高斯信念时，后验均值 μ(x) = k(x, X) v 就是一个高斯随机变量的线性变换。

因此，我们可以将关于 v* 的不确定性（即计算误差）边缘化到预测中。令人惊喜的是，在最终的预测协方差表达式里，所有需要 K⁻¹ 的项都相互抵消了。我们最终得到的组合协方差为：

k_*(x, x) = k(x, x) - k(x, X) C_i k(X, x) + σ²

其中 C_i 是我们的算法在迭代中自然构建的、对 K⁻¹ 的低秩近似矩阵。我们完全不需要显式计算 K⁻¹。

这个组合协方差包含两部分：

  1. 数学不确定性：源于数据有限，即使精确计算也无法消除的不确定性。
  2. 计算不确定性：源于我们使用近似算法而非精确解所引入的额外不确定性。

随着迭代次数增加，C_i 越来越接近 K⁻¹，计算不确定性逐渐减小，组合协方差收敛到精确的数学后验协方差。

算法：InterGP 与不同行动策略

我们将这个框架实例化为一个名为 InterGP 的算法。其步骤与共轭梯度法等类似，但多了一个显式维护不确定性估计 Σ_i（通过 C_i 体现）的视角。

不同的行动策略 s_i 会生成 InterGP 的不同变体：

  • 单位向量行动：恢复为部分 Cholesky 或子集数据方法。计算集中在选定的数据点上，不确定性降低是局部的。
  • 残差向量行动：恢复为共轭梯度法。计算更全局，通常能以更少的迭代次数接近真实解。
  • 核向量行动：行动 s_i 是一个核函数值向量，类似于稀疏高斯过程中诱导点的思想。这可以修复传统稀疏方法不确定性失准的问题，但每次迭代仍需 O(n²) 计算。

控制计算成本与哲学视角

一个重要的实践发现是：如果我们设计的行动向量 s_i 是稀疏的（即大部分元素为零），那么算法中所有涉及 K 的矩阵-向量乘法的计算成本，就只取决于行动向量非零元素的数量 k，复杂度变为 O(k²)。

这意味着我们可以通过设计行动策略，来自由控制每次迭代的计算预算，同时依然保持精确的不确定性量化特性。从某种意义上说，这模糊了“数据”和“计算”的界限：只有那些被我们的计算行动“处理”过的数据，才真正影响了我们的后验信念。这提供了一个更贴合实际计算场景的建模视角。

理论保证：精确的不确定性量化

InterGP 方法不仅直观，更有坚实的理论保证。它提供的组合标准差，可以作为预测误差的可靠上界。

具体来说，对于任何可能生成观测数据的潜在真实函数 g，其函数值 g(x) 与我们近似后验均值 μ_i(x) 之间的差距，几乎必然被我们的组合标准差所界定。这确保了：

  1. 真实函数几乎必然落在由“近似后验均值 ± 组合不确定性”构成的区间内。
  2. 精确的数学后验均值几乎必然落在由“近似后验均值 ± 计算不确定性”构成的区间内。

因此，我们实现了最初的承诺：在任意给定的计算预算下，进行精确的不确定性量化。我们可以像信任精确高斯过程一样，信任 InterGP 提供的不确定性估计。

总结

本节课中，我们一起学习了一种为高斯过程回归构建计算感知近似的方法——InterGP。核心要点如下：

  • 新视角：将迭代求解器重新诠释为对表征权重的概率化学习算法。
  • 误差量化：该框架能自然跟踪并量化由近似计算引入的误差（计算不确定性）。
  • 不确定性传播：通过边缘化，将计算不确定性融入最终预测，得到组合不确定性。
  • 灵活性：通过选择不同的行动策略，可以恢复经典算法（如 CG），也可以设计新算法，并能控制每次迭代的计算成本。
  • 理论保证：该方法提供的不确定性估计是校准良好的，为预测误差提供了可靠上界。

这种方法弥合了高效计算与可靠统计推断之间的鸿沟，是概率数值计算思想的一个成功范例。下一讲，我们将探讨如何将物理定律等先验知识融入到高斯过程建模中。

005：状态空间模型

在本节课中，我们将要学习一种用于描述和估计动态系统的强大工具——状态空间模型。我们将从基本概念开始，逐步构建一个完整的算法框架，用于从观测数据中在线估计系统的隐藏状态。

概述：什么是动态系统？

上一讲我们学习了高斯过程，它用于描述固定数据集上的函数。本节中，我们将开启新的篇章，探讨动态系统。动态系统是指状态随时间演化的系统。我们的核心任务是：从观测数据中，估计这种动态系统随时间的演化过程。

我们最终的目标是能够模拟这样的系统。为了实现这个目标，我们首先需要建立一套描述动态系统的语言和算法框架。

动态系统的例子与核心特征

以下是几个动态系统的例子，从中我们可以提取出状态空间模型需要具备的关键特征：

  • COVID-19疫情模型：我们每天可以观测到新增感染病例数（红色曲线），但无法直接测量人群的接触率（β参数）。这体现了系统的部分可观测性。
  • 智能手机方向估计：手机内置的加速度计等传感器会持续产生数据流，我们需要快速、在线地更新对手机方向的估计。这体现了在线估计的需求。
  • 流体动力学模拟：这是一个非常复杂的时空动态系统，其动力学可能是高度非线性的。

从这些例子中，我们总结出状态空间模型的关键特征：

  1. 时间演化：系统状态随时间变化。
  2. 部分可观测性：我们只能观测到系统状态的一部分。
  3. 在线估计：需要算法能够基于最新到达的数据流快速更新估计。
  4. 可能包含复杂非线性：系统的动态可能非常复杂。

状态空间模型：建模语言

现在，我们来看看如何用概率模型来描述这样的系统。我们将使用概率图模型来直观表示。

图中白色节点 X_k 代表系统在时刻 k 的状态。状态的定义是：一组物理量，其规格完全决定了系统的演化。例如：

  • 对于钟摆系统，状态可以是 [角度， 角速度]。
  • 对于疫情模型，状态可以是 [易感者数量， 感染者数量， 康复者数量， 接触率]。
  • 对于手机追踪，状态可以是 [X位置， Y位置， X速度， Y速度]。

我们通过观测变量 Y_k 来获取关于状态 X_k 的信息，但通常只能观测到状态的一部分（例如，只能测位置，不能直接测速度）。

从该图模型中，我们可以读出两个核心性质：

  1. 马尔可夫性质：当前状态 X_k 只依赖于前一个状态 X_{k-1}，与更早的历史无关。用公式表示为：
    p(x_k | x_{1:k-1}) = p(x_k | x_{k-1})
  2. 条件独立观测：当前观测 Y_k 只依赖于当前状态 X_k，与过去和未来的状态及观测无关。即：
    p(y_k | x_{1:k}, y_{1:k-1}) = p(y_k | x_k)

基于以上性质，我们可以正式定义概率状态空间模型。它包含三个部分：

  1. 初始分布：p(x_0)
  2. 动态模型（状态转移模型）：p(x_k | x_{k-1})。这定义了状态如何随时间演化。
  3. 测量模型（似然函数）：p(y_k | x_k)。这定义了给定状态下观测数据的分布。

这个三元组 {p(x_0), p(x_k | x_{k-1}), p(y_k | x_k)} 就是我们描述动态系统的建模语言。

贝叶斯滤波：在线状态估计

有了模型，我们想要解决什么问题？在在线估计的设定下，我们主要计算四种分布：

  1. 预测分布：p(x_k | y_{1:k-1})。基于过去所有观测，预测当前状态。
  2. 滤波分布：p(x_k | y_{1:k})。基于包括当前观测在内的所有过去观测，估计当前状态。这是我们在线估计的核心目标。
  3. 数据似然：p(y_k | y_{1:k-1})。计算当前观测数据的概率，在计算过程中自然得到。
  4. 平滑分布：p(x_k | y_{1:T})。基于全部数据（包括未来数据）重新估计过去某一时刻的状态。这通常在数据收集完成后进行，能得到更优的估计。

滤波和预测分布是递归计算的。我们从初始状态 p(x_0) 开始。

预测步骤：从上一时刻的滤波分布 p(x_{k-1} | y_{1:k-1}) 出发，预测当前状态。
p(x_k | y_{1:k-1}) = ∫ p(x_k | x_{k-1}) p(x_{k-1} | y_{1:k-1}) dx_{k-1}
这个积分方程称为 Chapman-Kolmogorov 方程。

校正（更新）步骤：当新的观测 y_k 到达时，我们用它来更新预测。
p(x_k | y_{1:k}) ∝ p(y_k | x_k) p(x_k | y_{1:k-1})
这本质上是贝叶斯定理：后验 ∝ 似然 × 先验。其中先验就是预测分布，似然就是测量模型。归一化常数就是数据似然 p(y_k | y_{1:k-1})。

将这两个步骤结合起来，就得到了贝叶斯滤波算法的伪代码：

初始化： p(x_0)
循环（对于每个新观测 y_k）：
    预测：计算 p(x_k | y_{1:k-1}) （使用 Chapman-Kolmogorov 方程）
    校正：计算 p(x_k | y_{1:k}) = η * p(y_k | x_k) * p(x_k | y_{1:k-1}) （使用贝叶斯定理）
    返回滤波分布序列 {p(x_k | y_{1:k})}

这是一个前向递归算法，只使用当前和过去的信息，非常适合在线估计。

卡尔曼滤波：线性高斯情形

上述贝叶斯滤波公式中包含积分，对于一般模型很难计算。为了使计算可行，我们引入一个具体的模型假设：线性高斯状态空间模型。

我们假设：

  • 初始分布：x_0 ~ N(μ_0, Σ_0)
  • 动态模型：x_k = A * x_{k-1} + b + q_k，其中 q_k ~ N(0, Q)。A 称为状态转移矩阵，Q 称为过程噪声协方差矩阵。
  • 测量模型：y_k = H * x_k + c + r_k，其中 r_k ~ N(0, R)。H 称为观测矩阵，R 称为观测噪声协方差矩阵。

在线性高斯假设下，所有的预测、滤波、平滑分布都仍然是高斯分布。我们只需要递归地计算它们的均值向量和协方差矩阵即可。这就是著名的卡尔曼滤波算法。

以下是卡尔曼滤波的方程：

预测步骤：
μ_k^- = A * μ_{k-1} + b
Σ_k^- = A * Σ_{k-1} * A^T + Q

校正步骤：

  1. 计算观测预测：ŷ_k = H * μ_k^- + c
  2. 计算新息协方差：S_k = H * Σ_k^- * H^T + R
  3. 计算卡尔曼增益：K_k = Σ_k^- * H^T * S_k^{-1}
  4. 更新状态估计（均值）：μ_k = μ_k^- + K_k * (y_k - ŷ_k)
  5. 更新不确定性（协方差）：Σ_k = Σ_k^- - K_k * S_k * K_k^T

卡尔曼增益 K_k 是一个关键矩阵，它决定了我们应该在多大程度上信任新的观测数据。残差 (y_k - ŷ_k) 代表了新观测与我们的预测之间的差异。

数据似然也自然地得到了：p(y_k | y_{1:k-1}) = N(ŷ_k, S_k)。

卡尔曼滤波算法伪代码如下：

初始化： μ_0, Σ_0
循环（对于每个新观测 y_k）：
    # 预测
    μ_k^- = A * μ_{k-1} + b
    Σ_k^- = A * Σ_{k-1} * A^T + Q
    # 校正
    ŷ_k = H * μ_k^- + c
    S_k = H * Σ_k^- * H^T + R
    K_k = Σ_k^- * H^T * inv(S_k) # 实践中应避免直接求逆，而应解线性系统
    μ_k = μ_k^- + K_k * (y_k - ŷ_k)
    Σ_k = Σ_k^- - K_k * S_k * K_k^T
    保存 (μ_k, Σ_k)
返回滤波轨迹 {(μ_k, Σ_k)}

实例演示：卡尔曼滤波用于车辆追踪

让我们通过一个具体例子来理解卡尔曼滤波。假设我们要追踪一辆在二维平面上行驶的汽车。

  • 状态：x = [位置_x, 位置_y, 速度_x, 速度_y]^T
  • 观测：y = [位置_x, 位置_y]^T （例如来自GPS，带有噪声）
  • 动态模型：假设汽车近似匀速运动（速度可能缓慢变化），用线性模型加噪声来描述。
  • 测量模型：观测矩阵 H 是一个从4维状态中选取前2维（位置）的矩阵。

我们按照卡尔曼滤波的步骤编写代码，使用模拟的汽车轨迹和带噪声的观测。运行滤波算法后，我们可以看到：

  • 在两次观测之间，滤波器的预测（基于运动模型）会形成一条直线轨迹，并且不确定性（协方差）会随时间增长。
  • 当新的GPS观测到达时，滤波器会迅速“修正”其估计，将状态拉向观测值，同时不确定性显著降低。
  • 最终，我们得到一条平滑的、考虑了所有历史观测的估计轨迹，它比单纯的观测数据更接近真实轨迹。

从滤波到平滑：利用全部数据

卡尔曼滤波给出了在线的最优估计，但它只使用了当前和过去的数据。如果我们已经收集完一段时间内的所有数据（例如汽车已到达目的地），我们就可以利用未来的数据来改进对过去某一时刻状态的估计。这个过程称为平滑。

平滑分布为：p(x_k | y_{1:T})，其中 T > k。

对于线性高斯模型，平滑也可以通过高效的递归算法（Rauch–Tung–Striebel smoother）完成，而且平滑分布也是高斯的。平滑是一个后向递归过程：

  1. 我们从最后一个时间点 T 开始，此时平滑分布等于滤波分布：p(x_T | y_{1:T})。
  2. 然后我们向后迭代（k = T-1, T-2, ..., 1），利用以下关系更新：
    p(x_k | y_{1:T}) 可以通过 p(x_k | y_{1:k})（滤波分布）、p(x_{k+1} | x_k)（动态模型）和 p(x_{k+1} | y_{1:T})（下一时刻的平滑分布）计算出来。

具体更新方程涉及一个平滑增益矩阵，其形式与卡尔曼增益类似。平滑算法的伪代码如下：

初始化： μ_T^s = μ_T, Σ_T^s = Σ_T # 从最终滤波结果开始
循环（k 从 T-1 到 0）：
    # 计算平滑增益
    G_k = Σ_k * A^T * inv(Σ_{k+1}^-)
    # 更新平滑均值和协方差
    μ_k^s = μ_k + G_k * (μ_{k+1}^s - μ_{k+1}^-)
    Σ_k^s = Σ_k + G_k * (Σ_{k+1}^s - Σ_{k+1}^-) * G_k^T
返回平滑轨迹 {(μ_k^s, Σ_k^s)}

对车辆追踪例子应用平滑算法后，我们可以看到估计轨迹变得更加平滑和准确，尤其是在观测稀疏的区域，平滑算法利用了整个时间序列的信息进行了更好的插值。

有趣的是，对于线性高斯状态空间模型，平滑后得到的序列 {x_0, x_1, ..., x_T} 的联合后验分布，恰好是一个高斯过程在离散时间点上的边际分布。这意味着我们通过一次前向滤波和一次后向平滑，以线性时间的复杂度，计算出了这个高斯过程后验。

扩展卡尔曼滤波：处理非线性

现实中的动态系统往往是非线性的。例如，钟摆的运动方程就是非线性的。当动态模型 f 或测量模型 h 是非线性函数时，f(x) 或 h(x) 的分布就不再是高斯分布，卡尔曼滤波的闭合形式解也就不再成立。

扩展卡尔曼滤波 是解决此问题的一种常用方法。其核心思想是：对非线性模型进行局部线性化，然后在线性化点附近应用卡尔曼滤波公式。

线性化通过一阶泰勒展开实现。对于非线性函数 g(x)，在点 χ 处的线性化近似为：
g(x) ≈ g(χ) + J_g(χ) * (x - χ)
其中 J_g(χ) 是 g 在 χ 处的雅可比矩阵（导数矩阵）。

EKF 的巧妙之处在于线性化点的选择：

  • 在预测步骤，我们在上一时刻的滤波均值 μ_{k-1} 处线性化动态模型 f。
  • 在校正步骤，我们在当前时刻的预测均值 μ_k^- 处线性化测量模型 h。

这样选择后，EKF 的方程与 KF 非常相似，只需做如下替换：

  • 将 KF 中的线性变换 A * μ_{k-1} + b 替换为非线性函数求值 f(μ_{k-1})。
  • 将 KF 中的矩阵 A 和 H 替换为相应非线性函数在对应线性化点处的雅可比矩阵 J_f 和 J_h。

EKF 的预测和更新方程如下：

预测：
μ_k^- = f(μ_{k-1})
Σ_k^- = J_f(μ_{k-1}) * Σ_{k-1} * J_f(μ_{k-1})^T + Q

校正：
S_k = J_h(μ_k^-) * Σ_k^- * J_h(μ_k^-)^T + R
K_k = Σ_k^- * J_h(μ_k^-)^T * S_k^{-1}
μ_k = μ_k^- + K_k * (y_k - h(μ_k^-))
Σ_k = Σ_k^- - K_k * S_k * K_k^T

同样，也存在扩展卡尔曼平滑器，只需在平滑增益计算中用雅可比矩阵替换原来的转移矩阵 A 即可。

实例演示：EKF用于钟摆跟踪

我们用一个阻尼钟摆的例子来演示EKF。

  • 状态：x = [角度 θ, 角速度 θ_dot]^T
  • 观测：y = θ + 噪声 （只能观测到带有噪声的角度）
  • 动态模型：f(x) 由物理方程给出，包含 sin(θ) 等非线性项。
  • 测量模型：h(x) = θ，这里是线性的。

我们编写代码实现EKF。在预测步骤，我们使用自动微分来计算动态模型 f 在 μ_{k-1} 处的雅可比矩阵 J_f。运行EKF后，我们可以看到：

  • 预测轨迹不再是一条直线，而是体现了钟摆的振荡行为。
  • 在缺乏观测的时间段，不确定性会迅速扩大。
  • 当新的角度观测到达时，估计值会被修正。
  • 应用平滑器后，我们得到了一条非常漂亮且准确的钟摆运动轨迹估计。

其他滤波方法简介

EKF通过线性化来近似非线性模型。还有其他方法：

  • 无迹卡尔曼滤波：UKF采用不同的策略。它不线性化模型，而是通过精心选择一组称为“Sigma点”的确定性样本点来近似状态的分布。将这些点通过真实的非线性函数传播，然后用传播后的点来拟合一个高斯分布。UKF通常比EKF更稳定，尤其对于高度非线性系统。
  • 粒子滤波：当系统不仅是非线性的，而且噪声也是非高斯的时候，我们需要更强大的工具。粒子滤波是一种序列蒙特卡洛方法。它用一组随机样本（粒子）来近似后验分布，每个粒子都有权重。通过重采样等步骤，粒子滤波可以处理非常一般的非线性非高斯模型，但计算成本也更高。

总结

本节课中我们一起学习了状态空间模型这一用于描述动态系统的强大语言。我们的学习路径是：

  1. 我们首先建立了概率状态空间模型的基本框架，明确了马尔可夫性和条件独立观测两个核心假设。
  2. 我们推导了贝叶斯滤波的一般原理，包括预测（Chapman-Kolmogorov方程）和校正（贝叶斯定理）两个递归步骤。
  3. 在线性高斯的特殊情形下，贝叶斯滤波简化为高效且解析的卡尔曼滤波算法，我们详细推导了其方程。
  4. 我们通过车辆追踪的实例演示了卡尔曼滤波的应用。
  5. 为了利用全部数据改进估计，我们介绍了平滑的概念和Rauch–Tung–Striebel平滑算法。
  6. 为了处理现实中的非线性系统，我们引入了扩展卡尔曼滤波，其核心是通过一阶泰勒展开进行局部线性化。
  7. 我们通过钟摆跟踪的实例演示了EKF的应用。
  8. 最后，我们简要介绍了无迹卡尔曼滤波和粒子滤波等其他高级方法。

下节课，Nathaniel将深入讲解如何用数学语言（如微分方程）来描述动态系统本身的物理规律。而在再下一节课，我们将把今天学习的滤波/平滑技术与下一节课的动力学描述结合起来，探索如何使用概率数值方法（贝叶斯滤波与平滑）来模拟复杂的动态系统。敬请期待！

006：常微分方程解法

概述

在本节课中，我们将学习常微分方程（ODEs）及其在机器学习中的应用。我们将探讨如何数值求解ODEs，并了解如何从数据中推断ODE模型的参数。课程内容从ODE的基本定义开始，逐步深入到数值求解算法，最后讨论参数推断问题。

ODE简介与机器学习应用

上一节我们概述了课程内容，本节中我们来看看常微分方程（ODE）的正式定义及其在机器学习中的几个关键应用场景。

一个常微分方程（ODE）通常写作：

ẋ(t) = F(x(t), t)

其中 x 是未知函数（可理解为状态），F 被称为向量场，它描述了状态 x 在时间 t 处的导数。我们通常关注的是在给定初始值 x(0) = x₀ 下的初值问题。

ODE在机器学习中有着广泛的应用，以下是几个核心例子：

  • 扩散模型：这类生成模型的核心依赖于（随机）微分方程来描述从噪声到数据的映射过程。
  • 标准化流：用于建模复杂概率分布，常通过ODE实现从简单分布到复杂分布的连续变换。
  • 神经常微分方程：将残差神经网络重新解释为一个连续动力系统的离散化，引发了大量相关研究。
  • 梯度流：为理解梯度下降等优化算法的行为提供了一个连续的视角。
  • 参数推断：当我们拥有描述数据生成过程的先验知识（例如物理、生物模型）时，可以通过拟合ODE参数来从数据中学习更有意义的结构。

如何数值求解ODE？

上一节我们介绍了ODE及其应用，本节中我们来看看如何数值求解它们。对于大多数ODE，我们无法得到解析解，因此需要数值方法。

求解初值问题的核心思想是逐步外推。假设我们在时间 t 有一个近似解 x̂(t)，目标是计算 t + h 时刻的解。根据积分形式：

x(t+h) = x(t) + ∫[t, t+h] F(x(τ), τ) dτ

数值方法的关键在于如何近似计算这个积分。

欧拉方法

首先，我们介绍两种最基本的数值方法：前向欧拉法和后向欧拉法。它们都源于对积分中被积函数的零阶泰勒近似。

前向欧拉法（显式）：
在积分区间内，用起点 t 处的函数值近似整个区间。这导出了显式更新公式：

x̂(t+h) = x̂(t) + h * F(x̂(t), t)

后向欧拉法（隐式）：
在积分区间内，用终点 t+h 处的函数值近似。这导出了隐式方程：

x̂(t+h) = x̂(t) + h * F(x̂(t+h), t+h)

该方程需要额外的数值方法（如求根算法）来求解 x̂(t+h)。

隐式方法（如后向欧拉）通常比显式方法（如前向欧拉）具有更好的稳定性。稳定性意味着当步长 h 较大或系统本身变化剧烈（刚性系统）时，算法解不会无界地偏离真实解。一个简单的测试方程 ẋ = λx (λ < 0) 可以清晰地展示这一特性：前向欧拉法在步长过大时会发散，而后向欧拉法则能保持稳定。

龙格-库塔方法

上一节我们介绍了基础的欧拉方法，本节中我们来看看更强大和通用的一类算法：龙格-库塔方法。

龙格-库塔方法的动机来源于数值积分（求积法）。它通过计算向量场 F 在多个中间点（求积节点）的值，并以特定权重（求积权重）进行组合，来更精确地近似积分。

一个通用的 s 级显式龙格-库塔方法步骤如下：

  1. 计算中间斜率：
    k₁ = F(x̂(t), t)
k₂ = F(x̂(t) + h * a₂₁ * k₁, t + c₂ * h)
...
k_s = F(x̂(t) + h * (a_{s,1}*k₁ + ... + a_{s,s-1}*k_{s-1}), t + c_s * h)
  2. 组合斜率进行更新：
    x̂(t+h) = x̂(t) + h * (b₁*k₁ + b₂*k₂ + ... + b_s*k_s)

系数 a_{ij}, b_i, c_i 通常用一个称为布彻表的三角形表格来简洁表示。

以下是几个重要的龙格-库塔方法：

  • 前向/后向欧拉法：它们都是一级龙格-库塔法。
  • 显式中点法：一个二级方法，先使用半步长的欧拉步来估计中点斜率，再用该斜率进行更新。
  • 经典四阶龙格-库塔法：一个广泛应用的四阶方法，在精度和计算成本之间取得了良好平衡。
  • Dormand-Prince 5(4) 方法：这是一个自适应步长方法的例子。它同时计算一个五阶和一个四阶的解，两者的差异用于估计局部误差，进而自动调整步长 h 以满足用户指定的精度要求。SciPy 中的 solve_ivp 等求解器默认使用此类方法。

方法阶数：一个方法的阶数 p 描述了其局部截断误差如何随步长缩小而减小，为 O(h^{p+1})。高阶方法通常单步更精确，但计算量也更大，且不一定在所有问题上都更稳定或高效。

参数推断：从数据中学习ODE

上一节我们讨论了如何模拟已知的ODE系统，本节中我们来看看一个更常见的问题：如何从观测数据中学习未知的ODE参数（甚至方程本身）。

我们假设数据由以下隐式生成模型产生：

  1. 状态 x 由含参数 θ 的 ODE 演化：ẋ(t) = F(x(t), t; θ), x(0) = x₀(θ)。
  2. 我们在时间点 t_i 观测到带有噪声的数据：y_i = H * x(t_i) + ε_i, ε_i ~ N(0, Σ)。

目标是基于观测数据 {y_i, t_i} 推断参数 θ。一种直接的方法是寻找最大似然估计，即最小化负对数似然。由于真实解 x(t_i; θ) 无法精确计算，我们用数值解 x̂(t_i; θ) 替代：

L(θ) = -log p({y_i} | θ) ≈ ∑_i || y_i - H * x̂(t_i; θ) ||_Σ²

这本质上是一个最小二乘问题，可以通过优化算法（如L-BFGS）求解。优化过程的内循环需要调用ODE数值求解器来计算当前参数 θ 下的轨迹 x̂ 和损失。

示例：SIRD流行病模型

我们用一个简化的流行病模型（SIRD：易感者-感染者-康复者-死亡者）演示参数推断。模型包含参数：接触率 β、康复率 γ、死亡率 η 以及初始感染人数 I₀。给定模拟生成的带噪声数据，我们可以成功优化这些参数，使模型轨迹很好地拟合数据，并且得到的参数具有明确的物理解释。

扩展：时间相关参数与神经ODE

为了建模更复杂的动态（如COVID-19疫情中随时间变化的接触率），我们可以将常数参数 β 替换为一个函数，例如一个神经网络 β(t; θ)。这时的推断问题就变成了同时学习神经网络的参数 θ。虽然这种方法能提供参数的点估计，但缺乏不确定性量化，且结果严重依赖于神经网络的结构和优化过程。

这引出了下节课的主题：概率数值ODE求解器。我们将结合高斯过程等概率模型，在求解和推断ODE时，不仅能得到预测，还能提供可靠的不确定性估计。

总结

本节课中我们一起学习了常微分方程在机器学习中的重要性。我们探讨了数值求解ODE的基本原理，从简单的欧拉法到更复杂的龙格-库塔法，并了解了显式与隐式方法在稳定性上的差异。最后，我们将ODE模拟与参数推断相结合，展示了如何从数据中学习具有物理意义的模型参数，并简要提到了更复杂的模型（如神经ODE）及其挑战。下节课，我们将从概率视角重新审视ODE求解问题。

007：概率数值常微分方程求解器

概述

在本节课中，我们将学习如何将常微分方程求解问题重新表述为一个贝叶斯状态估计问题。我们将结合前两周关于状态估计和经典ODE求解器的知识，构建一个概率数值ODE求解器。这种求解器不仅能提供近似解，还能给出解的不确定性量化。

数值ODE求解回顾

上一节我们介绍了状态估计的框架，本节我们首先回顾一下经典的数值ODE求解方法。

我们考虑一个常微分方程初值问题：
$$
\dot{x}(t) = f(x(t), t), \quad x(t_0) = x_0
$$
其中，( x(t) ) 是我们希望找到的解，( f ) 是向量场。

数值方法（如欧拉法、龙格-库塔法）的目标是：给定初始估计 ( \hat{x}0 )，通过评估向量场 ( f ) 来生成未来时间点的估计 ( \hat{x} )。

然而，这些方法本质上是估计真实解 ( x(t) )。它们总会引入误差，误差大小通常取决于步长 ( h )。即使步长很小，误差也总是存在。

既然数值ODE求解器本质上是在进行估计，我们很自然地会想：能否将其明确地构建为一个估计问题，并利用我们已有的强大估计工具（如卡尔曼滤波器）来解决它？

贝叶斯状态估计回顾

在构建ODE求解的状态估计模型之前，我们先简要回顾一下相关的工具。

我们之前研究的是非线性高斯状态估计问题。我们有一个状态 ( x ) 和观测 ( y )。生成模型包括：

  1. 描述状态如何随时间演化的动态模型（先验）。
  2. 描述观测如何依赖于状态的观测模型（似然）。

给定一些数据后，我们的任务是进行推断，即计算在给定观测数据条件下状态的后验分布。主要任务包括：

  • 滤波：估计当前状态，给定截至当前的所有信息。
  • 平滑：估计所有时间点的状态，给定全部数据。

我们可以使用扩展卡尔曼滤波器和扩展卡尔曼平滑器等算法来计算这些后验分布。这些算法包含预测步和更新步，概念简单且易于实现。

结合两种思路：概率数值ODE求解器

现在，我们将结合上述两种思路。我们将常微分方程视为一个状态估计问题，并使用贝叶斯滤波和平滑工具来解决它。

我们的目标不再是寻找一个精确满足ODE的 ( x(t) )，而是计算一个后验分布 ( p(x(t) | \text{数据}) )。这里的“数据”包括：

  1. 初始条件 ( x_0 )。
  2. ODE在离散时间点上的信息（即 ( \dot{x}(t_i) = f(x(t_i), t_i) )）。

因此，问题被重新表述为：在给定初始条件和ODE在离散点上定义的“观测”下，推断状态 ( x(t) ) 的后验分布。

为了使用高效的高斯滤波与平滑算法（其计算复杂度为 ( O(n) )），我们需要构建一个状态空间模型。这包括定义先验模型、似然模型和数据。

构建先验模型

首先，我们需要定义状态 ( X(t) ) 是什么。在ODE中，仅仅知道 ( x(t) ) 的当前值不足以完全描述系统；理想情况下，我们需要知道所有阶导数。受泰勒级数展开的启发，我们选择跟踪状态的前 ( q ) 阶导数作为我们的状态向量：
$$
X(t) = [x^{(0)}(t), x^{(1)}(t), ..., x{(q)}(t)]T
$$
其中 ( x^{(k)} ) 表示 ( x ) 的第 ( k ) 阶导数。

当我们用有限项 ( q ) 截断泰勒级数时，会引入一个误差项。为了得到一个概率模型，我们将这个误差项建模为一个随机变量。这就引出了一类称为 ( q ) 次积分维纳过程 的先验模型。

在离散时间设置下，这个先验具有高斯转移密度的形式：
$$
X_{t+h} | X_t \sim \mathcal{N}(A(h) X_t, Q(h))
$$
其中矩阵 ( A(h) ) 和 ( Q(h) ) 有闭合形式。例如，对于 ( q=2 )（跟踪函数值、一阶导和二阶导），( A(h) ) 是一个上三角矩阵，包含泰勒级数系数 ( 1, h, h^2/2 )。

这个先验是高斯-马尔可夫过程，易于处理。我们可以从中采样来可视化其性质：第 ( q ) 阶导数是一个维纳过程（布朗运动），对其积分得到低阶导数。

构建似然与数据模型

先验定义了状态的演化，但尚未包含ODE本身的信息。现在，我们构建似然模型来编码ODE。

原始的ODE是 ( \dot{x}(t) = f(x(t), t) )。利用我们定义的状态 ( X(t) )，我们可以将其等价地写为：
$$
x^{(1)}(t) - f(x^{(0)}(t), t) = 0
$$
我们定义一个测量函数（或信息算子）：
$$
\mathcal{H}(X(t), t) = x^{(1)}(t) - f(x^{(0)}(t), t)
$$
在理想情况下，对于真实解，应有 ( \mathcal{H}(X(t), t) = 0 ) 对所有 ( t ) 成立。

由于我们只在离散时间点 ( t_i ) 上进行计算，我们放松要求，只希望在这些点上满足：
$$
\mathcal{H}(X(t_i), t_i) = 0, \quad i = 1, ..., N
$$
这引导我们定义以下似然模型：我们引入“观测”变量 ( Z_i )，并定义
$$
Z_i | X(t_i) \sim \mathcal{N}(\mathcal{H}(X(t_i), t_i), R)
$$
其中 ( R ) 是观测噪声协方差。然而，由于我们想要精确满足ODE（而不是有噪声的观测），我们将其设为一个无噪声似然，即令 ( R = 0 )。数据点则是 ( z_i = 0 )。因此，我们是在以 ( \mathcal{H}(X(t_i), t_i) = 0 ) 为条件进行“观测”。

最后，我们还需要纳入初始条件 ( x(t_0) = x_0 )。这可以通过另一个无噪声的“观测”来实现：
$$
\mathcal{H}_0(X(t_0)) = x^{(0)}(t_0), \quad \text{数据为 } x_0
$$

完整的状态空间模型与算法

现在，我们拥有了完整的状态空间模型：

  • 先验（动态模型）： ( q ) 次积分维纳过程，定义了 ( X(t) ) 的演化。
  • 似然1（初始条件）： ( \mathcal{H}_0(X(t_0)) = x^{(0)}(t_0) )，数据为 ( x_0 )。
  • 似然2（ODE）： ( \mathcal{H}(X(t_i), t_i) = x^{(1)}(t_i) - f(x^{(0)}(t_i), t_i) )，数据为 ( 0 )。

有了这个模型，我们可以直接应用扩展卡尔曼滤波器（EKF） 和平滑器来进行推断。算法步骤大致如下：

  1. 初始化：使用初始条件似然更新先验分布。
  2. 循环每个时间步：
    • 预测步：使用先验动态模型将状态估计向前推进一步。
    • 更新步：使用ODE定义的似然模型（测量函数 ( \mathcal{H} )）更新状态估计，以“观测”到的 ( 0 ) 为条件。
  3. （可选）进行平滑以获得更优的轨迹估计。

通过实现上述算法，我们得到了一个概率ODE求解器。它不仅输出解的估计（均值轨迹），还输出表示不确定性的协方差。代码实现相对简洁，核心是EKF的预测和更新步骤。

不确定性校准

初始实现的概率求解器给出的不确定性估计通常过于保守（即置信区间过大）。这是因为先验模型中包含一个缩放超参数 ( \sigma^2 )（类似于高斯过程中的输出尺度），它直接影响过程噪声的大小，而我们尚未对其进行设置。

为了获得有意义的、校准良好的不确定性估计，我们需要优化这个超参数 ( \sigma )。这类似于高斯过程中的超参数学习。

在状态估计的框架下，我们可以通过最大化边缘似然来估计 ( \sigma )。边缘似然是在给定超参数下观测到所有数据（即所有 ( z_i = 0 )）的概率。对于扩展卡尔曼滤波器，我们可以近似计算这个量。

有趣的是，对于我们的特定模型，最大似然估计 ( \hat{\sigma} ) 可以有闭合形式的解。更重要的是，可以证明，在EKF推断后对协方差进行简单的重新缩放 ( \Sigma_{\text{new}} = \hat{\sigma}^2 \Sigma_{\text{old}} )，等价于使用最优 ( \hat{\sigma} ) 重新运行整个推断。

因此，校准步骤可以高效地完成：先运行一遍EKF，计算残差，得到 ( \hat{\sigma} )，然后直接重新缩放输出的协方差矩阵。经过校准后，不确定性带能够更紧密地反映真实误差，变得更有信息量。

方法的灵活性与扩展性

概率ODE求解框架的强大之处不仅在于提供不确定性量化，还在于其灵活性。我们可以通过修改状态空间模型，用几乎相同的算法解决各种相关问题。

以下是几个例子：

  1. 二阶ODE： 许多物理系统由二阶方程 ( \ddot{x} = f(x, \dot{x}, t) ) 描述。我们只需调整测量函数为 ( \mathcal{H}(X(t_i), t_i) = x^{(2)}(t_i) - f(x^{(0)}(t_i), x^{(1)}(t_i), t_i) )，并同时加入初始位置和初始速度的条件。算法核心无需改变。

  1. 微分代数方程（DAE）： 形式为 ( M \dot{x} = f(x, t) )，其中 ( M ) 是奇异矩阵。我们可以定义相应的测量函数来编码这个方程，然后同样使用EKF求解。

  2. 带守恒量的ODE： 如果系统存在守恒量（如能量）( g(x, t) = \text{常数} )，我们可以将其作为额外的“观测”添加到模型中。在算法上，这相当于在每个时间点进行两次更新步：一次针对ODE，一次针对守恒量。这有助于数值解更好地保持系统的物理特性。

  3. 隐力推断与数据同化： 我们可以将ODE求解与外部数据拟合结合起来。例如，在流行病学SIRD模型中，将传染率 ( \beta(t) ) 建模为一个高斯过程（潜变量）。状态 ( \tilde{X}(t) ) 现在同时包含ODE解 ( [S, I, R, D] ) 和潜变量 ( \beta )。测量函数同时编码ODE、初始条件以及实际观测数据（如感染人数）。通过一次EKF平滑，我们可以同时推断出模型的轨迹和时变的传染率。

这种灵活性表明，贝叶斯滤波与平滑提供了一个统一的框架，用于表述和解决结合了物理知识（ODE）与观测数据的动态系统推断问题。

总结

本节课中，我们一起学习了如何从概率视角看待常微分方程求解。

  • 我们认识到，经典数值ODE求解器本质上是进行点估计，而贝叶斯方法允许我们进行分布估计，从而提供不确定性量化。
  • 我们构建了一个状态空间模型，其中先验是 ( q ) 次积分维纳过程，似然编码了ODE和初始条件。
  • 我们使用扩展卡尔曼滤波器和平滑器作为推理算法，得到了一个概率数值ODE求解器。
  • 我们讨论了通过最大化边缘似然来校准不确定性估计的重要性，以获得有信息量的置信区间。
  • 最后，我们探讨了该框架的强大灵活性，它可以轻松扩展到求解二阶ODE、DAE、带守恒量的系统，以及进行复杂的潜力推断和数据同化问题。

总之，贝叶斯状态估计为我们提供了一个强大、灵活且统一的工具，用于处理涉及动态系统和不确定性的各种数值计算和推理任务。

008：偏微分方程

在本节课中，我们将学习偏微分方程，并探讨如何将基于PDE的模型整合到机器学习，特别是概率机器学习模型中。我们将通过一个实际的建模示例来使内容更加直观。

什么是偏微分方程及其重要性

首先，我们来明确什么是偏微分方程，以及为什么这些方程很重要。

PDE被用作描述机制性知识的语言。机制性模型描述了生成数据的底层机制。例如，在现实世界或金融市场中，我们可能不知道两个质子如何具体相互作用，但我们知道它们会相互排斥，并且知道这种力的强度。这就是机制性知识。

另一个由PDE描述的物理系统是流体力学领域，即对各种流体的描述。实际上，气候模型和天气模型也基于这些PDE。这里一个非常重要的PDE系统是纳维-斯托克斯方程，用于模拟天气、气候、海洋等系统，例如预测海啸何时会袭击海岸。

从这些情况可以看出，这些模型通常规模非常大，例如模拟整个地球气候或海洋，并且在实际中求解也非常困难。

PDE是一个有趣的话题，因为即使在经过一个世纪甚至更长时间的研究后，PDE的理论和实践仍然是一个高度活跃的研究领域。实际上，著名的千禧年难题之一就是关于纳维-斯托克斯方程，即仅仅证明这组方程是否有解以及解的平滑性如何。这或许已经让你对处理这些模型的难度有了一些了解。

由于非线性或一般的PDE通常很难理解，在本讲座中，我们将稍微限制自己，只讨论一个更简单的类别，即线性PDE。即使我们限制在这个类别，它仍然是一个非常强大的建模语言。例如，我们周围的许多物理过程仍然可以通过线性PDE相当准确地描述。

  • 热传导：热量在金属中的扩散由所谓的热方程描述，这是一个线性PDE。
  • 电磁现象：质子相互作用等由一组称为麦克斯韦方程组的线性PDE描述。
  • 波动力学：水波和空气中传播的波可以非常近似地由波动方程描述，这也是一个线性方程。
  • 布朗运动：布朗运动中粒子的速度分布由所谓的福克-普朗克方程或柯尔莫哥洛夫方程描述。
  • 金融市场：数学金融中著名的布莱克-舒尔斯方程也是这样一个方程。

最后，在实际中处理非线性偏微分方程时，我们可以使用这些非线性方程的线性近似，例如通过迭代线性化进行数值模拟。

将PDE模型与机器学习结合

有了这些动机，我们通常使用这些模型来描述现实世界系统的行为。我们事先并不知道该系统的确切行为，但正如我所说，我们知道系统背后的机制是如何工作的。

我特意将本讲座的目标表述为一个图模型，因为我们希望将已知的概率模型与这些机制性模型、机制性知识融合起来，以便在实践中获得机器学习的一些优势。这有时也被称为混合建模，因为我们有来自观测数据的经验知识，也有机制性知识，它们可以很好地结合在一起。

我们知道机制，并希望推断系统行为。然而，PDE通常有一组参数，我们通常事先并不知道这些参数。这里有几个例子：

  • 强度和分布：例如热源强度、电荷分布。
  • 材料参数：例如热量通过特定材料扩散的速度，或经典力学中的力。

我们通常不知道这些参数，但可以测量它们。这些测量通常是有噪声的。由此我们已经可以看到，我们的模型实际上需要能够处理观测噪声。而这些东西的经典描述通常处理得不太好。

我们拥有这些系统参数的测量值。为了求解系统行为，也需要知道这些系统参数，因为方程是由这些参数控制的。

有时，我们也有系统本身的测量值。例如，在热分布模拟中，我们可能在模拟的某个点放置一个温度计来测量温度。但我们可以在模拟的每个点都这样做。这就是为什么我们实际上需要机制性知识来在我们的测量点之间进行插值。

这取决于具体情况。有时可以直接测量力，有相应的测量设备。有时测量这些东西实际上更困难，在这种情况下，你可以使用这些测量值，并将知识反向传播到系统参数，这被称为反问题。Natan在ODE讲座中已经详细讨论过这个问题。所以，这实际上是我们在这里所做工作的另一面，即不仅仅是模拟系统做什么，还要推断其背后的因果机制是什么。

我们在这里的方法，实际上也是本讲座的目标，是使用贝叶斯统计估计来融合我们所知道的东西（例如我们知道的机制）、这些测量数据，以及可能包含在大多数系统参数中的不确定性。

为什么使用贝叶斯统计估计？因为我们有所有这些不确定性，并且实际上将我们不知道的所有东西传播到解中会是一件很好的事情，这样我们就可以为我们所做的预测提供一些置信度。

线性偏微分方程的定义

现在，让我们真正进入PDE的世界，首先回答一个问题：我们所说的线性PDE实际上是什么意思？

首先，考虑我们想要模拟一个物理系统。通常，物理系统具有某种空间范围。如果它是一个随时间演化的系统，我们还有一个想要模拟的时间跨度。首先，我们定义这个集合 D，它被称为定义域。然后我们寻找一个描述物理系统的函数，例如温度分布或一组电荷产生的力。这实际上就是未知函数 u，这是我们想要模拟的系统的描述。

在这些基于PDE的模型中，机制性知识实际上由这个方程给出：D u = f。其中 D 是一个所谓的线性微分算子，我稍后会展示一些例子。但一般来说，这只是未知函数 u 的偏导数的线性组合。我们不知道函数 u，但我们规定了一个固定值 f，称为方程的右侧项。事实证明，这是描述许多物理过程的一种非常优雅的方式。

现在，线性微分算子的一些例子：

  • 最著名的可能是拉普拉斯算子，它只是所有二阶偏导数的和，本质上是Hessian矩阵的迹。
  • 线性PDE的另一个例子实际上是一个一阶常微分方程。在函数只有一个输入变量而不是D个输入变量的情况下，我们可以构造这个微分算子，它只包含我们可以用这个函数构建的一个导数，以及另一个线性项。然后，如果你重新排列方程中的项，就会得到一个一阶常微分方程的向量场。从某种意义上说，一阶常微分方程是线性PDE的一个特例。一般来说，ODE是PDE的一个特例。这样思考并不总是有帮助，因为PDE往往更难模拟，但这只是一个闭合论证。

实际应用中的挑战

现在，谈谈如果我们要在实际中应用这些方程会遇到的一些问题。

首先，通常我们得不到解析解。这对于我们考虑的许多ODE来说已经如此。但在这里，从某种意义上说，因为它是一类更广泛的模型，我们实际上继承了这些问题。

因此，我们需要使用数值求解器来实际得到这里的未知函数。由于函数 u 实际上是一个无限维对象，这将固有地引入离散化误差，除非你对问题有所了解。我们稍后会讨论这个问题。

其次，我们已经讨论过，PDE有参数，我们现在可以更准确地指出这些参数是什么。右侧函数 f 是这些参数之一，例如，当我提到热源、电荷时，这些通常由方程的右侧描述。而材料参数等，则是函数 u 的偏导数线性组合中的系数。我们通常并不确切知道这些。

最后，在过去一个世纪中发展起来的经典求解器有时很难嵌入到计算流程中，原因在于参数通常并不确切已知，但我们需要一个点估计。因此，用这些求解器传播这些不确定性往往相当困难，而这正是我们在这里采取的贝叶斯方法旨在解决的问题。

边界条件与唯一性

这里还有一个问题，因为PDE本身通常不能唯一地确定其解。对于ODE，我们已经看到我们还需要表述初值问题，因为本质上涉及一个积分常数。对于PDE，情况没有太大不同，但我们需要施加的附加条件的类型比ODE情况要复杂一些。

但让我们先看一个例子。我们来看所谓的泊松方程，它只是函数的拉普拉斯算子等于一个规定值。现在，让我们考虑这个方程的一个解候选，它只是一个线性函数。我们可以看到，如果我们将拉普拉斯算子应用于这个线性函数，由于这是一个线性函数，Hessian矩阵是0，没有二阶项。所以这正好是0，这意味着线性函数实际上在这个微分算子的核中。因此，对于泊松方程的任何解，我们都可以直接添加一个线性项，仍然得到一个解。因为，我稍后会谈到，这个微分算子实际上是线性的。

因此，通常，也许本着你在ODE情况下所做的精神，唯一性可以通过要求一个额外的条件来实现，这通常是一个边界条件，或者在这种情况下，为了保持一切线性，是一个线性边界条件，其中我们再次有一个线性算子。我稍后会详细讨论这一点。

从物理直觉来看，为什么我们要规定模拟域边界上发生的事情？实际上，如果你考虑你正在模拟的问题，可能有来自模拟域外部的相互作用进入你的系统。如果你不模拟这些外部影响，即如果你没有在你的数学框架中对它们建模，那么本质上任何事情都可能发生。例如，如果你模拟一块金属的热分布，有一辆卡车冲进来，那么会发生什么？你不知道。因此，你本质上需要模拟模拟域外发生的一切，通过总结外部如何与你在边界上的模拟相互作用，因为这是你实际上可以施加影响的唯一方式。

一个PDE，连同一组边界条件，通常被称为一个边值问题。

现在，这种边界条件的一个具体例子是狄利克雷边界条件，它本质上规定函数在边界上的值。如果你考虑函数 u 是某块材料中的热分布，那么这只是说我们知道边界处的温度是多少。在物理类比中，这可能就像你有一个巨大的冰水桶，无论域内发生什么，它的温度都不会改变，所以你知道你的域边界将始终保持为0。

函数空间与线性算子

好的，我已经围绕这一点讨论了一段时间。但PDE实际上是关于函数的陈述。所以我们有一个未知函数 u，而函数通常是无限维对象。你可能已经知道，如果涉及无限，你需要小心一点。

但事实证明，函数具有相当方便的结构，因为事实证明，函数或函数空间（在特定条件下）实际上是向量空间。从线性代数课程中，你可能非常熟悉向量空间 R^N，本质上是N维列向量的空间。你知道，你可以将它们相加，可以用标量乘以它们。对于函数，你可以做本质上相同的事情，只需说两个函数的和以及函数与标量的乘积只是逐点定义的。这就是这里所设定的。

我们具有向量空间结构，但这也带来了向量空间的所有便利。在这些向量空间中可以有基。它只是一组基向量的线性组合。对于函数空间，你有一组基函数的线性组合。通常，至少对于我们这里关心的函数空间，这些空间实际上是无限维的，有时甚至不存在于这种形式中。这就是我说的需要小心的意思。但对于本讲座，这些空间存在，它们是无限维的，并且本质上作为一系列向量存在。

你也有线性映射，从有限维向量空间我们知道，我们可以用矩阵表示线性映射。线性性质意味着我们可以将矩阵拉入和式中，并且可以将标量从矩阵向量积中拉出来。

无限维向量空间，特别是函数空间中的线性映射，被称为线性算子。这实际上也是“线性微分算子”这个术语的由来。因为微分算子将一个函数映射到其导数。所以它是函数向量空间之间的一个映射，具有这种线性性质。本质上，因为你有求导的和法则，并且可以将标量从导数中拉出来。

这实际上就是我们谈论线性方程和线性微分算子的原因。通过应用这些知识，你可以看到线性PDE只不过是无限维向量空间中的一个线性系统。因此，在本讲座中，我们将实际应用从Unhan讲座中发展起来的求解线性系统的相当多的直觉来处理PDE的情况。这实际上是一个很好的类比，在处理这些系统时你总是可以想到。

还有两个细节。你可以在向量空间上定义范数。例如，这里我们有 R^N 上的最大范数或无穷范数。本质上，对于某些函数类，也存在这种范数的类比。你可以直接取函数的上确界，然后事实证明它是一个范数。如果我们有范数，我们也可以谈论收敛性。如果这些空间中的每个序列都收敛，我们可以谈论巴拿赫空间。我们知道带有无穷范数的 R^N 是一个巴拿赫空间，但实际上，某个集合上k次连续可微函数的空间，在这个特定范数下，实际上也是一个巴拿赫空间。所以我们再次可以将有限维情况下的许多直觉带到无限维情况下，但也许并不总是如此。所以我们必须再次在这些意义上小心。

实际上，对于内积和希尔伯特空间也是如此，但有一些注意事项，我现在不打算深入讨论。只需注意，通常这些转换涉及将列向量中的索引替换为函数求值，然后将求和替换为积分。这是一种非常直接的推导这些类比的方式。

一个实际建模示例：CPU热分布

好的，让我们转向一个更实际的话题，这实际上是一个玩具示例，将作为介绍我们将要使用的方法的主要动机示例。这是一个计算机中央处理器中热分布的简单模型。

现在，这是典型的CPU外观。实际上，顶部的这块金属片不是芯片本身，它只是一个盖子，用于从芯片中提取热量。实际的硅片只是这个小黑芯片。这些组件特别受其散发的热量限制。如果有电流通过，由于内部电阻等原因，芯片会产生大量热量。如果这些系统过热，它们可能会损坏或无法正常工作。我们希望避免这种情况。因此，在实践中了解芯片中的温度分布实际上是一件好事。

你可能已经猜到，这个东西实际上非常薄。因此，我们可以将其建模为一个二维物体。现在，我们实际上将自己限制在二维建模，即这个系统的数学建模。这是我们谈论的空间域。在这种情况下，它只是芯片长度和宽度的区间笛卡尔积。但实际上，由于这里有一些均匀的几何形状，为了本讲座的简单起见，我们甚至将自己限制在仅模拟这条穿过芯片的线上的温度分布。我们实际上将在讲座接近尾声时看到一个二维示例，并且我们将看到对于这个特定情况，这实际上是一个相当好的模型假设。

但现在，让我们转向从物理模型转向数学。首先，我们已经定义了我们将要处理的空间域，它只是从0到CPU长度的这个区间。现在，我们知道当我们想要模拟涉及热的问题时，我们实际上需要知道热量来自哪里，即什么在我们的系统中产生热量。现在，我们实际上假设CPU内置的GPU实际上处于空闲状态，它什么也不做。而计算核心实际上正在计算一些非常困难的东西。所以CPU的计算核心正在工作并产生热量。同时，这些系统通常以某种方式内置在计算机中，顶部有一个散热器，将芯片产生的所有热量带走，以防止其过热。

实际上，这个草图可能有点误导。散热器通过挤压所谓的导热界面材料（一种膏状物）来提取热量。在这些边缘周围，你可以假设它从芯片表面的所有部分均匀地提取热量，所以热量通过其表面离开芯片。

因此，如果我们实际上想要模拟芯片中的热源在哪里，以及散热器在哪里（热量离开系统的地方），我们可以看一个像这样的函数。我们在核心本身上放置了三个高斯斑点，它们是热源。你可以看到这里的单位是瓦特每立方毫米，所以本质上是每单位体积的热量。然后在这些高斯斑点上叠加了另一个负常数函数，它模拟了散热器，即所有从CPU中被拉出的热量。

现在你可能想知道，我们在这里谈论PDE。那么实际模拟这个系统的PDE是什么？我已经谈过热方程。它就在这里。它是一个线性PDE，也是二阶的。为什么是二阶？因为拉普拉斯项包含非混合二阶偏导数，正是如此。这就是二阶的由来。我稍后会谈谈这些单独项实际上意味着什么。但现在请注意，这个函数 u 将是我们芯片中的温度分布。

我们可以看到这里有一个时间变量的导数。为了使事情对我们来说更简单，一个合理的假设是假设温度分布随时间保持不变，这最终会在一个持续运行的CPU中发生，并且当散热器的风扇控制达到一个稳定点时。所以我们要做的是假设，在某个时刻，我们实际上将达到一个稳定的温度分布。你可以通过说，如果温度不变化，那么这个时间导数正好是0，温度没有变化，至少不再有温度变化，从而将其纳入这个模型。在这种情况下，我们得到了稳态热方程，你只需将整个第一项设置为0，就可以从这个方程得到它。

由于我们在这里将自己限制在CPU的这个一维子集，即这条线，这实际上变成了这个方程。在D=0的情况下，拉普拉斯算子只是一个二阶导数。现在，这是什么？我们实际上已经在本讲座中看到过这个方程。它让你想起了什么？它只是一个ODE。因为我们限制自己处理一维问题，它是一个ODE。你可能会说，在这里谈论PDE有什么意义？这只是为了让本讲座的内容更直观一些。我们接下来要讨论的一切实际上也适用于多维情况，我们在这里没有任何针对一维情况的专门化。

将机制性知识注入统计模型

现在，我们如何将从这个微分方程中获得的机制性知识注入到我们的统计模型中？

思考这些PDE的一种方式实际上是将其解释为对这个未知函数 u 的观测。这本质上与我们处理ODE情况的方式非常相似。有一个未知量。我们不知道关于这个函数 u 的任何信息，但我们可以观察到一个导出量。我们可以观察它在微分算子下的像。并且我们知道，因为我们希望PDE在我们的系统中成立，这个微分算子的像（暂时将其视为微分）必须是一个特定值，即右侧函数的值。

在物理学中，很多最基本的定律实际上都是守恒定律。它们描述了能量、质量、动量、电荷等一些基本物理量的守恒。事实证明，这些通常被表达为PDE。我们刚刚讨论的热方程（实际上我将一项移到了另一边）是一个关于能量守恒，特别是热能的陈述。这个左侧项与温度的变化（温度的时间变化）成正比。实际上，这也与热能的变化成正比。如果温度变化，热能变化。请注意，温度和热能实际上不是一回事，但在这种情况下，它们通过这些材料参数彼此成正比。我们说过，热能的每一次变化都必须通过热源（即 Q_dot_V 实际上是一个热源）来解释，本质上，已知的热量在任何点进入系统；或者必须通过热传导从周围环境流入某一点来解释。这实际上是有道理的，因为拉普拉斯算子计算了函数的某种曲率估计。在一维中，它实际上是曲率，就是二阶导数。如果你有一个类似抛物碗的情况，那么你会期望，因为中心点周围的温度比中心点本身更高，热量会流入那里，因为温度倾向于随时间达到平衡。这个陈述是，能量的每一次变化要么由来自周围环境的传导解释，要么由那里存在的热源解释。所以，除了我们能解释的之外，没有能量损失或增加。

这是一个局部陈述，因为涉及所有这些导数，它们是在该域的每个点计算的。

我已经说过，通常，或者抽象地说，这只是一个局部数学性质。所以导数在某个点的规定值。

我们实际上可以写出这个方程。首先，这只是一个符号上的事情。但我们可以将其写为微分算子减去右侧等于0。我们将其定义为函数 I 或算子 I，在概率数值方法中也被称为信息算子。对于仅仅是导数的情况，你实际上已经在ODE讲座中见过这些信息算子。这正是我们在概率ODE求解器中条件化的东西。所以这只是一个更一般的框架。它不仅适用于微分方程，而且你本质上可以用这样的信息算子表达任何信息。你可以将其视为数据点概念的扩展。PDE本身可以说不是一个数据点，但它仍然为你提供了关于你试图解决的问题的信息。所以它是数据或信息的广义概念。

使用高斯过程求解PDE

好的，现在我们实际上想要求解这个微分方程。我们有一个未知函数 u。那么，如果我们在贝叶斯统计中有一个未知量，我们该怎么做？我们只需在其上放置一个先验。在这种情况下，不出所料，它是一个高斯过程先验 u。你实际上可以在上面看到那个先验。它是一个Matérn 7/2核，以及一个常数均值函数。

而观测，它取代了常规普通高斯过程推断中的正常点观测，现在是信息算子。我们要求微分方程在域的每个点都成立，现在我们将解候选替换为高斯过程。

现在，我们实际上如何做到这一点？首先，我们如何将这样的微分算子应用于高斯过程？我们从那得到什么样的对象？换句话说，我们如何取高斯过程的导数？其次，我们现在如何实际计算后验？我们如何对这个信息进行条件化？

这就是，哦，实际上，事实证明这两个对象又是高斯过程。所以我们有一个在高斯过程下线性观测的某种闭合性质。

我们仍然有点问题，因为这实际上是一组无限的观测，因为我们希望PDE在域的每个点都成立。而域通常像一个区间，例如，它是不可数无限的。因此，在计算上，这将是一个相当大的挑战。对于一般情况，这基本上是不可能的，除非你对你要解决的问题有解析上的了解。

因此，我们通过本质上说，我们不希望它在域的每个点都成立，而只是在一组有限的训练点 X 上成立，来放松这个信息。这些点由于类似于此的经典方法而被称为配置点，我认为Natan在ODE讲座中也谈到过这个。所以这本质上是我们ODE讲座中采用的方法，只是我们这里不使用状态估计，而是使用实际的高斯过程。

现在，事实证明，这些对象非常类似于如果你对一个有限维高斯随机变量进行线性观测条件化所得到的形式。这实际上已经出现在关于高斯过程的讲座中，它本质上只是高斯推断定理在 R^D 上的应用。但现在回想一下，我们实际上说过，函数也是向量空间，所以我们可以将高斯过程视为函数空间上的概率测度，或者类似于函数空间上的随机变量。我们只是在这里使用某种高斯过程先验。然后和以前一样，我们选择一个线性算子。在这种情况下，它实际上是一个将高斯过程的样本路径映射到 R^N 的线性算子，即对该高斯过程的N个线性观测。我们在这里引入某种噪声变量，它独立于高斯过程，这只是独立的高斯噪声，就像有限维情况一样。

现在，事实证明，先验预测分布，即高斯过程在这个线性算子下的像，实际上由一个正态分布给出。其中均值由高斯过程均值函数通过该线性算子的像给出。协方差矩阵与我们实际在有限维情况下得到的协方差矩阵非常相似，只是在这里我们不能真正应用矩阵乘积，因为这是一个线性算子，而不是矩阵。而且这不是一个线性算子本身，也是一个函数。但在函数空间中，这个 A Σ A^T 类比的本质是首先将线性算子应用于核函数的第一个参数。所以固定核函数的第二个参数，然后你有一个关于 x 的单变量函数，它返回一个实数值。这实际上正是我们可以输入到那个微分算子的东西，因为，嗯，那是它定义的空间。然后之后，再次将其视为 x 的函数，现在它从空间域 X 映射到 R。这又是一个我们可以输入到微分算子的函数。所以我们现在再次应用它。你稍后会看到一个带有具体微分算子的具体例子，这将使一切更加清晰。

而这个事件的后验实际上也具有类似的结构。本质上，我们只是用将线性算子应用于函数来替换矩阵向量乘法。我们用我们已经在这里看到的那个东西来替换这个 A Σ A^T 对象，即那个基础矩阵。我在这里省略了很多理论细节。我们将在讲座结束时至少以某种模糊的方式回到这一点，但要知道，涉及的内容远不止写下这些方程。由于涉及所有这些无限性，在做这件事时你需要非常小心。

应用示例：配置点与边界条件

现在，如果我们实际上想要计算这个对象，它看起来是什么样子？对于这个线性算子的具体选择，有人有想法吗？它只是某个固定的 x，比如3。所以它只是，你知道，你对一个参数求导，然后你对另一个参数求导。就这么简单。只是用这种标准线性算子符号表达起来有点困难。所以首先，将第二个参数固定为某个值，然后将线性算子应用于核函数（现在它只是第一个参数的函数），意味着我对第一个参数求导一次，然后将 x 代入其中。然后现在我将其视为第二个参数的函数，然后我对其求导。对于微分算子来说，这实际上很容易做到。然后，如果你这里有多个 x 点，那么你只需构建所有这些点之间两两导数的矩阵，本质上是为了得到这个对象。它实际上是一个矩阵。

实际上，为了进入你刚才提到的情况，即我们实际上不进行点求值，而是有一个函数空间之间的适当线性算子，我们也定义这些协方差核。这是一个核函数，这些是互协方差函数。它们本质上与我们上面定义的相同，只是我们实际上将其视为我们之后想要对其求导的点的函数。你可以证明，通过应用你刚刚在前一张幻灯片上学到的知识，如果你有这样一个算子，并且它满足某些条件（我们稍后会讨论），那么该算子下的高斯过程的像就又是一个高斯过程。本质上，根据定义。这就是你如何得到高斯过程的导数，例如，如果这只是一个导数算子，例如对于标量高斯过程，只求导一次。然后你对均值求导，并对两边的核函数求导，本质上，以得到一个对称函数。例如，你可以这样看它。然后那又是一个高斯过程，在某些条件下。

现在，让我们将其应用到我们的PDE案例中。所以我们回到同样的高斯过程先验和观测。我们已经说过，我们已经看到将微分算子应用于高斯过程又是一个高斯过程，具有特定的后验矩形式。你实际上可以在这里看到那个高斯过程。所以在这种情况下，是那个二阶导数或缩放后的二阶导数，因为那是我们正在处理的微分算子。你可以看到这种情况的发生，因为样本实际上不那么平滑。通过对这些函数求导，我们失去了可微性的阶数。这也是你知道你正在使用Matérn核而不是像平方指数核这样的东西的原因，因为平方指数核的样本实际上是平滑的，所以是无限次可微的。

现在，既然我们已经看到了这一点，我们也可以将相同的定理或相同的表格应用于这个问题，其中我们的线性算子现在是，首先应用 D，然后在点集处求值。我们可以计算后验。在这种情况下，我们知道它现在是一个高斯过程。如果我们实际将其应用于这个问题，那么我们首先定义一组配置点。这里的蓝色虚线就是这些 x 点。然后，观测由这个黑色函数给出，因为本质上这是PDE的右侧空间，一旦我们应用了微分算子，高斯过程必须与右侧函数匹配。这就是为什么这些观测在这个变换空间中的 y 值只是由右侧函数的点求值给出。你可以看到，在这个空间中，我们本质上只有一个普通的高斯过程回归问题。但我们实际上将我们从那里获得的知识传播回这个原始的高斯过程，它通过微分算子连接。我们可以看到它对此有所反应，但似乎并没有真正起作用。仍然有很多不确定性。那么问题出在哪里？

想想我说的关于PDE不能唯一求解的问题。边界条件缺失了。

它实际上确实有效。从某种意义上说，它没有坏。你可以看到所有单独的样本，你也可以从这个图中看到，所有样本都近似地求解这个方程。只是存在PDE没有固定的自由度。这实际上正是我所说的线性自由度。所以这个高斯过程的每个样本，至少近似地，彼此之间相差一个添加的线性函数。只是，你知道，只是缩放或倾斜，实际上正是通过这一项。

既然这个后验现在又是一个高斯过程。如果我们，例如，施加狄利克雷边界条件，即我们规定这个温度分布在区间边界点（左和右边界点）的值，从物理上你可以将其解释为你可能在那里连接的温度计的测量值。那么，我们只有一个普通的高斯过程回归问题，因为前一个问题的后验是一个高斯过程，你可以直接将其作为普通高斯过程回归问题的先验。你只是在边界处观测两个值，这没有什么特别的。

所以，我们这里仍然有一些不确定性。实际上观察到右侧基本上没有变化，所以本质上和以前一样。但现在不确定性崩溃了。你在这里看到的剩余不确定性只是近似误差。我们没有实际要求PDE在域的每个点都成立，只是在这些配置点成立。但这实际上也反映在高斯过程的置信度估计中。它好像在说，我不确定，因为我没有获得所有信息。

模型改进：更现实的边界条件与参数不确定性

所以让我们稍微回顾一下。我们做了什么？我们已经看到，广义形式的高斯过程推断实际上可以产生我们在这里表述的这个边值问题的近似解。并且我们得到了一个通常很难处理的近似误差的估计。

现在，我们实际上回到了开头的图模型，如果你还记得的话。所以我们有一个物理系统的描述，我们最初并不知道，但我们对其条件化，基于我们拥有的一些关于系统的机制性知识的观测。在这种情况下，我们测量CPU在边界处的温度值，并且我们知道热量如何通过一块硅片流动。

现在，关于这一点有点不现实的是，边界值在部署中通常是未知的。如果你的系统中有CPU，这些温度计并不存在。CPU上有温度计，但不在边界处。所以我们需要摆脱这一点，以便实际拥有一个现实的模型。其次，这些热源，这些热源的值也不确切知道。

我们建模的方式是，好吧，如果核心实际上在计算什么，那么让我们用高斯斑点来近似热源分布。但你不知道它实际上是一个高斯斑点。我们没有测量过。

为了解决这个问题，实际上很酷的一点是能够对边界值和这个热源分布的确切值都添加不确定性，因为这些也不是确切已知的。

好的，让我们回到我们只对PDE进行条件化的情况。对于这个情况，什么是更现实的边界条件？至少是更实用的边界条件。事实证明，我已经说过，我们知道热量大致均匀地从CPU的整个表面提取，实际上这些边界点就是这个盒子的侧面部分，这个盒子是CPU的理想化表示。所以从物理上，你实际上可以对这些边界条件建模。所以关于有多少热量从表面提取的信息转化为所谓的诺伊曼边界条件，它不是设置边界点的值，而是说边界点的一阶导数。如果有大量热量被提取，那么你可以假设温度分布在边界处以一个相对陡峭的下降斜率移动。如果有大量热量进入，那么它是一个正斜率，至少在右侧斜率上是这样。

所以如果我们实际这样做，我们对该导数的未知值进行建模，因为我们不知道确切有多少热量通过边界离开，但我们可以添加一些正负不确定性。

所以我们用另一个高斯分布来建模，在这种情况下是一个高斯过程。但实际上，它只是一个高斯分布，因为它只是两个值而不是一个函数。我们添加一个额外的信息算子，它由边界处的那个导数给出。不要担心这个。这是一个方向导数，因为通常，如果你有一个多维域，例如像一个平面，它不仅仅是导数，而是边界法线方向上的导数。所以本质上，有多少热量通过那个法线方向流出表面。但现在就把它想成一阶导数。

这是一个描述诺伊曼边界条件的信息算子。在这种情况下，它也不再是一个普通的高斯过程回归问题，因为我们有通过另一个线性算子的观测，它只是与PDE的微分算子不同的一个算子，是我们之前看到的这个边界算子 B。

如果我们实际这样做，哦，对不起，我们看到我们摆脱了一个自由度。所以不是存在一个线性自由度，我们可以有一个额外的斜率，我们固定了斜率，但这个高斯过程样本中唯一剩余的自由度是偏移量，即平移自由度。这有点道理，因为我们只是说了热量流向哪里和从哪里流出，而不是系统在什么绝对尺度上运行。

所以这可能是一个在1000摄氏度下运行的系统，也可能是在零下200度。我们实际固定那个绝对尺度的方式是通过CPU上实际包含的温度计，只是不在边界上，而是在CPU核心的中心。所以我们有来自这些传感器的测量值，热读数，这些传感器被称为数字热传感器，沿着这条穿过CPU的线有三个点。我们实际上添加了这些。这使不确定性崩溃，但没有完全崩溃，因为这些温度计存在测量不确定性。所以我们只能做这么多。本质上，我们可以从这三个测量的绝对尺度中学到什么。

现在，我们再次说过，我们实际上并不知道热源项。应该对它有一些不确定性。这只是对可能发生的情况的粗略估计。那么我们如何解决这个问题？也许，甚至很琐碎，只需添加另一个高斯过程。所有东西都是高斯过程。

所以与其说这个PDE的右侧是某个固定函数，我们不如说，好吧，它只是某个高斯过程。它是一个概率测度，一个对我们实际上应该知道的函数的不确定估计。所以我们用另一个高斯过程来建模我们对该函数实际是什么的先验信念。我们可以使用相同的推断技术。但这里有点问题，因为问题的物理实际上，首先，我们之前使用的这个确定性右侧函数是精心选择的，使其实际上积分为0。这是有道理的，因为如果它不积分为0，那么总体上，进入系统的热量将多于离开的热量。在这种情况下，我们将永远达不到热平衡。系统只会持续加热，因为总是有更多的能量进入系统。所以我们需要选择这个右侧函数使其积分为0。

但是你如何用高斯过程做到这一点？如果你只是向那个单一估计添加高斯噪声，那么其中一些样本将持续位于我们用作右侧函数原始估计的均值之上。那么，它们将具有更大的积分。那么我们如何解决这个问题？我们如何本质上保证这个高斯过程在其所有样本中面积为1？积分是一个线性算子。积分是线性的，因为两个函数和的积分是两个函数积分的和。

所以我们实际上可以表述另一个线性观测，即线性信息算子，它正好表述了这个特定条件，我称之为稳态条件，因为它指的是热稳态。实际上，你还必须加入边界效应，因为虽然有热量通过边界离开，但通过边界离开的热量，加上内部离开的热量（由右侧函数中的那个负项建模），再加上CPU核心产生的所有热量，需要为0，本质上需要平衡。

因此，通过实际添加这个额外的约束（它现在实际上并不首先影响我们的函数 u，它只是右侧函数和边界函数 Q_V 和 Q_A），对这些函数的先验发生了变化。如果你仔细观察，你可以看到，开始时低于均值的样本（均值实际上已经积分为0），实际上最终会持续高于它。然后它们以这样的方式发生，使得这个积分为零的条件对所有样本路径都成立。

现在，我们可以将这些信息算子应用于这个联合高斯过程先验。这实际上现在是一个多输出高斯过程先验，以得到这个系统。

我认为，就其对问题所做的假设而言，这是一个比我们开始时更现实的模型。

你可以看到，仍然有一些不确定性，这是有道理的，因为我们没有确定的右侧函数，没有确定的边界条件，也没有确定的测量值。所以所有这些不确定性都贡献给了解的后验不确定性。但你可以看到，红色阴影区域，实际上，我意识到很难分辨，这里的红色阴影区域本质上持续位于这个蓝色阴影区域内。蓝色阴影区域是关于右侧函数的不确定性。所以本质上，我们的大多数估计，或者说我们的大多数估计，对于高斯过程在微分算子下的像的曲率，与关于右侧的不确定性是一致的。

扩展到多维与时间相关问题

所以这很有道理。好的，让我们再次回顾一下。我们已经看到，这种高斯过程方法，结合信息算子的概念，使我们能够整合关于系统行为的先验知识（表述为高斯过程，例如Matérn核所做的假设），注入机制性知识（以我们表述的这个线性PDE的形式），处理不确定的边界条件和右侧（这在实际中经常出现，特别是在反问题的背景下），并使用有噪声的经验测量来消除我们从不确定边界条件、右侧和近似不确定性中得到的一些不确定性。所有这些都在我们提供近似误差的量化时发生，这个误差来自于我们没有要求PDE在域的每个点都精确成立。我们得到了从这些系统参数（右侧、边界条件）的不确定估计中传播而来的误差。

这背后的更大故事是，所有这一切之所以可能，是因为我们不仅仅给出解的点估计，而是实际上稍微放松了这一点，说我们只是用一个由概率测度加权的无限解候选集来回答。本质上，这只是贝叶斯推断承诺的另一种形式。但它非常强大。与其使用该函数的点估计，不如使用函数上的概率测度，它可以给你置信度，可以给你样本，所有这些样本在这种情况下都精确满足你施加的条件。

而且它可能也更诚实，因为与其说我们实际上知道这是我们想要解PDE的右侧函数，你承认这些估计中几乎总是存在不确定性。你几乎从不知道确切的值。这是实际建模这一点的一种方式。

现在，我已经看到，你实际上也可以在二维中模拟完全相同的模型。这种方法完全相同，只是你将一维点网格替换为二维点网格，并稍微改变先验。你可以看到它有效。实际上，因为这里沿Y轴的温度几乎相同，你可以认为一维模型实际上一开始就是一个相当好的模型，因为没有太多变化。通过以我们所做的方式设定这个一维模型，我们说沿那个维度的温度本质上是相同的，所以可能这样建模是可以的，只是作为这里的一个概念证明。

现在，与ODE讲座相比，我们在这里并没有真正讨论时间。但时间除了只是另一个空间维度之外，还有什么？你可以争论，时空域，本质上。你可以将完全相同的方法应用于热方程的一维版本。所以现在是实际的热方程，不是我们将时间导数设置为0的那个，现在这个轴是时间轴。正如我所说，它只是另一个空间维度。这是我们的一个空间变量。现在你可以看到，从某种意义上说，这很像一个具有无限维状态空间的ODE。所以你可以将状态变量视为每个空间域点的状态变量。然后你有一个ODE部分，它描述了这种随时间演化。

如果你切片穿过这个函数随时间的变化，并实际动画显示发生了什么，那么你会得到一个热扩散的物理上合理的模拟，它也具有所有这些不确定性，因为它是一个高斯过程。这里没有太多不确定性，因为我没有向右侧和边界条件添加不确定性，并且我使用了一个相对密集的网格。但最终，你可以看到网格指定了，特别是在边缘有更多不确定性。

所以你也可以通过本质上说时间没有什么特别的，时间只是另一个输入维度，来用这种方法模拟时间相关问题。

与经典数值方法的联系

好的，所以我谈了一点关于这些PDE的经典数值方法，这些方法在过去一个世纪中发展起来。那么这种方法如何融入其中？当所有这些东西都已经发展起来时，使用高斯过程可能看起来有点奇怪。但事实证明，我们一直在发展的这种方法的后验均值，假设你不添加任何不确定性（即你有精确的边界条件和精确的PDE右侧），正是某些经典方法产生的点估计。在这种情况下，它被称为对称配置法，这也是为什么这些点被称为配置点。

所以这只是一种方法。那么其他众多方法在哪里？更一般地说，我们实际上已经表明（不是在本讲座中），但你可以实际证明，所有所谓的加权残差法都可以实现为高斯过程的后验均值，只是没有不确定性量化，显然，而一些配置法实际上是这些方法的一个实例。但也有有限体积法，你可能听说过。还有一大类所谓的伽辽金或彼得罗夫-伽辽金方法，其中包含可能最著名的求解PDE的方法，即有限元法，以及谱方法。

你得到这些方法的方式本质上是通过改变你离散化方程的方式。记住，在这里我们取这个PDE残差 D u - f，然后只在几个点处求值。但你为什么只要求值？你本质上可以投影到任何其他函数上，因为你通常有一个希尔伯特空间。所以你有一个函数上的内积，你可以直接在那里放入另一个函数。通过以这种方式实现它，通过仔细选择你投影的函数和你工作的先验，你可以得到所有这些其他不同的方法。

通过一个有点不同的方案，你实际上也可以证明PDE的有限差分离散化可以通过高斯过程推断实现。这是我们小组的另一篇论文。

现在你可能想知道，为什么你甚至要在这一系列其他方法中使用高斯过程，或者除此之外？我们得到了

009：蒙特卡洛方法

概述

在本节课中，我们将学习一种称为蒙特卡洛的数值积分方法。我们将从它的基本思想、历史起源开始，探讨其数学原理，并介绍其核心算法——马尔可夫链蒙特卡洛及其高级变体。我们还将反思蒙特卡洛方法在计算中的本质及其局限性。

历史起源与基本思想

上一节我们介绍了数值模拟的复杂性。本节中，我们来看看如何通过模拟来解决积分问题，这源于一个历史性的需求。

蒙特卡洛方法起源于第二次世界大战期间，由斯坦尼斯瓦夫·乌拉姆和约翰·冯·诺依曼等人在曼哈顿计划中提出。他们面临一个核心问题：如何计算核裂变链式反应中中子的预期密度？这是一个积分问题。

他们提出了一个开创性的想法：用求和来近似积分。具体来说，如果我们想计算某个函数 φ(x) 在分布 P(x) 下的期望值，即积分 ∫ φ(x) P(x) dx，我们可以通过从 P(x) 中独立抽取 S 个样本 x_s，然后用样本均值来近似：

φ̂ = (1/S) * Σ_{s=1}^{S} φ(x_s)

当时，他们甚至没有现代计算机，而是使用了一种叫做“费米模拟计算机”的模拟装置，配合骰子来生成随机数，手动模拟中子的随机路径，从而估算积分。

简单蒙特卡洛：原理与性质

从历史回到数学，我们来分析这个简单蒙特卡洛估计器的性质。

无偏性：该估计器是无偏的。因为样本 x_s 独立同分布于 P，所以估计器的期望值等于真实积分值：
E[φ̂] = (1/S) * Σ_{s=1}^{S} E[φ(x_s)] = ∫ φ(x) P(x) dx = φ

收敛速率：估计器的方差决定了其误差。可以证明，其均方误差为：
Var(φ̂) = Var(φ(x)) / S
这意味着误差以 O(1/√S) 的速率衰减。从统计学角度看，这是无偏估计器可能达到的最优速率。但从数值分析角度看，这个 1/√S 的收敛速度相对较慢。

维度影响：一个常被提及的观点是，蒙特卡洛的误差常数 Var(φ(x)) 表面上与问题维度 D 无关。然而，对于某些精心构造的高维函数，这个方差可能会随维度急剧增大（例如指数增长）。但在许多实际应用中，例如计算均值或方差时，蒙特卡洛方法对高维问题依然相对稳健，这是其相对于下周将介绍的方法的一个关键优势。

蒙特卡洛实践：示例与局限

为了直观理解，我们来看一个经典示例：用蒙特卡洛方法估算圆周率 π。

方法是：在单位正方形 [0,1]×[0,1] 内均匀随机采样点 (x, y)。判断点是否落在单位圆扇形（x² + y² ≤ 1）内。落在扇形内的概率是 π/4。因此，用频率估计概率，再乘以4即可得到 π 的估计值。

以下是核心代码逻辑：

# 伪代码示例
samples = uniform_random(size=(N, 2))
inside = sum(samples[:,0]**2 + samples[:,1]**2 <= 1)
pi_estimate = 4 * inside / N

绘制误差随样本数 N 变化的曲线，可以清晰地看到误差以 1/√N 的速率下降。

这个例子揭示了蒙特卡洛方法的本质：

  • 优点：对于难以直接计算的积分，它能提供一种相对容易实现的近似方案，快速给出一个粗略估计。
  • 局限：若要获得高精度结果，则需要极其大量的样本，计算代价高昂。

因此，在机器学习中进行贝叶斯推断时，如果只想要一个下午就能得到的粗略答案，蒙特卡洛是个好主意。但如果需要精确结果，则需要考虑其他方法。

从简单蒙特卡洛到马尔可夫链蒙特卡洛

简单蒙特卡洛的前提是能够直接从目标分布 P(x) 中采样。然而在实践中，P(x) 往往非常复杂，无法直接采样。因此，我们需要构建一个算法，来生成近似服从 P(x) 分布的样本。

标准方法是马尔可夫链蒙特卡洛。它通过构建一个马尔可夫链，使其平稳分布恰好是我们的目标分布 P(x)。当链运行足够长时间后，其状态就可以被视为来自 P(x) 的近似样本。

最著名的MCMC算法是Metropolis-Hastings算法。其步骤如下，我们用一个简短介绍来引入：

以下是Metropolis-Hastings算法的单步流程：

  1. 给定当前状态 x_t。
  2. 从某个提议分布 Q(x' | x_t)（例如以 x_t 为中心的高斯分布）中抽取一个候选点 x'。
  3. 计算接受概率 α = min(1, [P(x')Q(x_t | x')] / [P(x_t)Q(x' | x_t)])。若 Q 对称（如高斯），则简化为 α = min(1, P(x')/P(x_t))。
  4. 以概率 α 接受候选点，令 x_{t+1} = x'；否则拒绝，令 x_{t+1} = x_t。

该算法满足细致平衡条件：P(x_t) T(x_{t+1} | x_t) = P(x_{t+1}) T(x_t | x_{t+1})，这保证了 P(x) 是其平稳分布。在一定的遍历性条件下，马尔可夫链最终会收敛到该平稳分布。

然而，MCMC的收敛速度可能很慢。对于病态（长宽比大）的分布，链可能需要 O((L/ε)²) 量级的步数才能有效遍历分布，其中 L 和 ε 分别是分布最长和最短的长度尺度（类似于协方差矩阵的最大和最小特征值）。这导致有效样本量很低，估计方差很大。

进阶MCMC方法：哈密顿蒙特卡洛

为了克服随机游走行为导致的慢速混合问题，人们发展了更先进的MCMC算法。其中，哈密顿蒙特卡洛（HMC）是一个典范，它巧妙地利用了物理模拟和ODE求解器的思想。

HMC的核心思想是引入辅助变量——动量 p，将采样问题转化为在扩展的(x, p)相空间中模拟一个哈密顿动力学系统。假设目标分布为 P(x) ∝ exp(-E(x))（E(x) 称为势能）。我们定义哈密顿量 H(x, p) = E(x) + K(p)，其中动能通常取 K(p) = pᵀp / 2。相应的哈密顿动力学方程为：
dx/dt = ∂H/∂p = p
dp/dt = -∂H/∂x = -∇E(x)

这个系统有两个关键性质：

  1. 哈密顿量守恒：H(x, p) 在演化过程中保持不变。
  2. 相空间分布：联合分布 π(x, p) ∝ exp(-H(x, p)) = exp(-E(x)) exp(-K(p)) 恰好可分离为目标分布 P(x) 和一个标准高斯分布 N(p|0, I)。因此，从 π(x, p) 中得到的 x 边际分布就是 P(x)。

HMC算法交替进行以下步骤：

  1. 动量重采样：从 N(0, I) 中抽取新的动量 p。
  2. 动力学模拟：从当前 (x, p) 出发，沿哈密顿方程模拟一段轨迹（使用如蛙跳法等辛积分器），得到提议的 (x*, p*)。由于模拟误差，H 可能不严格守恒。
  3. Metropolis接受步骤：以概率 min(1, exp(-H(x*, p*) + H(x, p))) 接受新状态 (x*, p*)。

由于辛积分器能近似保持 H，接受率通常很高。更重要的是，动力学模拟允许提案在梯度方向上移动很远的距离，而不是随机游走，从而能更有效地探索分布，尤其对于具有强相关性的高维分布。其探索距离与模拟步数 L 成线性关系 O(Lε)，优于随机游走的 O(√L ε)。

HMC的性能依赖于两个关键参数：积分步长 ε 和模拟步数 L。设置不当会导致低效的“U型转弯”轨迹。No-U-Turn Sampler (NUTS) 是一种自适应调整 L 的HMC变体，它通过同时向前和向后运行轨迹来自动确定合适的模拟长度。

反思：蒙特卡洛中的“随机性”与算法选择

在介绍了主要算法后，让我们退一步思考蒙特卡洛方法的哲学基础。

蒙特卡洛方法的核心理论保证（如无偏性、收敛性）都建立在所使用的数是“真随机”的这一假设上。然而，在实际计算中，我们使用的是伪随机数生成器（PRNG）。PRNG产生的是确定性的、可重复的序列，只是看起来随机。一旦随机数种子被固定或公开，整个计算过程就完全是确定性的。

那么，我们为何还要使用基于“随机性”的论证呢？这引出了一个关键点：蒙特卡洛的“无偏性”是针对我们注入的随机性而言的。我们用一个确定性的PRNG来近似理想的随机性，并希望其性质足够好，使得理论论证大致成立。

这带来了蒙特卡洛方法，特别是MCMC的利弊权衡：

  • 优点：易于使用和实现；对于复杂的、高维的分布，通常能给出一个可运行的基准方案；理论框架成熟。
  • 缺点：收敛速度慢（O(1/√N)）；实际应用中难以诊断是否已收敛；其理论保证依赖于理想化的随机性假设。

因此，选择蒙特卡洛方法，通常意味着我们接受较慢的计算速度，以换取实现的简便性，将复杂的数学工作转化为计算机的耗时运行。对于许多问题，这是合理的权衡。但是，对于需要高精度结果或计算资源受限的情况，我们需要考虑其他数值积分方法，例如下周将要学习的那些具有多项式甚至指数级更快收敛速度的确定性或准蒙特卡洛方法。这些方法通常需要更多数学上的努力来构造，但一旦成功，计算效率会高得多。

总结

本节课中，我们一起学习了蒙特卡洛数值积分方法。

  1. 我们从其历史起源——曼哈顿计划中的物理模拟——开始，理解了用样本均值近似积分期望的基本思想。
  2. 我们分析了简单蒙特卡洛估计器的无偏性和 O(1/√N) 的收敛速率，并讨论了其对问题维度的敏感性。
  3. 针对无法直接采样的复杂分布，我们引入了马尔可夫链蒙特卡洛（MCMC）框架，详细讲解了经典的 Metropolis-Hastings算法 及其收敛挑战。
  4. 为了提升效率，我们介绍了哈密顿蒙特卡洛（HMC），它通过引入物理动力学模拟和ODE求解，实现了更高效的探索，并简要提到了其自适应变体NUTS。
  5. 最后，我们对蒙特卡洛方法中“随机性”的本质进行了反思，指出其理论保证与伪随机数实践之间的张力，并总结了蒙特卡洛方法易用但收敛慢的特点，为学习更高效的数值积分方法做好了铺垫。

蒙特卡洛是贝叶斯推断等机器学习任务中强大的工具，但理解其局限性能帮助我们在面对积分问题时，做出更明智的算法选择。

010：贝叶斯积分

概述

在本节课中，我们将学习如何将积分问题视为一个贝叶斯推断问题。我们将从一个具体的积分例子出发，探讨如何通过高斯过程先验来构建积分算法，并理解经典数值积分方法（如梯形法则）如何作为特定先验下的贝叶斯后验均值估计出现。我们还将讨论如何通过信息论准则选择评估点，以及如何通过引入更多先验知识来获得更高效、更自适应的算法。

问题引入：一个具体的积分

我们考虑一个具体的积分问题。函数 f(x) 定义为：

f(x) = exp(-sin²(3x) - x²)

这是一个实值函数，定义在实数域上。其图像是一个被高斯钟形曲线包络的振荡函数。

我们关心的量是定积分：

F = ∫₋₃³ f(x) dx

虽然函数 f(x) 被唯一确定，但其积分 F 的值对我们而言是未知的。我们无法从积分表中找到它，也没有现成的库函数可以直接计算它。因此，我们对 F 存在不确定性。

然而，我们并非一无所知。我们可以推断出 F 的一些性质：

  • 函数 f(x) 非负，因此 F ≥ 0。
  • 函数 f(x) ≤ 1，积分区间长度为6，因此 F ≤ 6。
  • 函数 f(x) 被高斯函数 exp(-x²) 上界约束，后者的积分（从-∞到∞）是 √π。在有限区间 [-3, 3] 上，高斯函数的积分是 √π * erf(3) ≈ 1.77。因此，我们有更紧的界：0 ≤ F ≤ 1.77。

这表明，即使对于一个未知的计算结果，我们通常也能获得一些先验信息。这启发我们可以将计算积分视为一个贝叶斯推断问题：为未知量（积分 F）和可计算量（函数值 f(x)）设置先验，然后进行贝叶斯推断。

回顾：蒙特卡洛方法

在深入贝叶斯观点之前，我们先回顾上节课讨论的蒙特卡洛方法。对于我们的积分问题，可以使用重要性采样。

方法：从一个已知解析积分的提议分布（例如高斯分布或均匀分布）中抽取样本 x_i，然后计算估计值：
F_MC ≈ G * (1/N) Σ [f(x_i) / g(x_i)]
其中 G 是提议分布 g(x) 在积分区间上的积分。

蒙特卡洛估计的误差以 O(1/√N) 的速率收敛，这是无偏估计器的最佳可能速率。然而，这个速率相对较慢。本节课的目标是探索能否显著超越这个速率。

贝叶斯积分：高斯过程框架

现在，让我们忘记传统的积分知识，尝试将积分构建为一个贝叶斯推断问题。我们将可计算的函数值 f(x) 视为数据，将待求的积分 F 视为潜变量。

我们首先在函数 f 上放置一个高斯过程先验。高斯过程由均值函数 m(x) 和协方差函数（核函数）k(x, x') 定义。

高斯过程有一个优美的性质：作用于该函数的任何线性算子，都会诱导出另一个高斯分布。积分正是一个线性算子。因此，如果 f 服从高斯过程，那么积分 F 也服从一个高斯分布（如果考虑不定积分，则是一个高斯过程）。

假设我们已经在一系列点 X = [x₁, ..., x_N] 上观测到了函数值 y = f(X)。那么，高斯过程的后验均值和协方差为：
μ(x) = m(x) + k(x, X) k(X, X)⁻¹ (y - m(X))
σ²(x, x') = k(x, x') - k(x, X) k(X, X)⁻¹ k(X, x')

我们对积分 F 的后验分布感兴趣。根据线性性质，F 的后验均值是后验均值函数的积分，后验方差是核函数的双重积分：
E[F|y] = ∫ μ(x) dx
Var[F|y] = ∬ σ²(x, x') dx dx'

注意，这些表达式只涉及对核函数 k 的积分，而不涉及对未知函数 f 的积分。因此，只要我们能够解析地计算所选核函数的积分，就能进行贝叶斯积分推断。

选择核函数：从维纳过程到梯形法则

我们需要选择一个能进行解析积分的核函数。一个简单的选择是维纳过程核（也称为线性核或最小核）：
k(x, x') = min(x, x') - c
其中常数 c 用于确保正定性，可以设为小于积分区间左端点的值（例如，对于区间 [-3, 3]，设 c = -3）。同时，我们设先验均值 m(x) = 0。

对于这个核函数，我们可以解析地计算其“核均值嵌入”和双重积分。具体公式如下（a 和 b 为积分上下限）：
∫ₐᵇ k(x, x') dx = 1/2 (x² - a²) + (x - b - a)c
∬ₐᵇ k(x, x') dx dx' = (b-a)³/12

现在，我们可以实现这个贝叶斯积分算法了。当我们在一系列点上评估函数并更新后验时，会发现一些有趣的现象。

核心发现：使用维纳过程先验和均匀网格评估点，积分 F 的贝叶斯后验均值估计，恰好等于经典的梯形法则的估计结果。

梯形法则是用分段线性函数来近似被积函数，并计算这些梯形面积之和。我们的贝叶斯框架不仅复现了这个经典算法，还额外提供了一个后验方差，作为积分估计的不确定性度量。

算法实现与计算效率：卡尔曼滤波视角

你可能会担心高斯过程推断的计算成本。然而，对于维纳过程这类特定的高斯过程，其推断可以通过卡尔曼滤波高效实现，计算复杂度与蒙特卡洛方法同阶。

我们可以将状态空间定义为包含函数值 f 和积分值 F 的向量。维纳过程对应于一个简单的随机微分方程。将其离散化后，我们得到卡尔曼滤波的预测模型和观测模型。

更新方程：在每一步，我们评估函数得到一个观测值 y_i，然后更新积分估计：
F_new = F_old + (Δx/2) * (y_i + y_{i-1})
这正是梯形法则的递推形式！同时，我们也能更新积分估计的不确定性（方差）。

因此，整个贝叶斯积分算法可以像蒙特卡洛一样高效运行，只是将随机求和替换为按特定规则（梯形法则）的求和，并附带不确定性估计。

更广泛的启示：经典算法作为贝叶斯估计

这个例子揭示了一个更深层次的模式：许多经典数值算法都可以解释为某种贝叶斯后验估计。

  • 线性代数求解器（如共轭梯度法）可以视为关于解或矩阵的最小二乘估计（高斯估计的特例）。
  • 优化方法（如拟牛顿法）可以视为关于损失函数海森矩阵的高斯估计。
  • 常微分方程求解器（如龙格-库塔法）可以视为单步内ODE解的特定高斯后验估计。
  • 偏微分方程求解器（如广义伽辽金法）也可以在高斯过程推断的框架下理解。

这意味着，经典数值计算在某种意义上已经是“贝叶斯式”的，它提供了点估计（后验均值）。而概率化数值计算的观点则更进一步，它免费提供了这些估计的不确定性（后验方差），因为我们计算后验均值时已经得到了所需的所有量。

评估点的选择：信息论准则

在之前的演示中，我们使用了均匀网格。为什么？这可以通过信息论来证明。

我们希望后验分布的不确定性（熵）尽快减小。对于高斯分布，熵由协方差矩阵的行列式决定。在我们的问题中，积分 F 的后验方差是：
Var[F] = (θ²/12) Σ (Δx_i)³
其中 θ 是核的尺度参数，Δx_i 是评估点之间的步长。

如果我们要求评估点必须覆盖积分区间的两端，那么通过最小化后验方差（即最小化熵），可以推导出最优策略是让所有内部的 Δx_i 相等，即使用均匀网格。

这展示了我们可以将计算过程视为一个智能体的决策过程：它遵循一个最大化信息增益的策略来选择下一步计算哪里。在这个简单例子中，策略的结果是均匀网格，但在更复杂或更高维的问题中，策略可能产生非平凡的自适应采样方案。

不确定性校准与先验假设的局限性

我们的贝叶斯积分算法提供了一个误差估计（后验标准差）。然而，在之前的图中，这个估计误差（红色虚线）远大于真实的积分误差（黑色点），并且以更慢的速率 (O(1/N)) 收敛，而真实误差以 O(1/N²) 收敛。

问题出在哪里？ 问题在于我们的先验假设。我们使用的维纳过程先验，其样本路径是几乎处处不可导的、非常粗糙的函数。我们告诉算法：“你要积分的函数可能极其不规则。” 因此，算法变得非常保守，给出了很大的不确定性估计。

然而，我们实际要积分的函数 f(x) = exp(-sin²(3x) - x²) 是无限次可导的，非常光滑。算法在每次评估时都发现函数比它预期的要光滑得多，但它的先验在每个尺度上都假设了粗糙性，因此误差估计始终过于悲观。

启示：如果一个算法被设计用于处理非常广泛的问题类别（如粗糙函数），那么它在任何一个具体问题（如光滑函数）上的表现通常不会是最优的。如果你有一个具体的、重要的计算问题，值得根据该问题的特定知识（如光滑性）来设计或选择算法。

改进：学习尺度参数与使用更合适的先验

我们可以通过共轭先验推断来在线估计核的尺度参数 θ。对于高斯分布，方差的共轭先验是伽马分布（作用于精度）。通过证据优化，我们可以更新 θ 的后验分布，而计算所需的关键量（yᵀ K⁻¹ y）在计算后验均值时已经得到。这使我们能获得一个学生t分布下的积分估计，其误差估计有所改善，但根本性的阶次不匹配问题依然存在，因为先验的平滑性假设没有改变。

要获得更快的收敛速度，我们需要使用更符合函数真实性质的先验。例如，对于光滑函数，可以使用基于正交多项式的高斯求积法则。

高斯求积法则（如高斯-勒让德、高斯-切比雪夫求积）是一类高效的积分方法，它们对应于在特定多项式函数空间上的最小二乘估计（即贝叶斯估计）。这些方法可以达到指数级或超多项式级的收敛速度，但前提是函数足够光滑。

这些经典求积法则也可以解释为贝叶斯求积，但需要针对特定的评估点数量专门设计核函数，使得在该数量点处后验方差恰好为零。这虽然在点估计上极其高效，但在不确定性量化方面不太自然。

高维挑战与自适应贝叶斯求积

一维积分问题已有成熟高效的解法。真正的挑战在于高维积分。基于网格或多项式展开的方法会遭受“维数灾难”，所需评估点数量随维度指数增长。

然而，机器学习中的积分通常具有特殊结构：它们常是某个函数在概率分布下的期望值。概率分布往往具有“集中性”，即大部分概率质量集中在空间的某个小区域内。

蒙特卡洛方法（如MCMC）通过在高概率区域花费更多时间来自适应这种结构，但其收敛速率 O(1/√N) 较慢，且该速率与维度无关。这看似是优点，实则意味着它对问题结构利用不足。

我们可以构建自适应的贝叶斯求积方法。例如，通过一个“扭曲”的高斯过程先验（如对高斯过程取平方以确保非负性），我们可以得到一个后验，它不仅在函数值上，也在其不确定性上依赖于观测数据。这允许算法计算一个依赖于数据的“效用函数”（如不确定性），并选择最大化该函数的点进行评估。这样的算法会主动探索可能具有高函数值的区域，并自适应地细化这些区域的估计。

这类自适应算法可以扩展到中等维度（如10-50维），并且被证明具有优于蒙特卡洛的收敛速率。在维度极高（如数千维）且问题缺乏结构时，蒙特卡洛仍然是可行的选择，但我们也将在后续课程中看到其他可能更优的确定性方法。

总结

本节课中，我们一起学习了如何将数值积分视为一个贝叶斯推断问题。

  1. 贝叶斯框架：我们通过为未知量（积分）和可计算量（函数值）设置联合先验，将计算问题概率化。高斯过程先验与积分的线性性质相结合，产生了贝叶斯求积算法。
  2. 经典算法的贝叶斯解释：我们发现，经典的梯形法则对应于使用维纳过程先验的贝叶斯后验均值估计。许多其他经典数值算法也都可以理解为某种贝叶斯估计。
  3. 不确定性量化：贝叶斯方法免费提供了点估计的不确定性（后验方差），尽管其校准程度取决于先验假设的合理性。
  4. 先验知识的重要性：算法的性能强烈依赖于先验假设。过于保守的先验（如假设函数粗糙）会导致性能下降和误差估计不准。利用问题的具体知识（如光滑性、非负性、集中性）可以设计出更高效、更自适应的算法。
  5. 作为智能体的计算：数值算法可以被视为在其计算环境中行动的自主智能体，它可以根据信息论准则（如最大化信息增益）来决定下一步计算什么。这种观点为设计自适应算法提供了原则性指导。

核心思想是：计算不应被视为机械的符号操作，而应被视为一个从有限信息中学习未知量的过程。 通过明确地使用先验知识并量化不确定性，我们可以构建出更强大、更可靠、有时甚至是革命性的数值方法。

011：深度学习优化 🧠

在本节课中，我们将要学习深度学习中的优化问题。我们将探讨为什么训练神经网络如此具有挑战性，当前有哪些训练方法，以及如何通过基准测试和调试工具来改进训练过程。

概述：为什么训练神经网络如此困难？ 🤔

上一节我们介绍了机器学习优化的基本概念，本节中我们来看看深度学习训练的特殊性。训练神经网络远非简单地选择一个优化器（如Adam）和一个默认学习率那么简单。现实世界中的例子，如Meta训练大型语言模型OPT的过程，表明这需要复杂的超参数调整、精心设计的学习率调度，以及大量的试错。目前，深度学习社区尚未就“最佳”训练方法达成共识，这使得训练过程充满了不确定性和挑战。

机器学习 vs. 优化：核心区别 🔄

在深入探讨训练方法之前，我们需要理解一个核心概念：机器学习不等同于优化。

在监督学习设置中，我们有一个由参数 $\theta$ 参数化的神经网络。给定输入 $x$，模型预测输出 $\hat{y}$。我们的目标是让预测 $\hat{y}$ 尽可能接近真实目标 $y$，这通过损失函数 $L$ 来衡量。

我们真正关心的是最小化模型在真实数据分布上的期望损失（真实风险）：
$$R(\theta) = \mathbb{E}{(x,y) \sim P{data}}[L(f(x; \theta), y)]$$

然而，我们无法访问无限的真实数据分布。我们只能访问一个有限的训练数据集 $D_{train} = {(x_i, y_i)}{i=1}^N$。因此，我们实际操作的是经验风险（经验损失）：
$$\hat{R}(\theta) = \frac{1}{N} \sum^{N} L(f(x_i; \theta), y_i)$$

更常见的是，我们使用小批量（mini-batch）进行训练，每次只计算一小部分数据（大小为 $B$）的梯度：
$$\hat{R}B(\theta) = \frac{1}{B} \sum^{B} L(f(x_i; \theta), y_i)$$

因此，我们优化的目标（经验风险）与我们最终关心的目标（真实风险）是不同的。这导致了过拟合、泛化等复杂问题。优化器（如SGD、Adam）旨在最小化经验风险，而我们的目标是获得一个在未见数据上表现良好的模型。因此，更准确的说法是，我们使用训练算法（Training Algorithms）来“训练”网络，而不仅仅是“优化”一个损失函数。

训练算法的现状：一片混乱的海洋 🌊

上一节我们明确了训练与优化的区别，本节中我们来看看当前可用的训练算法。由于上述的复杂性，研究人员提出了大量不同的训练算法。

以下是几种经典算法的更新公式示例：

  • 随机梯度下降（SGD）:
    $$\theta_{t+1} = \theta_t - \eta \nabla \hat{R}_B(\theta_t)$$
    其中 $\eta$ 是学习率。

  • 带动量的SGD（Momentum）:
    $$v_{t+1} = \rho v_t + \nabla \hat{R}B(\theta_t)$$
    $$\theta     = \theta_t - \eta v_{t+1}$$
    其中 $\rho$ 是动量参数。

  • Adam:
    $$m_t = \beta_1 m_{t-1} + (1-\beta_1) \nabla \hat{R}B(\theta_t)$$
    $$v_t = \beta_2 v     + (1-\beta_2) (\nabla \hat{R}_B(\theta_t))^2$$
    $$\hat{m}_t = m_t / (1-\beta_1^t), \quad \hat{v}t = v_t / (1-\beta_2^t)$$
    $$\theta     = \theta_t - \eta \cdot \hat{m}_t / (\sqrt{\hat{v}_t} + \epsilon)$$
    其中 $\beta_1, \beta_2$ 是衰减率，$\epsilon$ 是稳定项。

目前存在超过150种不同的训练算法变体。问题不仅在于选择哪个算法，还在于如何为每个算法设置其超参数（如学习率、动量参数等）。例如，Adam就有学习率、$\beta_1$、$\beta_2$ 和 $\epsilon$ 四个主要超参数需要调整。这种复杂性使得确定一个通用的、最优的训练“配方”变得极其困难。

基准测试的挑战：如何比较训练算法？ 📊

面对如此多的算法，一个自然的想法是通过基准测试来找出最优者。然而，对深度学习训练算法进行公平、全面的基准测试本身就是一个巨大的挑战。

以下是基准测试面临的主要问题：

  1. 定义“更好”： “更好”可能意味着更快的训练速度、更高的最终精度、对超参数更鲁棒、或需要更少的调参。
  2. 随机性： 深度学习训练具有固有的随机性（如权重初始化、数据打乱）。单次运行的结果可能是噪音。需要多次重复运行以获得统计显著性，这极大地增加了计算成本。
  3. 泛化性： 一个算法在某个特定任务（如CIFAR-10图像分类）上表现好，不代表它在其他任务（如语言建模、强化学习）上也表现好。需要在多样化的任务上进行测试。
  4. 超参数交互： 算法的性能严重依赖于其超参数设置。比较时必须考虑公平的调优预算。此外，最佳学习率调度（如常数、余弦衰减、热身重启）也会与算法产生交互。
  5. 算法家族 vs. 具体实例： 像“SGD”或“Adam”是一个算法家族，包含无数超参数设置的具体实例。我们只能比较具体的、完全定义的算法实例。

一项大规模基准测试研究尝试了15种流行算法，在8个不同任务上，使用4种调优预算和4种学习率调度进行组合测试，总共进行了超过5万次独立训练运行。结果显示，并没有一个算法在所有任务和设置下 consistently 表现最佳。像Adam这样的流行算法表现稳健，但某些任务上其他算法（如RMSprop）可能更好。这证实了当前缺乏一个明确的、通用的“最先进”训练协议。

利用额外信息：从梯度分布到调试工具 🛠️

既然没有“银弹”算法，我们如何改善训练体验呢？一个策略是挖掘训练过程中未被充分利用的额外信息，并用它来构建更好的工具或算法。

目前，训练算法只使用了小批量梯度的均值：
$$g_{batch} = \frac{1}{B} \sum_{i=1}^{B} \nabla L(f(x_i; \theta), y_i)$$

然而，我们实际上计算了每个样本的独立梯度 $\nabla L(f(x_i; \theta), y_i)$。这些独立梯度包含了关于梯度估计噪声或不确定性的宝贵信息（例如，它们的方差）。这些信息目前被框架（如PyTorch）在求平均后丢弃了。

通过工具（如Backpack库），我们可以高效地访问这些独立梯度及其统计信息（如方差）。这为理解和控制训练过程打开了新的大门。

一个关键的应用是构建深度学习调试工具。传统的调试器对于检查代码逻辑很有用，但对于诊断训练效率低下（如错误的学习率、数据未归一化、有问题的模型架构）帮助有限。我们需要能压缩高维训练状态（如数百万权重的轨迹）成为人类可理解的、有意义的“状态报告”的工具。

工具Cockpit就是这样一个概念验证。它像飞机驾驶舱一样，为训练过程提供了一系列“仪表盘”：

  • 步长信息： 例如，$\alpha$ 指标。它衡量当前步长相对于该方向局部最优步长是过小（负值，欠步）还是过大（正值，过冲）。$\alpha$ 是一个标量，高效可计算，能直接提示是否应调整学习率。
  • 梯度信息： 显示梯度范数、梯度元素分布直方图、以及基于样本梯度的噪声估计（如梯度方差），这有助于判断批量大小是否合适。
  • 曲率信息： 提供损失函数局部曲率的见解（将在后续讲座详细讨论）。

这些工具可以帮助发现不同类型的“训练Bug”：

  • 数据Bug： 例如，输入数据未正确归一化，这在像素级视觉检查中难以发现，但梯度直方图会显示出异常分布。
  • 模型Bug： 例如，使用某些激活函数导致梯度在深层网络消失，通过逐层的梯度直方图可以清晰看到。
  • 超参数Bug： 例如，学习率设置不当。$\alpha$ 指标显示持续负值（欠步）是学习率过小的明确信号。

有趣的是，基准测试和Cockpit的分析都表明，在每一步都追求局部最优（$\alpha \approx 0$）通常并不是最佳策略。由于梯度噪声的存在，适度的“过冲”（$\alpha > 0$）有时能带来更好的长期性能和泛化能力。这再次强调了“训练”与“优化”的根本区别：在随机环境下，短视的最优步长可能导致长期性能下降。

总结与展望 🎯

本节课中我们一起学习了深度学习优化的核心挑战与现状。

我们首先明确了机器学习训练不等同于损失函数优化，因为我们操作的是有噪声的经验风险，而目标是未知数据上的泛化性能。

接着，我们看到了当前训练算法的现状是一片“混乱的海洋”，存在海量算法且缺乏明确的最佳实践，训练过程仍需大量人工调参和“保姆式”监控。

然后，我们探讨了进行公平基准测试的巨大挑战，以及现有大规模测试表明尚无一个算法能通用地统治所有任务。

最后，我们提出了一条改进路径：利用训练过程中未被充分利用的额外信息（如梯度分布）。这既能用于构建像Cockpit这样的调试工具，帮助实践者更高效地诊断问题，也为未来开发更自动化、更智能的训练算法奠定了基础，以期最终将人们从繁琐的超参数调整中解放出来。

总之，当前神经网络的训练仍然是一个复杂且尚未完全解决的挑战，但通过更严谨的评估和更深入的系统性分析，我们正在朝着建立更可靠、更高效的训练方法学迈进。

012：深度学习的二阶优化

在本节课中，我们将要学习深度学习中二阶优化的基本概念、优势以及面临的挑战。我们将从一阶优化方法的问题出发，介绍牛顿法等二阶方法的核心思想，并探讨其在非凸、高成本和随机性环境下的应用策略。

概述

上周的课程和练习中，我们了解到使用一阶方法（如SGD和Adam）训练神经网络是一个耗时、昂贵且不稳定的过程。本节课，我们将学习另一类优化器——二阶优化方法。这类方法至少可以避免一阶方法的一些问题，但同时也带来了新的挑战和未解决的难题，特别是在深度学习领域。

回顾：深度学习优化问题

在深入探讨二阶方法之前，我们先快速回顾一下深度学习优化的基本框架。

我们拥有训练数据，即输入 x 和输出 y 的样本，它们来自一个数据生成分布 P_data。我们假设这些样本是独立同分布的。

我们有一个神经网络模型 f(θ, x)，它基于参数向量 θ 和输入 x 产生一个 C 维的输出向量（例如，C 是类别数）。

我们通过损失函数 ℓ(ŷ, y) 来比较预测值 ŷ 和真实标签 y 之间的差异。我们的目标是找到最小化期望风险的参数 θ：

R(θ) = E_{(x,y) ~ P_data}[ℓ(f(θ, x), y)]

由于我们无法直接访问真实分布 P_data，我们转而最小化其蒙特卡洛估计量，即经验风险：

R_emp(θ) = (1/N) Σ_{i=1}^{N} ℓ(f(θ, x_i), y_i)

梯度下降法通过沿负梯度方向更新参数来最小化目标函数：

θ_{k+1} = θ_k - α * ∇R(θ_k)

其中 α 是学习率。在深度学习中，我们通常使用其随机版本——随机梯度下降（SGD），它使用小批量数据计算梯度的无偏估计：

θ_{k+1} = θ_k - α * ĝ(θ_k), 其中 ĝ(θ_k) = (1/|B|) Σ_{i in B} ∇ℓ(f(θ, x_i), y_i)

一阶方法的局限性：条件数问题

梯度类方法的效率受到目标函数条件数的严重影响。条件数 κ 定义为最大曲率 L 与最小曲率 μ 之比：

κ = L / μ

当 κ 接近1时（各方向曲率相近），梯度方向指向最小值，收敛迅速。当 κ 很大时（存在平坦和陡峭方向），梯度方向不佳，收敛速度会变得非常慢。梯度下降法的收敛定理表明，其收敛速度与 (κ - 1)/(κ + 1) 相关，当 κ 很大时，这个因子接近1，意味着每一步的进展微乎其微。

为了消除对条件数的依赖，我们需要引入逆曲率信息来重新缩放梯度。这直接引向了二阶优化方法。

二阶方法的核心：牛顿步

二阶方法的核心思想是利用曲率或其逆来改进更新步长。其核心是牛顿步。

牛顿步基于在当前迭代点 θ_k 处的损失函数二次泰勒近似：

Q(θ) = R(θ_k) + ∇R(θ_k)^T (θ - θ_k) + 1/2 (θ - θ_k)^T H(θ_k) (θ - θ_k)

其中 H(θ_k) 是海森矩阵（二阶导数矩阵）。我们假设 H 是正定的。

牛顿法的思想是直接跳到这个二次近似函数的最小值点。通过求解 ∇Q(θ) = 0，我们得到牛顿更新公式：

θ_{k+1} = θ_k - H(θ_k)^{-1} ∇R(θ_k)

与线性收敛的一阶方法不同，在适当条件下（如初始点足够接近最优解，且海森矩阵 Lipschitz 连续），牛顿法可以实现局部二次收敛，速度远快于线性收敛。

牛顿法的优势与挑战

上一节我们介绍了牛顿步的推导，本节我们来看看牛顿法相比一阶方法的优势，以及将其应用于深度学习时面临的核心挑战。

优势总结：

  1. 对病态问题鲁棒：通过逆海森矩阵重新缩放梯度，消除了对条件数的依赖。
  2. 收敛速度快：在最优解附近可实现二次收敛。
  3. 自然的步长：理论上不需要像梯度下降那样手动调整学习率（1/L）。
  4. 广泛应用：是许多科学计算领域的默认选择。

核心挑战：
尽管牛顿法看起来很完美，但在深度学习中却很少使用，主要因为三个问题：

  1. 非凸性：深度学习损失函数通常是非凸的，海森矩阵可能不定（有正有负特征值），导致牛顿步无定义（二次模型无最小值）。
  2. 成本：海森矩阵的维度是参数数量 d 的平方（O(d^2)），存储和求逆（O(d^3)）对于百万级参数的神经网络来说计算上不可行。
  3. 随机性：我们只能基于有限数据（小批量）估计梯度和海森矩阵，这种随机性会带来偏差和不稳定性。

接下来，我们将逐一探讨这些挑战的应对策略。

应对挑战一：非凸性

上一节我们提到非凸性导致海森矩阵不定。本节中我们来看看如何通过寻找正定的曲率代理和使用阻尼技术来解决这个问题。

使用正定曲率代理：GGN 和 Fisher 矩阵

我们可以用其他保证正定（或半正定）的矩阵来替代海森矩阵。

1. 广义高斯-牛顿矩阵
GGN 矩阵基于将损失函数分解为模型输出 f(θ, x) 和损失函数 ℓ 的组合。通过链式法则推导海森矩阵，并忽略模型中二阶导数项（即假设模型在参数上是线性的），我们得到 GGN 矩阵 G：

G(θ) = E[J(θ)^T H_ℓ(f(θ, x), y) J(θ)]

其中 J(θ) 是模型 f 关于参数 θ 的雅可比矩阵，H_ℓ 是损失函数 ℓ 关于其第一个参数（模型输出）的海森矩阵。对于常见的损失函数（如均方误差、交叉熵），H_ℓ 是正定（或半正定）的，因此 G 也是正半定的。

2. Fisher 信息矩阵
当损失函数可以写成负对数似然形式时（例如使用 softmax 和交叉熵损失），我们可以定义 Fisher 信息矩阵 F：

F(θ) = E[∇ log p(y | x; θ) ∇ log p(y | x; θ)^T]

其中 p(y | x; θ) 是由模型参数定义的条件概率分布。Fisher 矩阵也是正半定的。使用 Fisher 逆乘以梯度的更新，可以解释为在分布空间（使用 KL 散度作为度量）的最速下降。

在许多情况下（如使用交叉熵损失的分类问题），GGN 和 Fisher 矩阵是相等的。它们都提供了有意义的曲率代理。

使用阻尼控制更新保守性

即使使用正定曲率矩阵，二次模型也可能在远离当前点时变得不准确。信赖域方法通过限制更新步长在一个可信半径 R 内来解决这个问题。这等价于在曲率矩阵中添加一个阻尼项 δI：

θ_{k+1} = θ_k - (G(θ_k) + δI)^{-1} ∇R(θ_k)

通过调整阻尼参数 δ，我们可以在牛顿步（δ=0）和梯度下降步（δ 很大）之间平滑插值。常用的启发式方法（如 Levenberg-Marquardt）会根据实际损失下降与模型预测下降的比率来动态调整 δ。

应对挑战二：高计算成本

上一节我们解决了非凸性问题，本节我们来看看如何应对海森矩阵（或其代理）存储和求逆的巨大计算成本。我们将介绍三种主流的近似方法。

1. 拟牛顿法：BFGS 和 L-BFGS

核心思想是逐步从梯度观测中学习逆海森矩阵的近似，而不显式构造或存储海森矩阵。

  • 基础：BFGS：它基于“割线方程” H_k s_k = y_k，其中 s_k = θ_{k+1} - θ_k，y_k = ∇R(θ_{k+1}) - ∇R(θ_k)。BFGS 通过迭代更新一个满足对称性、正定性和最小变化条件的逆海森近似矩阵 B_k。更新公式仅涉及向量 s_k 和 y_k。
  • 内存优化：L-BFGS：标准的 BFGS 需要存储稠密的 d×d 矩阵 B_k。L-BFGS 的改进在于，它发现 B_k 可以由初始矩阵 B_0（通常设为 βI）和最近的 m 对向量 {s_i, y_i} 递归构造。通过只保存最近的 m 对向量（m 通常很小，如 20），L-BFGS 大幅降低了内存消耗，并且有高效算法计算 B_k * g。

2. Hessian-Free 优化

核心思想是使用共轭梯度法来近似求解牛顿步 H^{-1}g，因为 CG 法只需要计算矩阵-向量乘积 H * v，而无需显式构造或存储 H。

  • 实践：在深度学习中，通常使用 GGN 矩阵 G 作为曲率矩阵。计算 G * v 可以通过高效的自动微分技巧实现，无需构造庞大的雅可比矩阵 J。技巧在于利用方向导数的定义：J * v = ∂/∂ε f(θ + εv) |_{ε=0}，这可以通过一次前向传播的高阶模式（如 torch.autograd.grad）高效计算。
  • 优势与局限：内存占用小，更接近“纯”牛顿步。但需要多次 CG 迭代，计算可能仍很耗时，并且通常需要较大的批量大小，这可能与批归一化层不兼容。

3. K-FAC：克罗内克因子分解近似曲率

核心思想是针对 Fisher 信息矩阵 F，利用神经网络的结构特性进行两种近似，使其易于求逆。

  1. 块对角近似：假设 Fisher 矩阵是块对角的，每个层对应一个块。这忽略了层间的曲率相互作用。
  2. 期望与克罗内克积交换近似：对于线性层，每个块可以精确写为两个小矩阵的克罗内克积的期望。K-FAC 将期望与克罗内克积运算交换，近似为两个期望的克罗内克积：E[A ⊗ B] ≈ E[A] ⊗ E[B]。

经过这两种近似，每个块 F_i 近似为 Ā_i ⊗ B̄_i。其逆矩阵可以轻松计算：(Ā_i ⊗ B̄_i)^{-1} = Ā_i^{-1} ⊗ B̄_i^{-1}。由于 Ā_i 和 B̄_i 的维度很小（与层的输入输出维度相关），存储和求逆成本极低。

方法对比总结：
以下是三种方法的简要对比：

  • BFGS/L-BFGS：动态低秩近似。是小规模确定性问题的默认选择，在参数少、噪声小的问题上可能很快。
  • Hessian-Free：更接近经典牛顿步，内存需求小，但计算可能更顺序化，且受批量大小限制。
  • K-FAC：轻量级 Fisher 矩阵表示，广泛使用。但有研究表明其性能可能很大程度上依赖于启发式阻尼，而非精确的曲率信息。

应对挑战三：随机性

上一节我们讨论了降低计算成本的方法，本节我们来看看最棘手的挑战：随机性。当我们只能使用小批量数据估计梯度和曲率时，即使估计是无偏的，牛顿步本身也会产生偏差。

随机牛顿步的偏差问题

假设我们有梯度的无偏估计 ĝ 和海森矩阵的无偏估计 Ĥ。我们希望的更新是 H^{-1}g。但我们实际计算的是 Ĥ^{-1}ĝ。即使 ĝ 和 Ĥ 各自无偏，Ĥ^{-1}ĝ 的期望也不等于 H^{-1}g。核心原因在于：求逆运算不是线性运算。

一维直观解释：
假设梯度估计完美（g=1），问题简化为比较 E[1/ĥ] 和 1/E[ĥ] = 1/h。由于 ĥ 是一个随机变量（分布如图），而函数 φ(x) = 1/x 在 x>0 时是凸函数。根据詹森不等式，E[φ(ĥ)] ≥ φ(E[ĥ])，即 E[1/ĥ] ≥ 1/h。这意味着逆曲率的估计在期望上是偏大的。特别是当 ĥ 的分布接近零时，求逆会严重扭曲分布，产生重尾，将期望值拉高。

多维影响：
这种偏差对不同方向的影响不同。在高曲率方向（h 值大），分布远离零，偏差较小。在低曲率方向（h 值小），分布可能接近零，求逆带来的偏差会被放大，导致该方向的更新步长在期望上过长。最终，整个牛顿步的方向和大小都可能被严重扭曲。

不稳定性问题

除了系统性偏差，随机性还会导致不稳定性。纯粹由于随机采样，我们可能偶然得到一个非常接近奇异的曲率估计 Ĥ。对其求逆会导致更新步长异常巨大，从而使优化过程崩溃。

阻尼的作用与局限

添加阻尼 δI 是缓解这些问题的一种方法。它将曲率估计的分布向右移动（远离零），从而减轻求逆带来的偏差放大效应和不稳定性。

然而，阻尼只是一个标量，平等地对待所有方向，这可能过于简化。此外，阻尼值通常通过启发式方法调整。理想情况下，我们需要更智能的算法，能够在求逆过程中就意识到估计的不确定性，并据此决定在不同方向上采取何种更新。

总结与展望

本节课我们一起学习了深度学习中二阶优化的核心内容。

我们首先回顾了一阶优化方法（如SGD）在训练深度神经网络时存在的耗时、昂贵和对病态条件敏感的问题。二阶方法（以牛顿法为代表）通过引入逆曲率信息，提供了对病态问题鲁棒、收敛速度更快的潜力。

然而，将二阶方法直接应用于深度学习面临三大挑战：

  1. 非凸性：我们介绍了使用正定曲率代理（如GGN和Fisher矩阵）和阻尼技术来获得定义良好且保守的更新步长。
  2. 高计算成本：我们探讨了三种主流近似方案：BFGS/L-BFGS（动态低秩近似）、Hessian-Free优化（使用CG和矩阵-向量积）以及K-FAC（基于网络结构的块对角克罗内克因子近似）。它们各自在内存、计算精度和适用性上有不同的权衡。
  3. 随机性：我们深入分析了基于随机估计的牛顿步存在的固有偏差和不稳定性问题。阻尼是当前主要的缓解手段，但这问题在根本上尚未完全解决，可能需要能够利用估计值分布信息而不仅是点估计的新算法。

尽管存在挑战，曲率信息仍然是理解优化问题、设计新算法以及进行不确定性量化（下节课主题）的宝贵工具。二阶优化方法为改进深度学习训练效率提供了重要的思路和部分解决方案。

013：深度学习中的不确定性

概述

在本节课中，我们将要学习如何为深度学习模型引入不确定性。我们将探讨为什么在深度神经网络中考虑不确定性至关重要，并介绍一种高效、低成本的方法来实现这一点，即拉普拉斯近似。我们将看到，这种方法不仅能解决深度分类器中的过度自信问题，还能为模型选择和参数调整提供新的工具。

为什么深度学习需要不确定性？

上一节我们介绍了在计算中引入不确定性的概念。本节中，我们来看看为什么在深度学习中，尤其是在深度神经网络中，考虑不确定性是必要的。

深度神经网络存在一些令人惊讶的具体问题，其中一个关键问题是过度自信。

以下是一个典型的例子：当向一个在ImageNet等数据集上训练良好的分类网络输入一些看起来像白噪声、与训练数据分布截然不同的图像时，网络仍会以极高的置信度（例如99.5%）将其归类为某个特定类别（如“知更鸟”或“猎豹”）。这显然不是我们期望的行为。

这种现象被称为渐近过度自信。2019年，Matthias Hein等人的一篇论文从理论上证明了，对于使用ReLU激活函数的标准分类网络，这是一个根本性问题。

其论证基于两个关键点：

  1. ReLU网络是分段线性的：由ReLU激活函数构成的深度神经网络，其最终输出是一个分段线性函数。
  2. 远离数据区域的行为：在输入空间中，当点远离训练数据时，最终会进入一个线性区域。在这个区域内，网络的输出（即softmax层的输入）是线性的。由于不同线性函数的“增益”不同，随着距离无限增大，增益最高的那个输出会无限地超过其他输出。当这个结果输入到softmax函数时，就会导致对某一个类别的预测置信度无限趋近于1。

公式描述：对于一个ReLU分类器，其输出函数 ( f_{\theta}(x) ) 是分段线性的。在远离数据的某个线性区域内，softmax的输入是 ( W^T \phi(x) )，其中 ( W ) 是最后一层权重，( \phi(x) ) 是倒数第二层的输出。当 ( |x| \to \infty ) 时，最大分量的值会主导softmax，导致预测概率趋近于1。

核心结论：要根本性地解决这个问题，就需要对神经网络的权重引入不确定性。

贝叶斯视角下的深度学习

在探讨如何引入不确定性之前，我们首先需要建立一个贝叶斯解释框架。

标准的深度神经网络训练是通过最小化经验风险（或损失函数）来进行的：

公式：
[
\mathcal{L}(\theta) = \sum_{i=1}^{N} \ell(y_i, f_{\theta}(x_i)) + \lambda R(\theta)
]

其中，( \theta ) 是网络权重，( \ell ) 是损失函数（如交叉熵），( R(\theta) ) 是正则化项（如权重衰减）。

这个最小化过程可以重新解释为贝叶斯推断：

  • 负对数似然：损失项 ( \sum \ell(y_i, f_{\theta}(x_i)) ) 可以看作数据在模型下的负对数似然。
  • 对数先验：正则化项 ( \lambda R(\theta) ) 可以看作权重的负对数先验。

因此，最小化总损失 ( \mathcal{L}(\theta) ) 等价于最大化后验概率 ( p(\theta | \mathcal{D}) )（最大后验估计，MAP）。公式：
[
p(\theta | \mathcal{D}) \propto \exp(-\mathcal{L}(\theta))
]

所以，深度学习已经在进行一种贝叶斯推断，只不过它只计算后验的众数（mode），而忽略了后验的完整形状。

为了获得完整的不确定性，我们需要考虑整个后验分布 ( p(\theta | \mathcal{D}) )，并在做预测时对权重进行积分（边缘化）：
公式：
[
p(y^* | x^, \mathcal{D}) = \int p(y^ | x^*, \theta) p(\theta | \mathcal{D}) d\theta
]

然而，对于深度神经网络，这个后验分布是高维且非线性的，上述积分是难以处理的。传统的马尔可夫链蒙特卡洛方法计算成本过高，不适用于大规模深度学习。

拉普拉斯近似：一种低成本方案

我们需要一种既高效又能引入不确定性的方法。本节介绍的方法结合了自动微分和线性代数，是成本最低的积分方法之一。

一个关键的发现是：即使只是对网络的权重引入一点点不确定性（例如一个高斯分布），也能修复过度自信的问题。具体来说，任何对网络最后一层权重的高斯近似度量，都能部分解决这个问题。

这引出了我们的核心工具：拉普拉斯近似。

拉普拉斯近似的思想：在找到损失函数的最小值点（MAP估计 ( \theta_{MAP} ) ）后，在该点对负对数后验（即损失函数 ( \mathcal{L}(\theta) ) ）进行二阶泰勒展开。

公式推导：

  1. 在 ( \theta_{MAP} ) 处展开：
    [
    \mathcal{L}(\theta) \approx \mathcal{L}(\theta_{MAP}) + \frac{1}{2} (\theta - \theta_{MAP})^T \mathbf{H} (\theta - \theta_{MAP})
    ]
    其中 ( \mathbf{H} = \nabla^2 \mathcal{L}(\theta_{MAP}) ) 是损失函数在 ( \theta_{MAP} ) 处的海森矩阵。由于在极值点梯度为零，一阶项消失。
  2. 代入后验公式：
    [
    p(\theta | \mathcal{D}) \propto \exp(-\mathcal{L}(\theta)) \approx \exp\left(-\mathcal{L}(\theta_{MAP}) - \frac{1}{2} (\theta - \theta_{MAP})^T \mathbf{H} (\theta - \theta_{MAP})\right)
    ]
  3. 忽略常数项，我们得到一个高斯近似：
    [
    p(\theta | \mathcal{D}) \approx \mathcal{N}(\theta_{MAP}, \mathbf{H}^{-1})
    ]

这样，我们就用了一个以 ( \theta_{MAP} ) 为均值、以海森矩阵逆为协方差的高斯分布来近似真实的后验。

优势：

  1. 低成本：主要成本是训练网络（本来就需要做）。训练完成后，只需额外计算一次海森矩阵。
  2. 事后处理：可以对任何已训练好的网络应用此方法，无需重新训练。
  3. 保持点估计：均值就是原网络的权重，因此不会改变模型原有的点预测性能。

实践中的挑战与解决方案：广义高斯-牛顿矩阵

上一节我们介绍了拉普拉斯近似的理想形式。然而，在实践中直接应用会面临几个挑战：

  1. 优化器可能并未收敛到真正的局部极小点。
  2. 损失函数的曲率（海森矩阵）可能不是正定的。
  3. 海森矩阵的尺寸是权重数量的平方，计算和存储成本极高。

解决方案与上一讲中处理优化问题的方法一致：使用广义高斯-牛顿矩阵作为海森矩阵的近似。

广义高斯-牛顿矩阵：
[
\mathbf{G} = \mathbf{J}^T \mathbf{H}f \mathbf{J}
]
其中，( \mathbf{J} ) 是网络输出 ( f(x) ) 相对于权重 ( \theta ) 的雅可比矩阵，( \mathbf{H}_f ) 是损失函数相对于网络输出的海森矩阵（对于分类问题常为常数）。

GGN矩阵总是半正定的，并且计算起来比完整的海森矩阵更高效。通常我们还会进一步对其做近似，例如只考虑对角线元素（对角近似）或分块对角结构。

使用GGN矩阵与网络线性化的概念紧密相关。线性化是指在权重空间中对网络函数 ( f_{\theta}(x) ) 进行一阶泰勒展开：
公式：
[
f_{\theta}(x) \approx f_{\theta_{MAP}}(x) + \mathbf{J}(x) (\theta - \theta_{MAP})
]
这表示网络输出在权重空间中是线性的，但在输入空间 ( x ) 中仍然是非线性的（因为 ( \mathbf{J}(x) ) 依赖于 ( x )）。这保留了深度网络的表达能力。

结合线性化和拉普拉斯近似（使用GGN），我们得到了一个完全可处理的模型。

线性化拉普拉斯近似带来的预测

当我们结合线性化和拉普拉斯高斯近似后，对于不同的任务，预测变得可处理。

1. 回归任务（二次损失）：
预测分布是高斯分布，可以直接解析计算。
公式：
[
p(y^* | x^, \mathcal{D}) = \mathcal{N}(f_{\theta_{MAP}}(x^), \mathbf{J}(x*)T \mathbf{G}^{-1} \mathbf{J}(x^*) + \sigma^2)
]
其中 ( \sigma^2 ) 是观测噪声。

2. 分类任务（交叉熵损失/Softmax输出）：
需要计算Softmax在高斯分布下的期望，这个积分没有闭式解。但可以使用一个简单而有效的近似（由David MacKay提出）：
公式：
[
p(y^=c | x^, \mathcal{D}) \approx \text{softmax}\left( \frac{\mu_c}{\sqrt{1 + \frac{\pi}{8} \sigma_c^2}} \right)
]
其中 ( \mu_c = [f_{\theta_{MAP}}(x^)]_c )，( \sigma_c^2 = [\mathbf{J}(x*)T \mathbf{G}^{-1} \mathbf{J}(x^)]_{c,c} )。

这个近似公式直观地显示了不确定性（( \sigma_c^2 )）如何“压平”原始的Softmax输出，从而降低置信度。当远离数据时，( \mu_c ) 线性增长，而 ( \sigma_c^2 ) 二次增长，其平方根线性增长，因此比值趋于常数，避免了置信度趋近于1。

代码描述：在实践中，流程如下：

# 1. 正常训练一个神经网络，得到 theta_MAP
model.train()
# ... 训练循环 ...

# 2. 训练完成后，计算GGN矩阵（例如使用PyTorch的自动微分）
loss = criterion(model(X_train), y_train)
hessian = compute_hessian(loss, model.parameters()) # 实际中常用GGN近似
jacobian = compute_jacobian(model(X_train), model.parameters())

# 3. 对于新输入 x_star，进行预测
with torch.no_grad():
    mean_pred = model(x_star)
    jac_star = compute_jacobian_single(model, x_star)
    variance = jac_star.T @ torch.inverse(hessian) @ jac_star
    # 对于分类，应用上述近似公式
    adjusted_logits = mean_pred / torch.sqrt(1 + (math.pi/8) * variance)
    prob_pred = torch.softmax(adjusted_logits, dim=-1)

拉普拉斯近似的应用：超越不确定性量化

线性化拉普拉斯近似不仅提供了预测不确定性，还解锁了其他有用的功能。

1. 模型选择与证据最大化
在贝叶斯框架中，我们可以计算模型证据（边缘似然）( p(\mathcal{D} | \mathcal{M}) ) 来比较不同模型。在线性化拉普拉斯近似下，这变得可计算：
公式：
[
\log p(\mathcal{D} | \mathcal{M}) \approx -\mathcal{L}(\theta_{MAP}) + \frac{1}{2} \log \det(\mathbf{G}) - \frac{d}{2} \log(2\pi)
]
其中 ( d ) 是权重的维度。

我们可以利用这个近似证据来选择神经网络的超参数，例如：

  • 网络层数或宽度
  • 权重先验的精度（正则化强度）

通过最大化证据，我们可以以概率化的方式自动选择这些参数，而无需单独的验证集。

2. 实现渐近校准置信度
基本的线性化拉普拉斯近似可以防止置信度趋近于1，但远离数据时，置信度会趋于一个非均匀的常数，而不是最大熵状态（各类别等概率）。为了达到更好的校准，我们可以添加一个非参数化的“后备”高斯过程先验，该先验的方差随着到训练数据中心的距离立方而增长。这样，在渐近区域，预测的不确定性将趋于均匀分布。

3. 在输入域内训练不确定性
我们可以通过向网络中添加特殊的“不确定性单元”来主动塑造不确定性。这些单元的权重均值设为零（不影响点预测），但其方差（通过拉普拉斯近似的曲率体现）是可训练的。我们可以定义额外的损失函数，迫使网络在指定的“分布外”区域（例如，与训练类别不同的图像集）具有高不确定性。这提供了一种将先验知识关于“何处应不确定”注入模型的方法。

总结

本节课中我们一起学习了如何为深度学习引入不确定性，并重点介绍了一种高效实用的方法——线性化拉普拉斯近似。

核心要点：

  1. 不确定性至关重要：忽略权重不确定性会导致如ReLU分类器渐近过度自信等病理问题。
  2. 贝叶斯解释：标准深度学习可视为寻找最大后验估计。
  3. 拉普拉斯近似：通过在MAP点进行二阶近似，将复杂的后验转化为高斯分布 ( \mathcal{N}(\theta_{MAP}, \mathbf{G}^{-1}) )，其中 ( \mathbf{G} ) 是广义高斯-牛顿矩阵。
  4. 网络线性化：在权重空间对网络进行一阶展开，结合拉普拉斯近似，可将深度网络转化为一个高斯过程，从而实现预测不确定性的解析或近似计算。
  5. 广泛应用：此框架不仅提供了预测不确定性，还可用于模型选择（证据最大化）、实现更好的置信度校准，以及根据需求在特定输入区域塑造不确定性。

最终启示：严格的贝叶斯推断（追求精确后验）可能成本高昂，但概率化的思维（引入一个合理的概率度量）本身可以解决许多实际问题（如过度自信），且成本可控。线性化拉普拉斯近似正是这种思想的一个成功体现，它让我们能够在保持深度学习强大性能的同时，获得宝贵的不确定性信息。

014：结论 🎓

在本节课中，我们将回顾整个课程的核心内容，从线性代数到优化问题，探讨数值计算在机器学习中的核心地位。我们将看到，机器学习中的计算本质上是信息的管理与推断。

课程回顾与核心思想

这是本课程的最后一讲。我将快速回顾整个课程的内容，从十月中旬至今，涵盖了几乎所有讲座的核心思想。本次回顾旨在帮助大家为周一的考试做准备，同时也让我们退一步思考，在这个学期里我们究竟学到了什么。

机器学习本身就是数值计算。机器在学习过程中所执行的操作，本质上是在解决一个数值问题。这与经典AI（如算法与数据结构中的算法）不同。在机器学习中，我们需要解决那些没有闭式解的数学问题，这些问题正是经典数值分析所研究的。

这些问题主要包括（并非完全按课程顺序）：

  • 积分：用于贝叶斯推断。当需要衡量假设空间中剩余的不确定性体积时，就需要解决积分问题。
  • 优化：用于训练深度神经网络或任何通过最小化正则化经验风险来构建点估计的模型。
  • 微分方程：出现在智能体与数据源交互（如强化学习、机器人学）或科学定律建模等场景中，用于预测系统演化。
  • 线性代数：因为线性代数无处不在，它是解决上述所有问题（如高斯积分、二次优化、线性微分方程）的基础。

因此，我们逐一探讨了这些问题。我们首先从线性代数开始。

线性代数：计算的核心 🧮

上一节我们介绍了机器学习中数值问题的多样性。本节中我们来看看作为这一切基础的线性代数。

对我而言，线性代数部分的高级要点可以概括为：进行线性代数计算（这个动词）发生在你的处理单元（CPU或GPU）上，而你的数据存储在磁盘上。因此，当你考虑涉及数据集的线性代数计算的复杂度时，你实际上需要考虑的是从磁盘读取了多少数据并在芯片上进行了多少处理，而不仅仅是磁盘上有多少数据。

这一点之所以重要，是因为完整读取一次数据的成本与数据大小成线性关系。而线性代数中昂贵的操作与读取数据无关，而与你在芯片上对读取的数字所做的处理有关。这意味着我们实际上并不总是需要加载整个数据集。

机器学习中的一个基础案例是高斯过程回归，它在本课程中扮演了双重角色：一方面，它是最简单的监督机器学习算法之一；另一方面，它也是我们为课程后续部分构建数值算法的模型类。

以下是第2讲（由Marvin讲授，Unnaan准备）中的幻灯片。为了得到这个预测分布（包含均值μ和协方差Σ），我们需要求解一个线性方程组，通常写作矩阵求逆乘以向量的形式：(K(X*, X) * (K(X, X) + σ²I)^(-1) * y)。

这种表示法隐藏了许多有趣的结构。它只是告诉我们“需要什么”，而没有说明“如何计算”。在像深度学习这样的其他机器学习算法中，这一点尤为重要。因为可以写出这样的方程，人们就认为这就是高斯过程回归所需的全部，从而认为它计算昂贵（因为教科书上说求解此类问题是O(n³)的）。而对于深度学习，我们无法写出这样的方程，人们甚至不会开始思考训练深度神经网络的复杂度到底有多高。

我们注意到，同样的问题（矩阵求逆）会出现两次，只是每次乘以的向量不同。在第2讲和随后的第3讲中，Marvin和Jonathan指出，求解此类线性方程组有多种算法。

经典的算法是Cholesky分解，通常被当作一个“黑箱”呈现。但实际上，我们可以将这个算法视为一个迭代过程：它逐行/逐列加载大矩阵，每一步的操作成本最初是线性的，但随着已处理步骤的增加，会有一个二次的修正步骤。最终，因为运行了n步，总成本是立方的O(n³)。

这个幻灯片揭示的一个关键点是：你可以将Cholesky分解不仅视为最终构造一个矩阵，实际上在过程的每一步都可以停止，并返回一个对所需逆矩阵的估计（即矩阵C）。这引发了一些有趣的讨论。

为什么这很有趣？首先，你可能会想，我可以提前停止算法并使用这个近似。当然，你需要担心这个近似是否足够好。在那次讲座中我们看到，如果数据是按某种顺序处理的，并且数据本身是有序的，那么你在某些区域会得到好的近似，在其他区域则不然，因为你还没有读取那些数据。你也可以随机顺序遍历数据，最终可能得到一个相当不错的后验分布近似。

一种理解方式是，这本质上就像“忘记”了磁盘上还有其他数据，只是不断加载数据块，而不关心其他部分。显然，这种方法的成本最多是数据量的线性关系，加上一个与我们已查看数据量相关的二次项。

这听起来可能有点平淡，因为你为什么要把所有数据留在磁盘上呢？也许你想做更聪明的事。但有趣的是，这确实是可行的。

线性代数的高级洞见

这里有一个更深入的洞见，它在论坛上引发了一些讨论。当你在概率机器学习课程中学习时，你会知道要求解高斯过程回归，需要解决两个问题：求后验均值和后验协方差。如果你有中间矩阵的Cholesky分解，你就可以高效地解决这两个问题。所以，某种意义上你只需要做一件事：计算Cholesky分解。

然后，有人可能会说Cholesky太贵了，对于大数据集，你应该使用迭代线性求解器，比如我们在第3讲讨论的共轭梯度法。但共轭梯度法在教科书中通常被呈现为一种解决二次优化问题的方法，而不是矩阵分解。

如果你使用这种方法，似乎这里有两个不同的优化问题（Ax1 = b1 和 Ax2 = b2），你需要运行两次算法。实际上，共轭梯度法和Cholesky一样，也是一种在估计特定线性方程组解的同时，也估计矩阵逆的方法。只是Cholesky的性质在谈论它试图解决的那个线性问题的解时特别好。

这意味着我们实际上可以使用共轭梯度法同时完成这两件事。我们只需在其中一个问题上运行它，然后利用它构造的数字来直接估计一个逆矩阵，进而用于估计后验不确定性。

当我们这样做时，会出现一个非常有趣的结构：这个算法同时在做两件事。它从数据中学习，最终给出关于真实函数的有限不确定性；同时，它也在学习这个有限数据集上矩阵的逆。如果你在早期停止，这两者共同告诉你，你应该对整个问题有多不确定。因为你对最终试图估计的函数的误差（图中黑色实线和虚线之间的误差）来自两个不同的不确定性源：一是只有有限数据集，二是只进行了有限次计算（可能少于整个数据集）。

有趣的是，这两者可以同时被追踪，使用的算法看起来很像Cholesky，但更通用。如果策略是按输入空间位置逐个加载数据，你就得到类似Cholesky的算法。如果使用更智能的策略，你可能会收敛得更快，但代价是：智能策略需要规划它要探索数据的所有方向，这意味着该方法的每次迭代成本将是数据集大小的二次方，而不是Cholesky的线性成本。

线性代数部分总结

本节课我们一起学习了高斯过程推理这个基础案例的实际成本，远比教科书上的陈述要微妙。我们发现了非常细致的层次划分：

  1. 如果你决定永远不看整个数据集：那么你只需加载K个数据点，计算估计的成本是O(K³)。
  2. 如果你决定加载整个数据集，但只进行K次迭代：你可以构造像Cholesky这样的算法，成本是数据量的线性关系加上迭代次数的二次方，并伴随一个有意义的、能追踪已加载数据信息量的不确定性。
  3. 如果你想要一个能高效加载正确数据的算法以快速收敛：那么每次迭代会更昂贵，因为每一步都需要查看所有其他数据点来决定加载哪些以及加载多少，每次迭代成本是二次的。
  4. 最后，你也可以直接做完整的Cholesky分解：不关心迭代，直接加载全部数据并计算，成本是O(n³)。

为什么一定要做出查看整个数据集的强硬决定呢？毕竟之后你仍然会因为只有有限数据集而具有有限的不确定性。数据集的大小通常与你试图解决的问题没有直接关系。因此，当处理非平凡规模的数据时，思考如何操作数据、以何种顺序加载、在每次迭代中投入多少时间和计算资源是有意义的。

这引出了下一个场景：在某种意义上，我们在一个无限维的数据集上操作，我们需要决定评估的频率，并谨慎处理这些信息。

模拟与微分方程：作为滤波的求解器 ⏱️

上一节我们探讨了有限数据集上的计算。本节中我们来看看对随时间演化的系统进行模拟，这通常与常微分方程相关。

对于这类系统，Jonathan在第5讲中告诉我们，就像高斯过程是线性问题的推断机制一样，对于线性时变或时不变问题（即状态演化是先前状态的线性函数，观测是状态空间的线性变换，且观测具有高斯噪声），也有一个算法，称为卡尔曼滤波与平滑。这实际上是高斯过程推断的一个特例。

这个算法非常容易写下来，并且成本是时间步数的线性关系（尽管涉及矩阵求逆，所以状态空间元素数量的成本最多是立方的，我们可以应用之前线性代数的所有工具）。在数值计算中，我们经常通过嵌套算法来构建复杂问题的解决方案。因此，了解计算层次结构的底层是有用的，可以用来加速更高层次的计算。

直观上，滤波算法从链的左端开始，带着对初始状态的信念，随时间向前工作，遇到数据点时将其纳入。平滑算法则是一种簿记算法，确保链中所有早期变量与未来所做的观测保持一致。

然后，Nathanael在第7讲中指出，我们可以将这类算法应用于观测不是状态空间的线性高斯变换、状态转移也不是线性高斯函数的场景。这就是扩展卡尔曼滤波器，它通过在任何需要的地方进行泰勒近似（计算非线性函数的值及其雅可比矩阵）来“假装”关心的函数是线性的。

我们利用这个机制，将关于状态空间内状态之间非线性关系的信息（特别是当状态空间由函数值及其导数组成时）纳入进来。这些关系就是微分方程。

因此，事实证明，存在一种看起来、工作起来、表现起来都像滤波器的算法，可以求解微分方程。Nathanael的第7讲展示了这一点（第6讲是Nathanael介绍经典ODE求解器如何工作的插曲）。经典方法（如Runge-Kutta）是写在80年代的、性能极佳、非常稳健的“黑箱”代码，包含了无人敢改的魔法数字（Butcher表）。而这里介绍的滤波器框架则超级灵活，它是一种非常“计算机科学”的模拟方式。

我们首先写下特定的线性高斯系统（一个状态演化的马尔可夫链），然后通过说“在多个时间点上，我观测到状态之间的某些代数关系成立”来纳入关于微分方程的知识。例如，如果我说x在时间t的一阶导数与其零阶导数的某个函数之差为0，那就编码了一个常微分方程。

实际上，我可以纳入任何代数隐式方程，不仅仅是ODE，还有代数方程、连续对称性、哈密顿量，甚至偏微分方程。我们还可以纳入其他信息，例如在多个时间点观测到系统的测量值。这意味着我不需要知道非常复杂的工具箱（Runge-Kutta、几何积分器等）以及何时调用哪个，我只需要使用卡尔曼滤波器即可。

特别是，如果状态空间X中有某些部分没有出现在方程中，我可以使用其他类型的观测来告知算法这些值可能是什么。我们在COVID设置的例子中看到了这一点：如果动力系统中存在一些潜在力（如未知的接触率），那么观测系统轨迹并结合微分方程可以告诉你潜在力是什么（当然，取决于微分方程和潜在力与观测轨迹的交互方式）。

这里的要点是，我们非常明确地将信息视为计算的核心对象。这些方法所做的是管理信息。我们称这些项为信息算子，它们是更新你对动力系统所知知识的算法部分。信息是来自磁盘、传感器还是程序员写下的代数方程，这并不重要，它只是一个数字。信息算子也充当了算法用户和设计者之间的接口。

偏微分方程：高斯过程推断的视角 🌊

上一节我们看到了如何用滤波方法解决常微分方程。本节中我们来看看偏微分方程。

求解偏微分方程本质上是相同的事情。我们有一个方程，通常写作 D[V(x, t)] = f(V(x, t)) 的形式，其中D是包含许多偏导数的微分算子，V是我们关心的解（依赖于空间和时间的函数）。天气预报方程就长这样。

这其实就是另一个形如 G(x) = 0 的方程，只不过现在状态X是随时间变化的函数V。因此，我们必须更小心地进行计算。滤波器不能立即解决所有问题，因为如果我们对空间进行网格离散，然后在时间上使用滤波器求解，当网格点很多时，卡尔曼滤波器中的大逆矩阵仍然会出现并使一切变慢。

Marvin指出，在一个特殊情况下，即当这里的算子是线性的时候，这可以作为高斯过程回归的一个实例来解决。我们只是观测了函数V在各种不同算子下的线性投影。这本质上就是高斯过程推断，只不过对象是函数。函数在某种程度上与向量不同，函数空间可能很复杂。Marvin用一张包含复杂定理的幻灯片提供了“许可”，允许我们在某些前提下将函数视为具有无限多项的超长向量。

因此，你可以将此类微分方程视为只是在进行高斯过程推断，只不过有一个非常有趣的观测模型。如果它们不是线性的呢？目前还没有完整的答案，但你可以想象，就像ODE情况一样，尽可能地进行线性化。

你可以写下函数的生成模型（一个高斯过程先验），然后用各种信息源（如某个偏微分方程成立、边界值、传感器测量值）来告知它。所有这些信息算子组合在一起，形成了一个大型高斯过程推断方案，最终归结为一个大型最小二乘线性代数问题。然后，我们又回到了第2讲和第3讲，重新进行线性代数计算，最终得到关于这个动力系统的、带有不确定性量化的表示。

所以，未来当有人在机器学习中提出微分方程时，你不必再害怕它们。它们只是一种信息源。从观测一个函数，到观测一个函数的投影，这只是前进了一步。有时它是非线性的，但这不应该成为问题，因为我们在分类等任务中已经经常非线性地观测函数了。

积分：从蒙特卡洛到贝叶斯求积 📊

上一节我们探讨了微分方程的求解。本节中我们来看看机器学习中的积分问题，特别是在概率推断中。

我们在第9讲开始讨论积分，通常人们会想到蒙特卡洛方法。蒙特卡洛是一套优美的算法，它通过从概率分布P中抽取随机数并计算和，将确定性的推断问题转化为随机的统计问题。这种方法的好处是它给出一个无偏估计，并且适用于几乎任何积分问题。但它的代价是收敛速度慢：积分估计值收敛到真值的最佳速率是 O(1/√N)，其中N是样本数。

在第10讲（圣诞节前），我们研究了另一种截然不同的方法：不是对问题一无所知地随机抽样，而是花尽可能多的时间来精确建模被积函数，然后从精心选择的函数评估中学习其积分。这连接了贝叶斯求积的思想。

我们从一个高斯过程先验开始（现在这对你们来说很自然了），它不仅定义了函数，还通过核函数的可积性，定义了关于积分的高斯边际分布。当我们用函数在某些点上的评估值来条件化这个分布时，我们不仅得到了函数的估计，也得到了积分的估计。

这种方法要求我们能够计算核函数及其组合的积分。除此之外，它是一个允许我们进行积分的算法。我们看到，在某些设置下，特别是低维问题中，这类算法可以表现得非常好，收敛速度可以比蒙特卡洛快得多（多项式快，而非平方根快），甚至存在收敛速度超多项式的特殊算法类。

诚然，这主要适用于低维问题，将其扩展到高维问题并保持快速收敛是棘手的。构建在高维上与蒙特卡洛效果相当的版本并不难，但那样做意义不大。

这部分课程的主要收获不是“如果你需要积分，就用贝叶斯求积”，而更像是一个哲学观察：当你使用一个适用于所有问题的算法时，它通常在几乎所有问题上都表现不佳。你为通用性付出了代价。而如果你构建一个只适用于很小一类问题的算法，那么它可能在这类问题的每个实例上都表现得非常好，但一旦超出该类，就可能完全失效。

这让人联想到贝叶斯推断：如果一个先验将大量概率质量放在一个小的约束空间上，那么你就不需要太多数据在该空间内进行推断。但如果先验假设是错误的，算法可能完全失效。

这对计算意味着什么？这意味着，从构建算法的理论家的社会角度看，他们受学术体系激励去构建适用于大类问题的方法。但在实践中，如果你是为一个商业或科学上有价值的问题构建解决方案，你通常只有一个问题要解决。因此，如果你了解情况并构建一个针对你的问题的定制方法，有时可以提高性能。

深度学习优化：未解的挑战与机遇 🤖

圣诞节后，我们回归课程，步调有所改变，开始探讨当代机器学习中那些尚无完美解决方案、但又极其重要的数值问题。训练深度神经网络就是最主要的例子。

Frank在圣诞后的讲座中向你们展示了超过120种优化器的列表，并论证了它们中的大多数效果都差不多，你大可以使用Adam来训练你的深度神经网络。

这乍一听似乎问题已经解决了。但实际上，这些算法从根本上令人沮丧。它们真的很难用，需要“精心照料”、超参数调优、不一定总是有效、需要良好的网络初始化、你甚至需要决定何时停止训练。这与不久以前优化领域研究的情况截然不同。那时，优化讲座会告诉你，有现成的完美算法可以下载使用，你只需以正确的形式写出问题，然后它就能解决。

而现在，在机器学习中，我们面对的是另一个极端。我们有上百人的工程团队训练大语言模型，每周开会决定是否要降低学习率、改变批次大小，或者因为梯度爆炸而回滚到上周的安全状态。我们正处于一个非常有趣的时间点：这些方法并不真正“有效”，它们只是“足够好”到能引起公众注意并获得大量资金，但这并不是对资源的高效利用。

为什么会这样？可能的原因有三个：

  1. 问题规模大：大语言模型有数十亿参数。但这并非根本原因，因为优化在30-40年前就已处理百万维问题，而且芯片变得更强大。
  2. 问题非凸：经典优化理论主要适用于凸问题。但早在70年代，人们就已将拟牛顿法应用于非凸问题并使其工作。
  3. 随机性：这才是问题的核心，也回到了本课程关于计算不确定性的中心。

训练深度神经网络时，你并不是在最小化总体风险（因为数据有限），甚至也不是在最小化经验风险（因为计算也有限）。你实际上是在评估小批量梯度，这是一个围绕批次平均梯度分布的对象。这极大地降低了优化算法的性能，因为它引入了非常显著的噪声，实际上噪声可能比信号还大。

典型的批次平均梯度甚至不指向正确的方向。我们不能将其视为小扰动，因为你想计算的东西和你实际计算到的东西之间的扰动，相对于你想计算的东西来说并不小。你真正获得的是一个具有少量结构的随机比特源，游戏的目标是从这些随机变量中获取尽可能多的信息。

Frank在他的讲座中强调了这一点：我们不必对这个设置感到绝望。事实上，它提供了许多有趣的新途径。因为，我们不仅能访问那个红色箭头（平均梯度），还能访问所有金色箭头（每个样本的梯度），它们在训练过程中已经被计算出来了。典型的深度学习库只返回平均梯度，但如果你愿意，你可以访问所有这些信息。突然之间，你有了一个更丰富的信息源——一个数组，而不仅仅是一个向量。你还可以计算这个对象的梯度（曲率估计），使用自动微分。

随后，Michael在第12讲中指出，你或许可以利用这些曲率估计来构建有趣的算法，比如经典的基于牛顿的优化方法。只是需要小心，因为我们现在处于这个“深度”设置中，之前提到的三个问题都适用：模型非常大（曲率矩阵更大）、问题非凸（Hessian不正定）、一切都是随机的。

目前，整个社区（不仅仅是这个讲堂）对此还没有结论性的答案。我个人的预测是，当这样的答案出现时，将会非常激动人心。你会知道它出现了，因为它将使训练深度神经网络的速度难以置信地加快，带来数量级上的性能变化，从而改变我们进行机器学习的方式。

与此同时，我们也可以利用这些量来做训练之外的其他事情。例如，一旦网络训练好，我们想知道它到底知道什么，权重空间的哪些部分受到约束，模型有多不确定。我（基于Ats的工作）做了一次讲座，展示了我们可以使用这些曲率估计，通过拉普拉斯近似的思想，以各种轻量级方式为深度神经网络构建不确定性估计。

其思想是：你用你喜欢的方式训练深度神经网络，直到在权重空间找到一个点（即你的最佳猜测）。然后，在该点计算损失函数相对于网络权重的曲率近似。这告诉你在哪个方向上损失是平坦的（在训练数据上），在哪个方向上是陡峭的。平坦的方向意味着不确定性高，陡峭的方向意味着置信度高。我们通过取这个曲率矩阵的负逆，将其转化为一个高斯协方差矩阵，然后可以将其插入到损失函数的简单线性化中，以获得预测不确定性。

这种功能可以用于各种不同的用例，特别是可以与权重空间中的网络线性化相结合，得到一个机制，将任何深度神经网络近似地转化为一个参数化的高斯回归算法后验，仅使用线性代数和自动微分，无需花哨、复杂、昂贵的蒙特卡洛方法，并且可以在训练后完成而不影响算法性能。

总结：计算即推断 🧠

本节课我们一起回顾了整个课程，从线性代数到深度学习优化。现在，让我们来总结一下核心的高级要点。

  • 计算即推断：所有的计算都是信息的管理和操作。
  • 信息管理：在芯片上，通过从磁盘加载数据，或者由人类设计者写下可以在不同点评估的自然定律来获取信息。无论数字来自磁盘还是作为信息算子可用，算法都必须决定如何处理它们。
  • 主动代理：数值计算是智能体与数据源交互的行为。在求解微分方程时，算法必须决定在哪里离散化问题、进行多少次计算，这很明显是主动决策。对于数据也是如此：进行高斯过程回归或深度学习时，思考实际想使用多少数据、加载多少、在其上进行多少计算是有益的。
  • 新模式与新工具：如果我们认真对待这种联系，就可以构建新的、以数据为中心的计算方式。如果做得对，它们可能更灵活、更易用、更容易推广到不同设置，并提供如信息算子、高斯过程回归、卡尔曼滤波与平滑、拉普拉斯近似等设计模式，可以非常灵活地用于纳入来自各种来源的信息。
  • 数值计算是核心：最终，这一切都是数值计算。对于构建学习机器的人来说，机器学习就是数值计算。因此，了解数值算法（无论是经典的还是我们重新诠释的）如何工作，能让你成为更好的机器学习工程师。如果不了解，你就只能受限于别人为你制造的工具。

希望这次回顾对你们有所帮助。请记得填写最终反馈表。当然，我们还需要谈谈考试的事情。以上就是本次讲座的全部内容。祝大家考试顺利！

posted @ 2026-03-26 13:19  布客飞龙V  阅读(0)  评论(0)    收藏  举报