《Geophysical Journal International》：Transport-Map Proposals for Efficient Markov Chain Monte Carlo
编辑推荐：
本文针对高维复杂后验分布中MCMC采样效率低下的难题，提出了一种基于可逆神经网络构建传输映射（Transport Map）的新型Metropolis-Hastings提案框架。研究人员通过将目标分布映射到简单参考分布，在参考空间设计提案并利用逆映射回目标空间，结合改进的接受准则，在保留精确采样的同时将自相关时间降低2.5-6倍。该方法为计算成本高昂的地球物理反问题提供了更高效的贝叶斯不确定性量化工具。
在地球科学领域，揭示地球内部结构严重依赖于地表收集的间接观测数据。这些数据往往不完整且含有噪声，使得通过求解不适定反问题来关联地下属性变得极具挑战。传统上，地球物理学家通过寻找满足特定最优性准则的单一“最佳拟合”模型来解决此类问题。近几十年来，贝叶斯方法日益普及，因为它能以观测数据为条件，刻画模型的后验概率分布，从而内在地量化不确定性。
在此背景下，马尔可夫链蒙特卡洛（MCMC）方法已成为关键工具，它能够从任意复杂的后验分布中进行采样。许多MCMC方法的核心是Metropolis-Hastings算法，该算法结合提案分布和概率接受准则来构建马尔可夫链，使其稳态分布为所需的目标分布。该算法通用性强，对提案分布的要求宽松，因此在科学领域广泛应用。然而，当参数空间庞大或正演模型计算成本高昂时（这两种情况在地球物理学中都很常见），该算法可能变得不切实际。高维性通常导致链混合不佳（即参数空间探索效率低下）以及向目标分布收敛缓慢，而昂贵的似然评估则限制了通常可承受的马尔可夫链步数。
这些挑战催生了多个研究方向。一种思路是彻底放弃MCMC，承认精确采样可能成本过高。这种理念支撑了变分推理方法，该方法用优化代替采样，并寻求目标分布的近似。第二种方法保留采样，但通过用替代模型（通常是神经网络）替代正演模型来加速似然评估，从而降低成本。这产生了近似但显著更快的MCMC运行。本研究追求的第三条研究路线是保留精确的MCMC和真实的似然，转而专注于提案机制，以实现参数空间的高效探索。
法布里齐奥·马格里尼和马尔科姆·桑布里奇在《Geophysical Journal International》上发表的研究，引入了一种基于自适应传输映射来设计Metropolis-Hastings提案的新方法。该框架受到计算统计学和最优传输最新进展的启发，其核心思想是找到一个可逆映射，将（复杂的）目标分布重铸为更简单的参考分布，后者更适合标准的MCMC步骤。与先前要求映射分量单调且具有三角雅可比矩阵的方法不同，本研究放宽了这些约束，转而采用基于可逆神经网络的参数化方法，这些网络通过设计保证可逆性，同时提供更大的灵活性和表达能力。
该方法迭代进行。首先使用标准采样策略（例如随机游走Metropolis）从目标分布获取初始样本。这些样本用于训练一个可逆神经网络，该网络学习将目标分布映射到一个更简单的参考分布（本研究中使用标准多元高斯分布）。随后的提案在参考空间中生成，可以是全局独立移动或局部扰动，并根据适当修改的Metropolis-Hastings准则进行接受或拒绝。随着更多样本的积累，更新网络参数，从而提高整体采样效率。最终，该方法不仅产生代表所需目标分布的样本集合（与标准MCMC一样），还通过网络的权重获得了目标分布的紧凑学习表示。
为开展研究，作者主要应用了以下几个关键技术方法：首先，利用可逆神经网络（特别是仿射耦合层构成的网络）参数化传输映射，确保映射的可逆性和雅可比行列式严格为正；其次，基于从目标分布获得的初步MCMC样本，通过最小化Kullback-Leibler（KL）散度的样本近似来训练传输映射；接着，在参考空间（标准高斯分布）设计局部高斯提案，并通过逆映射将提案转换回目标空间；最后，采用改进的Metropolis-Hastings接受准则，确保采样过程精确收敛到目标后验分布。研究涉及的数值实验包括二维Rosenbrock分布、三维地震定位问题以及高达16维的高斯混合模型。
2 数学背景
研究首先提供了支撑所提出方法的基本数学背景，该方法寻找可逆变换以将目标分布重新参数化为参考分布。该方法的思想源于更广泛的最优传输理论，但不需要解决完整的Monge问题，后者需要归一化的概率密度。这使得该方法非常适用于地球物理中的贝叶斯反问题，因为后验分布通常只能通过未归一化的密度评估或样本来获得。
2.1 传输映射的构建
目标是构建一个传输映射T，使得如果θ服从分布π，则T(θ)服从分布ρ。如果这一要求被精确满足，则映射诱导的密度?π将等于目标密度π。一个自然的选择是最小化两者之间的Kullback-Leibler（KL）散度。通过推导，优化目标简化为最小化关于T的期望项。
2.2 优化问题
在实践中，由于目标分布π通常只能通过样本或未归一化密度评估间接已知，优化目标中的期望一般无法以闭合形式计算。为了解决这个困难，作者采用样本平均近似来替代期望。具体地，假设已从目标分布获得K个样本，则可以用样本平均来近似期望。如果选择参考分布ρ为标准多元高斯分布，则可以进一步简化目标函数。最终，将映射参数化为由向量γ参数化的映射，得到具体的优化问题。
2.3 备注
从优化问题中推导出的传输映射?T至少包含两种不同的近似。第一种是前面讨论的样本平均近似。第二种是用有限维子集Tγ替代无限维函数空间T。与变分推理不同，变分推理是从变分族中选择一个分布来近似后验本身，而本框架中的近似映射?T用于构建MCMC采样的提案分布。因此，映射的准确性在本框架中并非严格关键。只要?T在整个域上保持可逆，并且使用它提出的样本通过适当的Metropolis-Hastings准则被接受或拒绝，基于传输映射提案的MCMC采样仍然是精确的。
3 作为可逆神经网络的传输映射
先前的研究采用了基于多元多项式或径向基函数的参数化，这些参数化因为生成的映射在展开系数γ上是线性的，从而导致优化问题是凸的而具有吸引力。然而，三角约束固有地限制了映射的表达能力，并带来了困难的模型选择问题。此外，这些展开方式无法直接评估逆映射?T-1，而逆映射对于在每个马尔可夫链步骤中将样本从参考空间变换到目标空间至关重要。
受这些限制的启发，本研究采用了不同的方法。以面临更富挑战性的（非凸）优化问题为代价，使用由仿射耦合层组成的可逆神经网络来参数化传输映射。仿射耦合层定义了一种可逆变换，其雅可比矩阵是三角矩阵且对角线元素严格为正。通过组合多个这样的层（一些具有下三角雅可比矩阵，另一些具有上三角雅可比矩阵），可以获得表达力强的可逆变换，其整体雅可比矩阵不再是三角矩阵，但仍然具有严格正的行列式。
3.1 仿射耦合层
仿射耦合层通过将输入向量分割为两个子向量，并保持其中一个子向量不变，同时使用依赖于第一个子向量的元素仿射函数来变换另一个子向量。该变换的可逆性很容易验证，其雅可比矩阵具有分块下三角形式，且行列式严格为正。
3.2 互补仿射耦合层
为了增强表达力，引入了“互补”仿射耦合层，其中两个子向量的角色被互换。该层保持第二个子向量不变，同时变换第一个子向量。其雅可比矩阵具有分块上三角形式，行列式同样严格为正。
3.3 仿射耦合块
通过组合标准仿射耦合层和互补仿射耦合层，形成仿射耦合块。该块的雅可比矩阵通常是稠密的，但行列式保持严格为正。通过组合多个仿射耦合块，可以构建表达力更强的传输映射。
3.4 传输映射参数化
在本研究中，每个缩放或平移函数都由具有单个隐藏层（含h个神经元）的前馈神经网络参数化。每个仿射耦合块包含四个神经网络。对于由M个块组成的传输映射，参数向量的维度随维度n线性增长，而不是像基于多项式或径向基展开的方法那样指数增长。在数值实验中，发现M=3和h=16能在广泛的问题范围内产生具有足够表达力的传输映射。
3.5 训练目标与停止准则
通过最小化一个正则化的成本函数来获得参数化映射?T的最佳权重和偏置。该成本函数包含一个正则化项，用于惩罚与恒等映射的偏差。使用基于自适应矩估计的随机优化器进行迭代最小化。训练前，将传输映射初始化为恒等映射，并对目标样本进行标准化。训练最多进行1000次迭代（周期）。为了监控收敛，考虑KL散度中被积量的方差，将其作为训练期间的停止准则。
4 使用传输映射提案的MCMC
Metropolis-Hastings方案迭代运行：从当前状态θ开始，从提案分布q(θ′|θ)中抽取新状态θ′。然后以一定概率接受提议状态。在温和的技术假设下，重复应用此过程保证产生的样本的经验分布收敛到目标分布。
在获得传输映射?T后，可以通过首先从参考分布ρ中抽取r′，然后应用逆变换得到θ′ = ?T-1(r′)来提议新状态。考虑在参考空间中定义的提案分布qr(r′|r)。通过变量变换公式，目标空间中的相应提案密度可知。代入Metropolis-Hastings准则，得到修改后的接受概率。在本研究中，在参考空间中使用局部提案，即具有对角协方差矩阵和恒定标准差σr=0.5的多元高斯分布。由于该提案是对称的，接受概率仅取决于当前和提议状态的目标密度和雅可比行列式。
该方法将提案设计与目标分布的具体特征解耦，减少了手动调整的需要。研究发现，固定的提案标准差σr=0.5能在不同问题中产生稳健的性能。
4.1 从Rosenbrock分布采样
作为第一个例子，考虑一个二维问题，其中目标分布π是Rosenbrock分布。采样策略迭代进行。首先使用标准Metropolis-Hastings算法从目标分布中抽取一组样本。然后使用初始样本集合训练一个可逆神经网络，表示传输映射?T1。从此时起，提案在参考空间中构建，每个提议状态根据修改后的准则进行接受或拒绝。使用基于?T1的提案收集更多样本后，从组合样本集合中推导出第二个传输映射?T2。这个更新的映射取代?T1用于采样，并用于生成更多样本。类似地，训练第三个映射?T3。随着采样的进行，连续的传输映射产生越来越接近真实目标分布的拉回密度。在此示例中，训练最终映射?T4会产生特别低的诊断值，相应的拉回密度?π4与π几乎无法区分。这表明在合适的情况下，拉回密度可以作为目标分布的替代，从而可以丢弃最终的样本集合而不会造成显著信息损失。
随着采样过程的进行，MCMC接受率稳步提高，当首次引入传输映射提案时，接受率从约25%显著跃升至65%。这一趋势反映了随着更多样本的获得和传输映射变得更加准确，参考空间提案被转换为越来越反映目标分布形状的目标空间提案。
为了系统评估改进效果，比较了两组各20个独立马尔可夫链的性能：一组使用上述基于传输映射的采样范式，另一组使用传统的Metropolis-Hastings算法。对于每组链，收集2500个样本。结果表明，使用传输映射提案的链的积分自相关时间（IACT）中值明显低于标准MCMC（16.1 vs. 40.6）。这意味着标准MCMC需要大约2.5倍的样本才能达到相同的有效样本大小。然而，在这种情况下，训练连续传输映射所产生的开销相对较大，导致配备参考空间提案的链的运行时间明显长于标准Metropolis-Hastings。
4.2 地震定位
作为第二个例子，考虑了熟悉的地球物理问题：从地震台阵记录的震相到达时间推断震源位置和发震时间。为了简化，将地震表示为具有横向速度变化的二维介质中的点源，并使用射线理论近似波传播。采样问题是重建后验分布。将每个源参数视为独立且均匀分布的随机变量。使用高斯似然模型。
类似于Rosenbrock例子，使用两种MCMC采样器对后验进行采样：一种是在目标（后验）空间中使用以当前状态为中心的高斯提案的标准Metropolis-Hastings采样器；另一种是使用参考空间提案的传输映射驱动采样器。对每种采样器运行20个独立的马尔可夫链，每个链初始化在真实的源参数处以便于公平比较采样器效率。对于标准采样器，采用特定的标准差。基于传输映射的采样器依赖于在采样过程中逐步更新的六个映射序列。
该实验的结果表明，传输映射采样器的采样效率更高，IACT中值约为标准采样器的四分之一（9.1 vs. 35.3）。这表明标准Metropolis-Hastings需要大约四倍的样本才能达到后验估计的可比精度。这一结论得到了单个马尔可夫链运行结果的支持：在传输映射的情况下，直方图紧密跟随后验边缘分布，表明收集的样本可以支持后验统计量的准确推断。相比之下，标准采样器的直方图显示出欠采样分布特有的锯齿状轮廓，某些区域的概率质量被显著低估，而其他区域则被高估。
关于训练开销与额外采样的讨论，在此问题中，单次后验评估的成本（平均0.1秒）远大于Rosenbrock例子。训练六个传输映射耗时约80秒，这仅占表中马尔可夫链平均运行时间的约10%。由于正演问题的相对较高成本，这种开销远小于标准MCMC达到可比有效样本大小所需的额外时间。
4.3 高维采样
为了证明传输映射驱动的MCMC适用于高维问题，考虑了n=4、8和16维的目标分布。这些目标分布由高斯混合模型给出。对于每个维度，生成20个这样的分布，并使用标准采样器和传输映射采样器对每个进行采样。后者涉及在大小线性增加的样本批次上训练的六个传输映射。对于标准采样器，使用根据随机生成的目标凭经验选择的局部提案（从标准多元高斯分布中抽取）。对每个采样器运行单个马尔可夫链。
与前面章节一致，这些实验揭示了两种采样器效率的明显差异。基于传输映射的MCMC在所有维度上都实现了显著更低的中值IACT——在4维中降低约3倍，在8维和16维中降低超过6倍。对于标准采样器，IACT随维度的急剧增加反映了随机游走提案在高维中的性能恶化：探索效率越来越低，因为提案未能与目标分布的各向异性或弯曲几何结构对齐。相比之下，基于传输映射的提案内在地适应目标的几何结构。当映射具有足够的表达能力来捕捉目标的主要特征时，即使维度增加，提案也能保持良好的尺度和方向。这解释了传输映射采样器中值IACT增长缓慢的原因，从4维的5.9到16维的8.1。
需要注意的是，与Rosenbrock例子类似，此处每个马尔可夫链步骤评估目标密度的成本可以忽略不计。因此，训练多个传输映射所需的时间大大超过了标准MCMC采样器。
4.4 收敛性考虑
上述示例中阐述的采样策略将方法置于自适应MCMC算法类别中，其中提案机制随着采样的进行而演变。对于任何定义提案的固定传输映射，由修改后的接受概率给出的转移规则或核保持固定，从而产生理论上保证从真实目标分布π产生样本的马尔可夫链。然而，每当发生自适应（即更新传输映射）时，核会相应改变，产生一个时间非齐次的马尔可夫链。在这种状态下，标准的收敛结果不能直接应用。
在自适应MCMC理论中，两个条件对于正式证明收敛至关重要： diminishing adaptation（意味着链核的连续变化变得越来越小）和containment（确保自适应核序列保持良好行为，不会导致诸如链逃逸到无穷远等病态行为）。作者提到，他们的工作与Parno和Marzouk（2018）的研究密切相关，后者通过验证这两个条件为他们基于映射加速的自适应采样器建立了收敛性。然而，作者的框架利用了由可逆神经网络参数化的传输映射，这带来了重要的技术优势，但通过最小化非凸成本函数进行训练。因此，Parno和Marzouk（2018）获得的理论保证不能直接转移到作者的设置中，正式的收敛性分析仍然是一项不平凡的任务。
在实现中，随着采样的进行，传输映射的更新变得越来越不频繁，成本函数中的正则化项以及将映射初始化为恒等映射原则上应阻止更新之间的突然变化。与此观点一致，作者观察到，一旦有适量样本可用，连续映射诱导的拉回密度是逐渐变化的，而不是突变的，后期的更新仅导致微小的调整。虽然这些考虑
