深度学习与神经网络10

2.9 图与集合

2.9.1 GNN 图神经网络：GCN、GAT、消息传递、图池化

2.9.1 GNN

一句话

每个节点反复从邻居那里收集信息、更新自己的表示；同一套更新函数用在所有节点上，所以能处理任意大小、任意连接方式的图。

消息传递框架（Gilmer 2017，几乎所有 GNN 都能写成这个形式）：

$\mathbf{h}_v^{(l+1)} = \phi\left(\mathbf{h}_v^{(l)},\ \bigoplus_{u\in\mathcal{N}(v)}\psi\left(\mathbf{h}_v^{(l)}, \mathbf{h}_u^{(l)}, \mathbf{e}_{uv}\right)\right)$

ψ 为每条边算一条「消息」，⊕ 是对邻居消息的聚合（求和、平均或最大值，与邻居顺序无关），φ 用聚合结果更新节点。堆 L 层，每个节点就能看到 L 跳以内的信息。

GCN（Kipf & Welling 2017）：消息就是邻居特征的线性变换，按度数归一化后求和：

$H^{(l+1)} = \sigma\left(\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}H^{(l)}W^{(l)}\right),\qquad \tilde{A} = A + I$

GAT（Veličković 2018）：邻居的权重不再由度数固定，用注意力算出来：

$\alpha_{vu} = \operatorname{softmax}_{u\in\mathcal{N}(v)}\left(\operatorname{LeakyReLU}\left(\mathbf{a}^\top[W\mathbf{h}_v\,\Vert\,W\mathbf{h}_u]\right)\right),\qquad \mathbf{h}'_v = \sigma\Big(\sum_{u}\alpha_{vu}W\mathbf{h}_u\Big)$

图池化（读出）：做整图级别的任务（分子性质预测、图分类）时，要把所有节点的表示汇成一个图向量。最简单的是全局求和 / 平均 / 最大池化；层级池化（DiffPool、TopK 池化）在中间层逐步把节点合并或删减，得到由细到粗的多层图表示。

几个常见问题

过平滑：层数多了，所有节点的表示趋于一致，所以 GCN、GAT 这类标准消息传递网络通常只用 2–4 层（经验范围，加残差、归一化等技巧可以更深）
表达能力上限：标准消息传递 GNN 区分不同图结构的能力不超过 1-WL 图同构测试（高阶 GNN、加子图特征等变体可以超过），有些结构不同的图在它看来一模一样
长距离依赖：两个节点相隔很多跳时，信息要传很多层才能到；图 Transformer 用全局注意力补这一点
几何等变：分子、物理系统里节点有三维坐标，等变 GNN（EGNN 等）保证整体旋转平移时输出跟着正确变换

import torch
import torch.nn as nn
import torch.nn.functional as F

def normalize_adj(A):
    A = A + torch.eye(A.size(0))                       # 加自环：节点也看自己
    d = A.sum(1).pow(-0.5)
    return d[:, None] * A * d[None, :]                 # D^-1/2 (A+I) D^-1/2

class GCNLayer(nn.Module):
    def __init__(self, din, dout):
        super().__init__()
        self.lin = nn.Linear(din, dout)
    def forward(self, x, A_norm):                      # x: (N, din)
        return A_norm @ self.lin(x)                    # 每个节点 = 邻居（含自己）特征按对称度数归一化系数的加权求和（每行系数之和一般不为 1）；偏置也一起被聚合了，和标准写法 Â X W + b 略有不同

class GATLayer(nn.Module):
    def __init__(self, din, dout, heads=4):
        super().__init__()
        self.W = nn.Linear(din, heads * dout, bias=False)
        self.a_src = nn.Parameter(torch.randn(heads, dout) * 0.1)
        self.a_dst = nn.Parameter(torch.randn(heads, dout) * 0.1)
        self.h, self.d = heads, dout
    def forward(self, x, A):                           # A: (N, N) 邻接矩阵（0/1）
        N = x.size(0)
        h = self.W(x).view(N, self.h, self.d)          # (N, 头数, d)
        e = (h * self.a_dst).sum(-1)[:, None] + (h * self.a_src).sum(-1)[None, :]   # (N, N, 头数)
        mask = (A + torch.eye(N)) == 0
        e = F.leaky_relu(e, 0.2).masked_fill(mask[..., None], float("-inf"))        # 只在邻居之间算注意力
        alpha = torch.softmax(e, dim=1)                # 对每个节点的邻居归一化
        return torch.einsum("vuh,uhd->vhd", alpha, h).reshape(N, -1)   # 多头拼接

class MessagePassing(nn.Module):
    """稀疏边表写法：只遍历存在的边，适合大图（PyG 的 MessagePassing 就是这个模式）"""
    def __init__(self, din, dout):
        super().__init__()
        self.msg, self.upd = nn.Linear(2 * din, dout), nn.Linear(din + dout, dout)
    def forward(self, x, edge_index):                  # edge_index: (2, E)，每列一条边 src → dst
        src, dst = edge_index
        m = F.relu(self.msg(torch.cat([x[dst], x[src]], -1)))            # 每条边一条消息
        agg = torch.zeros(x.size(0), m.size(1)).index_add_(0, dst, m)    # 汇总到目标节点
        return F.relu(self.upd(torch.cat([x, agg], -1)))

# 合成数据：3 个社区，社区内连边多、社区间少；节点特征带强噪声，单看自己很难分
torch.manual_seed(0)
N, K = 300, 3
labels = torch.randint(0, K, (N,))
same = labels[:, None] == labels[None, :]
A = (torch.rand(N, N) < torch.where(same, torch.tensor(0.05), torch.tensor(0.005))).float()
A = torch.triu(A, 1); A = A + A.T                      # 无向图
X = F.one_hot(labels, K).float() + 1.5 * torch.randn(N, K)
train = torch.rand(N) < 0.1                            # 只有 10% 的节点有标签

def fit(params, forward):
    opt = torch.optim.Adam(params, lr=0.01, weight_decay=5e-4)
    for _ in range(200):
        loss = F.cross_entropy(forward()[train], labels[train])
        opt.zero_grad(); loss.backward(); opt.step()
    return (forward().argmax(1)[~train] == labels[~train]).float().mean().item()

A_norm = normalize_adj(A)
mlp = nn.Sequential(nn.Linear(K, 16), nn.ReLU(), nn.Linear(16, K))
g1, g2 = GCNLayer(K, 16), GCNLayer(16, K)
gat1, gat2 = GATLayer(K, 8, heads=4), GATLayer(32, K, heads=1)
print("MLP（只看自己）:", round(fit(mlp.parameters(), lambda: mlp(X)), 3))
print("GCN（看邻居）  :", round(fit(list(g1.parameters()) + list(g2.parameters()),
                                    lambda: g2(F.relu(g1(X, A_norm)), A_norm)), 3))
print("GAT（注意力加权）:", round(fit(list(gat1.parameters()) + list(gat2.parameters()),
                                      lambda: gat2(F.elu(gat1(X, A)), A)), 3))

# 稀疏边表写法 + 图池化（读出）：真正的多图 batch 用 batch 向量记录每个节点属于哪张图；这里只有一张图，借同样的写法把节点分成 3 组分别读出
edge_index = A.nonzero().T                             # (2, E)
h = MessagePassing(K, 16)(X, edge_index)               # (300, 16)
batch = torch.arange(N) // 100                         # 把这张图的 300 个节点按编号分成 3 组分别读出（组间有边，不等于 3 张独立的图；真正的多图 batch 图之间不能有边）
counts = torch.bincount(batch).float()[:, None]
graph_mean = torch.zeros(3, 16).index_add_(0, batch, h) / counts     # 全局平均池化
graph_max = torch.full((3, 16), -1e9).scatter_reduce(0, batch[:, None].expand(-1, 16), h, "amax")
print("3 个节点组的读出表示:", tuple(graph_mean.shape), tuple(graph_max.shape))

常见说法与实测

常见说法：GAT 用注意力给不同邻居不同的权重，比按度数固定加权的 GCN 更灵活，效果更好。

实测：这个三社区的合成图上，GCN 0.86、GAT 0.77（只看节点自己的 MLP 0.39）。另外换 5 个随机种子：GCN 平均 0.853，GAT 平均 0.776，GAT 只在 1 个种子上略好。

怎么理解：一种可能的解释：这张图里同一社区的邻居都一样有用，按度数归一化后求和已经是合适的聚合方式，注意力能多做的事不多，反而多了参数；有标签的节点只有约 30 个，可能更容易过拟合。训练预算、学习率等也会有影响，这里没有做消融逐一排除。GAT 的优势更可能在邻居的重要性确实不同的图上体现。在标准基准上公平调参之后，GCN 这类简单模型也经常和更复杂的 GNN 持平（Shchur 等 2018，《Pitfalls of Graph Neural Network Evaluation》）。

典型例子

分子性质与药物发现：原子是节点、化学键是边；2020 年用 MPNN 筛选出了新抗生素 halicin
推荐系统：Pinterest 的 PinSage 在数十亿节点的图上做推荐
交通：Google 地图的到达时间预测用了 GNN
物理仿真与天气：DeepMind 的 GNS 学习粒子动力学，MeshGraphNets 学网格仿真，GraphCast 做全球天气预报
机器人：多体系统、运动学链、场景中物体之间的关系；不规则排布的触觉单元（taxel）阵列可以建成图，相邻传感点连边

优势：直接处理不规则连接结构；对节点编号顺序置换等变；参数量与图的大小无关，能泛化到不同大小的图；能把「谁和谁相连」这种关系先验写进结构
局限：过平滑，层数受限；远距离依赖要很多跳；表达能力受 1-WL 上限约束；大图需要邻居采样和专门的批处理
适合的数据：图结构数据：分子、社交与引用网络、交通路网、仿真网格、知识图谱、机器人运动学链、不规则排布的传感器阵列

2.9.2 PointNet 与 DeepSets

一句话

每个元素（点）独立地过同一个 MLP，再用一个与顺序无关的池化（求和或取最大值）汇总；输入是一个无序集合时，输出不受元素顺序影响。

DeepSets 定理（Zaheer 2017）：在一定条件下，「对集合的置换不变函数」都可以写成

$f(\{\mathbf{x}_1,\dots,\mathbf{x}_n\}) = \rho\left(\sum_{i=1}^{n}\phi(\mathbf{x}_i)\right)$

条件是这样的：元素取自可数集时成立；元素是连续的标量、集合大小最多为 M 时，Wagstaff 2019 证明求和之后的隐空间要有 M 维，才能表示所有连续的置换不变函数。点云这类向量元素也常用这个结构，但这个维度阈值不能直接照搬；而且这说的是表示能力，维度够了也不保证训练能学会。

PointNet（Qi 2017） 用的是取最大值：

$f = \rho\left(\max_{i}\ \phi(\mathbf{x}_i)\right)$

每个特征维度取所有点里的最大值，相当于每个维度「挑一个最有代表性的点」。原版还有一个小网络 T-Net 先预测一个变换矩阵把点云对齐。

PointNet 的缺点和后续：每个点独立处理，没有利用邻域的局部几何。PointNet++ 先采样中心点、在每个中心点的邻域里用小 PointNet 提特征，逐层抽象；DGCNN 在特征空间建 k 近邻图做边卷积；Point Transformer 在邻域内做注意力。

import torch
import torch.nn as nn

class PointNet(nn.Module):
    def __init__(self, n_cls, d=128):
        super().__init__()
        self.phi = nn.Sequential(nn.Linear(3, 64), nn.ReLU(), nn.Linear(64, d), nn.ReLU())   # 每个点共享
        self.rho = nn.Sequential(nn.Linear(d, 64), nn.ReLU(), nn.Linear(64, n_cls))

    def forward(self, pts):                       # pts: (B, N, 3)，N 个点，没有顺序
        feat = self.phi(pts)                      # (B, N, d) 逐点特征
        return self.rho(feat.max(dim=1).values)   # 对称池化：与点的顺序无关

def sample_shapes(n, n_pts=256):
    """球面、立方体表面、圆柱面三类点云，随机缩放"""
    labels = torch.randint(0, 3, (n,))
    p = torch.randn(n, n_pts, 3)
    sphere = p / p.norm(dim=-1, keepdim=True)
    cube = torch.rand(n, n_pts, 3) * 2 - 1
    face = torch.randint(0, 3, (n, n_pts, 1))
    cube.scatter_(2, face, torch.sign(torch.randn(n, n_pts, 1)))       # 把某一个坐标推到 ±1 → 落在表面上
    ang = torch.rand(n, n_pts) * 2 * torch.pi
    cyl = torch.stack([ang.cos(), ang.sin(), torch.rand(n, n_pts) * 2 - 1], -1)
    shapes = torch.stack([sphere, cube, cyl], 1)[torch.arange(n), labels]
    return shapes * (0.5 + torch.rand(n, 1, 1)), labels

torch.manual_seed(0)
model = PointNet(3)
opt = torch.optim.Adam(model.parameters(), lr=1e-3)
for step in range(300):
    x, y = sample_shapes(32)
    loss = nn.functional.cross_entropy(model(x), y)
    opt.zero_grad(); loss.backward(); opt.step()

x, y = sample_shapes(300)
with torch.no_grad():
    print("形状分类准确率:", (model(x).argmax(1) == y).float().mean().item())
    perm = torch.randperm(x.size(1))
    print("打乱点的顺序，输出不变:", torch.allclose(model(x), model(x[:, perm]), atol=1e-5))
    print("只给一半的点也能分类:", (model(x[:, :128]).argmax(1) == y).float().mean().item())

典型例子

三维物体分类与分割：ModelNet40、ShapeNet
自动驾驶：激光雷达点云的检测和分割
机器人抓取：从深度相机点云预测抓取位姿（Contact-GraspNet 基于 PointNet++）；3D Diffusion Policy（DP3）用一个轻量的点云编码器给扩散策略提供条件
无序集合：一组检测框、一群粒子（高能物理的喷注）、多智能体的状态集合、稀疏的接触点

优势：输入顺序无关、点数可变，直接吃原始点云，不需要体素化；逐点共享参数，快；最大池化只依赖每个维度上取到最大值的那几个点，缺的不是这几个点时输出不变
局限：原版不利用局部邻域结构；对整体旋转没有不变性（靠数据增强或等变网络）
适合的数据：集合与点云：激光雷达、深度相机点云、粒子集合、无序的检测结果、多智能体状态、接触点集合

2.10 不确定性与多峰输出

2.10.1 深度集成与 MC Dropout（估计预测的不确定度）

2.10.1 深度集成与 MC Dropout

一句话

不只给一个预测，还给「有多确定」：多个模型（或同一个模型开着 Dropout 多次前向）的预测分歧越大，通常说明这里越缺数据；分歧也受模型结构、训练和近似方式影响，远离数据时不保证变大。

两类不确定性

类型	来源	能不能靠更多数据消除	怎么估
认知不确定性（epistemic）	模型没见过这一片输入	能	贝叶斯神经网络、集成、MC Dropout
偶然不确定性（aleatoric）	数据本身的噪声或多解	不能	让网络输出方差或一个分布（2.10.2）

贝叶斯神经网络：每个权重用一个分布 p(θ | 数据) 表示，预测是对所有可能权重的平均：

$p(y\mid x, \mathcal{D}) = \int p(y\mid x,\theta)\,p(\theta\mid\mathcal{D})\,d\theta \approx \frac{1}{M}\sum_{m=1}^{M}p(y\mid x,\theta_m)$

θ_m 是从后验里采的样本。真实的后验算不出来，下面三种方法都能给出不确定度，但和贝叶斯的关系各不相同：

方法	做法	成本
深度集成（Lakshminarayanan 2017）	不同随机初始化训练 M 个模型（常用 5 个），用预测的分歧当不确定度。原论文明确把它作为非贝叶斯方法提出：M 个模型不是从后验里采的样本，和上面的平均只是形式相似	训练和存储 M 份；效果通常最好
MC Dropout（Gal & Ghahramani 2016）	推理时不关 Dropout，前向 T 次；可以解释成对后验的一种变分近似	几乎零额外训练成本
变分推断（Bayes by Backprop，Blundell 2015）	每个权重学一个高斯分布，显式近似后验	参数翻倍，训练复杂

import torch
import torch.nn as nn

torch.manual_seed(0)
x = torch.cat([torch.rand(100, 1) * 2 - 3, torch.rand(100, 1) * 2 + 1])   # 训练数据只在 [-3,-1] 和 [1,3]
y = torch.sin(x) + 0.1 * torch.randn_like(x)

def make():
    return nn.Sequential(nn.Linear(1, 128), nn.ReLU(), nn.Dropout(0.1),
                         nn.Linear(128, 128), nn.ReLU(), nn.Dropout(0.1), nn.Linear(128, 1))

def train(model):
    opt = torch.optim.Adam(model.parameters(), lr=1e-2, weight_decay=1e-4)
    for _ in range(1000):
        loss = ((model(x) - y) ** 2).mean()
        opt.zero_grad(); loss.backward(); opt.step()
    return model

x_test = torch.tensor([[-2.0], [0.0], [2.0], [5.0]])       # 见过、中间空洞、见过、外推

model = train(make())
model.train()                                               # MC Dropout：推理时故意保持 Dropout 打开（这个模型没有 BatchNorm；有 BN 时只把 Dropout 层设成 train）
with torch.no_grad():
    preds = torch.stack([model(x_test) for _ in range(200)])
print("MC Dropout 的标准差（x = -2, 0, 2, 5）:", preds.std(0).squeeze().numpy().round(3))

ensemble = [train(make()).eval() for _ in range(5)]        # 深度集成：5 个随机初始化
with torch.no_grad():
    preds = torch.stack([m(x_test) for m in ensemble])
print("深度集成的标准差（x = -2, 0, 2, 5）:", preds.std(0).squeeze().numpy().round(3))
# 集成在外推处（x=5）的分歧明显变大；MC Dropout 在两段数据之间的空洞（x=0）处变大，
# 到了外推处反而变小。远离数据时 MC Dropout 给出的不确定度不可靠，这是它常被批评的一点

MC_Dropout与深度集成的不确定度

灰色区域里没有训练数据：左边是两段数据之间的空洞，右边是外推。深度集成的区间在空洞处略微变宽、在外推处明显变宽；MC Dropout 的区间处处都很窄，到了外推处也没有变宽，说明在这个例子里它远离数据时给出的不确定度不可靠。两者的均值在外推处都偏离了真实函数。图里的色带是预测之间的分歧（±2 个标准差），不是校准过的置信区间，也没把数据本身的噪声算进去。

常见说法与实测

常见说法：MC Dropout 近似贝叶斯神经网络，推理时多次前向，离训练数据越远，预测之间的分歧越大。

实测：训练数据只在 [−3, −1] 和 [1, 3]。MC Dropout 的标准差在两段数据之间的空洞（x = 0）处是 0.115，比有数据的地方（0.065–0.082）大；在外推处（x = 5）只有 0.031，比有数据的地方还小。深度集成的标准差在空洞处是 0.033（有数据处 0.017–0.027），外推处 0.073，随着远离数据而变大。

怎么理解：Dropout 的随机性只来自随机关掉的神经元，它并不知道哪里缺数据。这个例子里 MC Dropout 在外推处的预测几乎是一条平线（见上图），那里的输出主要由偏置和少数单元决定，随机关掉单元带来的变化很小。文献里对 MC Dropout 的不确定度质量也有不少批评（如 Osband 2016）。深度集成的几个模型来自不同的随机初始化，在没有数据约束的区域往往走向不同的函数，分歧变大，这是它在实践中通常更可靠的原因。这也不是保证：几个模型在某些远离数据的区域可能恰好外推得差不多，分歧就不会变大。

典型例子

模型预测控制与基于模型的强化学习：PETS 用概率集成学动力学模型，规划时考虑模型的不确定度，避免利用模型的错误
主动学习：优先标注模型最不确定的样本
安全攸关场景：医疗诊断、自动驾驶，在不确定时交给人或采取保守动作
分布外检测：输入和训练数据差别很大时报警
贝叶斯优化：用不确定度权衡探索与利用

优势：知道自己不知道什么；深度集成实现简单、效果稳定；MC Dropout 几乎免费
局限：真实后验难算，各种近似的校准质量不一；集成要训练和存储多份模型；推理要多次前向；在远离数据的地方，不确定度估计本身也可能不可靠
适合的数据：数据稀疏、可能遇到分布外输入、错误代价高的任务：医疗、自动驾驶、机器人安全控制、需要探索的强化学习

2.10.2 MDN 混合密度网络

2.10.2 MDN

一句话

网络输出一个高斯混合分布的参数（每个分量的权重、均值、方差），用负对数似然训练；同一个输入可以对应几个都合理的输出。

$p(\mathbf{y}\mid\mathbf{x}) = \sum_{k=1}^{K}\pi_k(\mathbf{x})\,\mathcal{N}\left(\mathbf{y}\mid\boldsymbol\mu_k(\mathbf{x}),\ \sigma_k^2(\mathbf{x})I\right),\qquad \mathcal{L} = -\log p(\mathbf{y}\mid\mathbf{x})$

π 用 softmax 保证和为 1，σ 用 exp 保证为正。Bishop 1994 年提出。

它解决的问题：1.4 节讲过，MSE 回归学到的是条件均值。逆问题（已知结果求原因）常常一对多：比如机械臂逆运动学，同一个末端位置可以由「肘部向上」和「肘部向下」两种关节角实现，取平均会得到一个够不到目标的姿态。MDN 用几个分量一起表示这个条件分布，学得好时不同的解落在不同的分量上。不过损失只看混合后的密度，多余的分量可能权重变小，也可能和别的分量重合、一起表示同一个峰，训练不保证每个分量正好对应一个解。

代价：分量数 K 要事先定；训练时 σ 可能塌到 0、似然趋于无穷，要对 σ 设下限；输出维度高时（比如一整段 16 步 × 7 维的动作）混合高斯不够灵活，这是扩散策略取代它的原因之一。

import torch
import torch.nn as nn
import torch.nn.functional as F

# 一对多的反问题：已知 x，求 y；y = +√x 和 y = −√x 都对（类似逆运动学里「肘部朝上」和「肘部朝下」两种解）
torch.manual_seed(0)
x = torch.rand(3000, 1)
branch = torch.where(torch.rand(3000, 1) < 0.5, -1.0, 1.0)
y = branch * x.sqrt() + 0.03 * torch.randn_like(x)

class MDN(nn.Module):
    def __init__(self, K=3, hidden=64):
        super().__init__()
        self.body = nn.Sequential(nn.Linear(1, hidden), nn.Tanh(), nn.Linear(hidden, hidden), nn.Tanh())
        self.pi, self.mu, self.log_sigma = nn.Linear(hidden, K), nn.Linear(hidden, K), nn.Linear(hidden, K)

    def forward(self, x):
        h = self.body(x)
        return F.log_softmax(self.pi(h), -1), self.mu(h), self.log_sigma(h).clamp(-5, 2)   # σ 设下限

def mdn_nll(log_pi, mu, log_sigma, y):
    log_prob = -0.5 * ((y - mu) / log_sigma.exp()) ** 2 - log_sigma - 0.5 * torch.log(torch.tensor(2 * torch.pi))
    return -torch.logsumexp(log_pi + log_prob, dim=-1).mean()       # 混合分布的负对数似然

mdn = MDN()
reg = nn.Sequential(nn.Linear(1, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1))
opt = torch.optim.Adam(list(mdn.parameters()) + list(reg.parameters()), lr=3e-3)
for step in range(3000):
    loss = mdn_nll(*mdn(x), y) + F.mse_loss(reg(x), y)
    opt.zero_grad(); loss.backward(); opt.step()

x0 = torch.tensor([[0.64]])            # 正确答案是 +0.8 和 -0.8
with torch.no_grad():
    log_pi, mu, log_sigma = mdn(x0)
    print(f"MSE 回归给出 y = {reg(x0).item():.3f}（两个解的平均，哪一个都不是）")
    for k in range(3):
        print(f"  MDN 分量 {k}：权重 {log_pi[0, k].exp():.2f}，均值 {mu[0, k]:.3f}，标准差 {log_sigma[0, k].exp():.3f}")
    k = torch.multinomial(log_pi.exp().expand(8, -1), 1).squeeze(1)  # 从混合分布采样：先选分量，再在分量里采
    samples = mu[0, k] + log_sigma[0, k].exp() * torch.randn(8)
    print("从 MDN 采样 8 次:", samples.numpy().round(2))

# 对照：Bishop 的经典反问题 x = y + 0.3·sin(2πy)，三个解离得不远、各自的概率也不同
yb = torch.rand(3000, 1)
xb = yb + 0.3 * torch.sin(2 * torch.pi * yb) + 0.05 * torch.randn_like(yb)
reg_b = nn.Sequential(nn.Linear(1, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1))
opt = torch.optim.Adam(reg_b.parameters(), lr=3e-3)
for step in range(3000):
    loss = F.mse_loss(reg_b(xb), yb)
    opt.zero_grad(); loss.backward(); opt.step()
g = torch.linspace(0, 1, 100001)
roots = g[1:][((g + 0.3 * torch.sin(2 * torch.pi * g) - 0.45).sign().diff() != 0)].round(decimals=3).unique()   # 无噪声方程在 x = 0.45 处的三个近似根：在 10 万格的网格上找变号位置，保留三位小数
with torch.no_grad():
    print(f"Bishop 反问题，x = 0.45：三个近似根 {roots.numpy().round(3)}，MSE 回归给出 {reg_b(torch.tensor([[0.45]])).item():.3f}")

MDN与回归

每个 x 对应两个都正确的 y（+√x 和 −√x）。MSE 回归学到的是两条分支的平均，几乎贴着 y = 0 走，哪一条分支都不在上面；MDN 的两个分量分别贴着两条分支。

常见说法与实测

常见说法：一对多的问题用 MSE 回归，会得到几个解的平均，落在解与解之间的空白处。

实测：y = ±√x 这种对称的两分支问题上确实如此：x = 0.64 时回归给出 −0.011，两个解是 ±0.8。但在 Bishop 的经典反问题上，x = 0.45 的三个解（无噪声方程的近似根）是 0.179、0.559、0.750，回归给出 0.537，紧挨着中间那个解。

怎么理解：回归输出的是按概率加权的平均。几个解对称、概率相近时，平均落在空白处；解的位置偏向一边、或者某个解的概率更大时，平均可能恰好落在某个解附近，看上去像答对了，它仍然只是一个平均值，输入稍有变化就可能落进空白。所以不能靠「回归结果看起来合理」来判断数据是不是多峰。

典型例子

手写生成：Graves 2013 用 LSTM + MDN 输出下一笔的位移分布，生成逼真的笔迹
World Models（Ha & Schmidhuber 2018）：MDN-RNN 预测下一时刻隐变量的混合分布
机器人模仿学习：robomimic 里的 BC-RNN-GMM 用混合高斯动作头处理多峰演示
自动驾驶轨迹预测：MultiPath 等输出几条可能轨迹及各自的概率
逆运动学、语音合成的声学参数生成

优势：一个输入可以对应多个合理输出，不会把多峰平均成无效值；给出完整的条件分布，包括不确定度；只改输出头和损失，主干不变
局限：分量数要预设；训练数值不稳定（σ 塌缩、分量坍缩）；高维输出时表达力有限
适合的数据：一对多的映射：逆运动学、轨迹预测、手写和语音参数生成、低维的多峰机器人动作

2.11 连续与隐式表示

2.11.1 NeRF 与隐式神经表示

一句话

用一个 MLP 表示一个连续的信号：输入坐标，输出该点的值（颜色、密度、距离）。一张图、一个三维场景就存成一小组网络权重。

隐式神经表示（INR）：图像是「(x, y) → 颜色」，三维形状是「(x, y, z) → 到表面的有符号距离」，音频是「t → 振幅」。用 MLP 拟合这个函数，可以在任意坐标上查询，存储量不随采样网格变大；但能恢复多少细节仍受数据、坐标编码的频率和网络容量限制，不等于无限分辨率。

频谱偏差：普通 MLP 天然偏向学低频，直接输入坐标会拟合出一张糊图。两种解决办法：

傅里叶特征 / 位置编码（NeRF、Tancik 2020）：先把坐标映射成一组不同频率的正余弦，再送进 MLP

$\gamma(p) = \left(\sin(2^0\pi p),\ \cos(2^0\pi p),\ \dots,\ \sin(2^{L-1}\pi p),\ \cos(2^{L-1}\pi p)\right)$

SIREN（Sitzmann 2020）：激活函数直接用 sin，配合特殊初始化

NeRF（Mildenhall 2020）：MLP 输入三维位置和观察方向，输出颜色 c 和体密度 σ。渲染一个像素时，沿着相机射线采样 N 个点，按体渲染公式累加：

$C(\mathbf{r}) = \sum_{i=1}^{N} T_i\left(1 - e^{-\sigma_i\delta_i}\right)\mathbf{c}_i,\qquad T_i = \exp\left(-\sum_{j<i}\sigma_j\delta_j\right)$

T_i 是光线走到第 i 个点时还剩多少没被挡住，δ_i 是相邻采样点的间距。整个渲染过程可微，所以只用多张已知位姿的照片、以「渲染出的像素和照片一致」为目标，就能把场景学进 MLP。之后的 Instant-NGP 用多分辨率哈希编码把训练从小时级压到秒级；3D Gaussian Splatting（2023）改用显式的高斯点云，渲染快得多。

import torch
import torch.nn as nn

# 目标「图像」：64×64 的高频图案；坐标 (x, y) ∈ [0,1]² → 灰度
H = 64
g = torch.linspace(0, 1, H)
coords = torch.stack(torch.meshgrid(g, g, indexing="ij"), -1).view(-1, 2)
target = (torch.sin(20 * coords[:, :1]) * torch.cos(15 * coords[:, 1:]) + torch.sin(40 * coords[:, :1] * coords[:, 1:])) / 2

def fourier_features(p, L=8):
    freqs = 2.0 ** torch.arange(L) * torch.pi             # π, 2π, 4π, …
    ang = p[..., None] * freqs                             # (N, 2, L)
    return torch.cat([ang.sin(), ang.cos()], -1).flatten(1)   # (N, 4L)

def fit(encode, in_dim, steps=1500):
    torch.manual_seed(0)
    mlp = nn.Sequential(nn.Linear(in_dim, 128), nn.ReLU(), nn.Linear(128, 128), nn.ReLU(), nn.Linear(128, 1))
    opt = torch.optim.Adam(mlp.parameters(), lr=1e-3)
    x = encode(coords)
    for _ in range(steps):
        loss = ((mlp(x) - target) ** 2).mean()
        opt.zero_grad(); loss.backward(); opt.step()
    return 10 * torch.log10(4 / loss).item()               # PSNR（取值范围 [-1, 1]，峰值差 2）

print(f"直接输入坐标：PSNR {fit(lambda p: p, 2):.1f} dB")
print(f"傅里叶特征  ：PSNR {fit(fourier_features, 32):.1f} dB")

def render_ray(sigmas, colors, deltas):
    """NeRF 的体渲染：沿一条光线把 N 个采样点的颜色按透明度加权求和"""
    alpha = 1 - torch.exp(-sigmas * deltas)                                    # 每段的不透明度
    T = torch.cumprod(torch.cat([torch.ones(1), 1 - alpha[:-1] + 1e-10]), 0)   # 光线走到第 i 点时的剩余透射率
    weights = T * alpha
    return (weights[:, None] * colors).sum(0), weights

# 一条光线上 64 个采样点：前面是空气，中间有一个红色的不透明物体，后面有一个被它挡住的蓝色物体
sig = torch.zeros(64); sig[30:34] = 50.0; sig[50:54] = 50.0
col = torch.zeros(64, 3); col[30:34] = torch.tensor([1., 0., 0.]); col[50:54] = torch.tensor([0., 0., 1.])
rgb, w = render_ray(sig, col, torch.full((64,), 0.05))
print("渲染出的颜色:", rgb.numpy().round(3), "（红色，蓝色物体被挡住）")

傅里叶特征与频谱偏差

典型例子

新视角合成：NeRF 从几十张照片重建场景，渲染任意视角
三维重建：神经有符号距离场（NeuS 等）重建高质量表面
机器人：用 NeRF 表示场景做抓取规划（包括透明物体）；Neural Contact Fields 用触觉信号推断物体与环境的接触位置
信号压缩：一张图、一段视频存成一个小网络（COIN 等）

优势：连续表示，可以在任意坐标查询；存储量不随采样网格变大（细节受容量限制）；完全可微，便于和渲染、物理过程结合求解逆问题
局限：每个场景要单独优化，原版 NeRF 训练要几个小时；渲染要对每条光线做很多次 MLP 前向，慢；推广到新场景需要额外的条件化设计
适合的数据：定义在连续坐标上的信号：图像、三维形状、辐射场、音频、物理场（温度、压力、接触场）

2.11.2 Neural ODE 神经常微分方程

2.11.2 Neural ODE

一句话

用神经网络定义隐状态随时间变化的速度 dh/dt = f(h, t)，输出由数值积分得到，离散的层变成了连续的积分时间；ResNet 是它的欧拉法离散化。

$\frac{d\mathbf{h}(t)}{dt} = f_\theta(\mathbf{h}(t), t),\qquad \mathbf{h}(t_1) = \mathbf{h}(t_0) + \int_{t_0}^{t_1} f_\theta(\mathbf{h}(t), t)\,dt$

和 ResNet 的关系：残差块 h_{l+1} = h_l + F(h_l) 就是步长为 1 的欧拉法。Neural ODE（Chen 2018，NeurIPS 最佳论文）把层数推到连续极限，用任意 ODE 求解器（如四阶龙格-库塔）积分。

几个特点

时间连续：可以在任意时刻取值，天然处理采样不均匀的时间序列
伴随法：反向传播时再解一个反向的 ODE 求梯度，理想情况下显存和积分步数无关（实际有数值误差，常要配合检查点）。下面的代码没用伴随法，是直接对每个积分子步反传，显存随步数增长
精度可调：求解器的容差就是「深度」，精度和速度可以在推理时权衡
局限：解是唯一的，同一时刻两条轨迹不能经过同一个状态，所以从 h(0) 到 h(T) 的映射连续可逆，有些映射表示不了（比如一维里把 −1 和 1 对调）；f 显式依赖 t 时，轨迹可以在不同时刻经过同一个状态，但这条限制仍在。系统变「刚性」时步数暴涨

相关：Latent ODE 在隐空间里做 ODE，处理不规则采样的医疗数据；Liquid Time-constant Networks（LTC）和它的闭式近似 CfC 是由 ODE 定义的循环单元，在无人机和自动驾驶的小规模控制任务里表现出较好的鲁棒性；连续正规化流（FFJORD）把 2.2.6 的标准化流做成连续时间版本。

import math
import torch
import torch.nn as nn

def rk4_solve(f, h0, ts, max_dt=0.1):
    """四阶龙格-库塔，沿任意（可以不等间隔）的时刻积分；整个求解过程直接可以反向传播"""
    hs, h = [h0], h0
    for t0, t1 in zip(ts[:-1], ts[1:]):
        n = max(1, math.ceil((t1 - t0).item() / max_dt))     # 间隔大时细分，保证精度
        dt = (t1 - t0) / n
        t = t0
        for _ in range(n):
            k1 = f(t, h); k2 = f(t + dt / 2, h + dt / 2 * k1)
            k3 = f(t + dt / 2, h + dt / 2 * k2); k4 = f(t + dt, h + dt * k3)
            h = h + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
            t = t + dt
        hs.append(h)
    return torch.stack(hs)

# 真实系统：阻尼单摆 θ'' = -sin θ - 0.2 θ'，状态 (θ, ω)
true_f = lambda t, s: torch.stack([s[..., 1], -torch.sin(s[..., 0]) - 0.2 * s[..., 1]], -1)
torch.manual_seed(0)
s0 = torch.tensor([2.0, 0.0])
fine_t = torch.linspace(0, 15, 151)
with torch.no_grad():
    truth = rk4_solve(true_f, s0, fine_t)
idx = torch.sort(torch.randperm(100)[:30]).values            # 只在前 10 秒里不等间隔地取 30 个观测
idx[0] = 0
ts, obs = fine_t[idx], truth[idx] + 0.01 * torch.randn(30, 2)

class ODEFunc(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(2, 64), nn.Tanh(), nn.Linear(64, 2))
    def forward(self, t, h):
        return self.net(h)                                    # 学习向量场 dh/dt

def evaluate(f):
    with torch.no_grad():
        pred = rk4_solve(f, obs[0], fine_t)
    return ((pred[:101] - truth[:101]) ** 2).mean().item(), ((pred[101:] - truth[101:]) ** 2).mean().item()

# 方式 1：从 0 秒一口气积分到 10 秒，和全部 30 个观测比
torch.manual_seed(0)
f = ODEFunc()
opt = torch.optim.Adam(f.parameters(), lr=5e-3)
for it in range(1, 1001):
    pred = rk4_solve(f, obs[0], ts, max_dt=0.25)
    loss = ((pred - obs) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()
    if it in (300, 1000):
        fit, ext = evaluate(f)
        print(f"整条轨迹一次积分，训练 {it:4d} 步：拟合区间误差 {fit:.4f}，外推区间（没有观测）误差 {ext:.4f}")

# 方式 2：每次从随机一个观测出发，只积分到后面 8 个观测（多次打靶法的思路）
torch.manual_seed(0)
f = ODEFunc()
opt = torch.optim.Adam(f.parameters(), lr=3e-3)
for it in range(3000):
    i = torch.randint(0, 30 - 8, (1,)).item()
    pred = rk4_solve(f, obs[i], ts[i:i + 9], max_dt=0.25)
    loss = ((pred - obs[i:i + 9]) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()
fit, ext = evaluate(f)
print(f"随机短片段，训练 3000 步：拟合区间误差 {fit:.4f}，外推区间（没有观测）误差 {ext:.4f}")

NeuralODE单摆外推

常见说法与实测

常见说法：从一个初始状态积分很长一段再和观测比，梯度要穿过很长的积分链，很难优化；常用办法是把轨迹切成短片段训练（多次打靶法）。

实测：这个 10 秒的阻尼单摆上，整条轨迹一次积分、训练 300 步时，拟合误差还卡在 0.53（另外换 5 个随机种子，有 4 个卡在 0.4–0.5）；训练到 1000 步，误差降到 0.003（5 个种子全部降到 0.001–0.004）。随机短片段训练 3000 步，用时相近，误差 0.0035，5 个种子也都成功。

怎么理解：在这个规模上，长链积分的困难表现为一段很长的平台期，训练够久就能走出来，没有到学不会的程度。轨迹更长、动力学更混沌时，长链积分的梯度会爆炸或消失，分段训练才变得必要。作为默认做法，短片段训练更稳。

典型例子