Cauchy difference priors for edge-preserving Bayesian inversion

Markku Markkanen; Lassi Roininen; Janne M. J. Huttunen; Sari Lasanen

doi:10.1515/jiip-2017-0048

开放式访问发布人：德古意特出版社 2019年1月30日

保边贝叶斯反演的柯西差分先验

马尔库·马尔卡宁，拉西·罗伊宁， Janne M.J.Huttunen女士和萨里·拉萨南

来自日志反问题与病态问题杂志

https://doi.org/10.1515/jiip-2017-0048

摘要

我们考虑未知目标包含尖锐边缘的逆问题，例如不同的材料。这类问题在图像重建、层析成像和其他反问题算法中是典型的。保边反演的一个常见解决方案是使用总变分（TV）先验值。然而，如Lassas和Siltanen 2004所示，电视节目主持人非离散化变量：边缘保持当计算网格变得越来越密集时，属性将丢失。本文提出了另一类边缘保持先验知识贝叶斯反演，柯西差分先验。我们从连续一维开始构造Cauchy先验柯西运动，并表明其离散化版本Cauchy random walk可用作保边缘贝叶斯反演的非高斯先验。我们将该方法推广到二维柯西场，并简要考虑柯西先验到Lévyα-稳定随机场先验的推广。我们开发了一种适用于单分量Metropolis–Hastings条件均值估计的后验分布抽样算法。我们将该方法应用于一维反褶积和二维X射线层析成像问题。

关键词：贝叶斯统计反问题;柯西分布;先验信息;X射线断层扫描

MSC 2010年：49号45;2015年1月62日;6020万;60J22型

1引言

我们建议使用由Cauchy行走构造的非高斯随机场先验值来处理保边贝叶斯统计逆问题[18，4]。例如，在地震和医学成像中，需要进行边缘保护反演[11，39]。我们用于保边贝叶斯反演的基本构造块是连续柯西运动及其离散化版本的柯西随机游动，其每个增量都是独立且同分布的柯西噪声。Lévyα-稳定随机游动是Cauchy的推广(α=1)和高斯(α=2)随机游走[37]。因此，我们将一维Lévyα稳定随机游动先验展开为二维随机场先验。我们将主要关注Cauchy先验函数的数值，例如，我们将对有限维和无限维后验分布之间的相互作用（即离散化方差）进行数学上严格的讨论留给未来的论文。关于无限维贝叶斯反演框架中Cauchy噪声先验的讨论，请参见Sullivan 2016[40]。侯赛尼2016年也考虑了类似的方法[17]其中研究了有界变差的Lévy过程先验。在本文中，我们关注具有无界变化的更粗糙样本路径。我们注意到，已知高斯差分先验可以提高平滑度，即在反演中平滑尖锐特征，因此它们不能用于保边反演。

Cauchy先验和保边贝叶斯反演的一个例子需要的是X射线断层扫描，也称为计算机断层扫描。在医学成像中，未知的，比如器官、牙齿和周围组织之间有接口，可以被视为锋利的边缘。贝叶斯框架内的X射线层析成像已经报道，例如[39，30，41]。处理X射线层析成像的标准方法是使用滤波反向投影[6，7，38]。在我们有足够投影数据的情况下，解决层析成像问题在数值上是稳定的[29]，我们可以重建平滑的部分以及边缘。然而，当测量次数变少、测量角度受限或测量噪声增加时，这些方法具有局限性[10，27，32]. 贝叶斯反演可以吸收不确定性和限制，为滤波反向投影提供了一种强大的选择。

保边贝叶斯反演的一个常见选择是使用总变分（TV）先验。然而，在离散化极限（Lassas）下，有限维TV先验收敛到连续参数高斯先验和Siltanen 2004[25]).这意味着，当我们使计算网格越来越密集时，贝叶斯统计逆问题的解，即后验分布，对于未知的离散化是不不变的[21，22]。在低维反演问题中，TV先验函数可能表现出保边反演特性。然而，当离散化变得更加密集时，高维TV先验函数开始失去其边缘保持特性，从而平滑了尖锐特征。基于Donsker中心极限定理[19]和稳定随机变量吸引域[12]可以看出，离散化方差的缺乏源于随机游动增量的分布。也就是说，增量具有有限的二阶矩，这通常导致分布极限。当允许增量使用其他分布时可以保持边缘保持特性高维先验。通过选择增量的柯西分布，我们可以证明离散柯西游动X（X）收敛于连续柯西运动𝒳当离散化步骤小时→0+，即它具有非高斯离散化限制。这与柯西随机变量的有限阶矩不大于1有关。

Lassas、Saksman和Siltanen 2009年[24]提出使用非高斯贝索夫空间先验，这些先验是基于小波构造的。这些可以被证明可以促进离散化-非变反演，例如，Vänskä、Lassas和Siltanen 2009将该方法应用于X射线层析成像[41]。然而，贝索夫先验通常通过小波展开定义，对于保边反演，通常使用哈尔小波基。然而，由于Haar基的结构，基函数的不连续性给出的基本并矢网格上首选不连续性。例如，在域上(0，1)，在x个=14结束x个=1三。在大多数实际情况下，这是一个强烈而不切实际的假设。Cauchy Prior没有此限制。Helin和Lassas 2011[16]，被认为是Mumford–Shah之前的等级制度。处理边缘的其他方法包括例如水平集贝叶斯反演，其中在实际中，一种方法使用高斯马尔可夫随机场先验重构[35]未知对象，然后将重建阈值设置为预定义的级别集[11，20，28]。

通过应用超参数模型，高斯先验也可以扩展为非高斯先验[三，5]。有关保边重建的层次正则化的讨论，请参见[1]。有关带Matérn域的离散不变超模型的讨论，请参见[33]，其中超验函数是用连续参数高斯和柯西随机场建模的。这些模型的基础是通过离散某些随机偏微分方程来构造逆协方差矩阵，并将长度尺度建模为超验函数中的一个随机场。这一思想类似于概率深度学习算法中的纵向建模[9，31]。结合Lasanen 2012年的结果[21，22]，这些超模型可以被证明是离散化不变的。有关贝叶斯反演中高斯马尔可夫随机场的一般参考，请参见[36，34]，在空间插值中，请参见[26]。Bolin 2014研究了具有非高斯噪声的Matérn-style随机场的构造[2]。

为了对所提方法进行数值验证，我们将该方法应用于一维反褶积和二维X射线层析成像问题。为了从具有非高斯先验的后验分布中得出估计，我们使用马尔可夫链蒙特卡罗（MCMC）方法，因为MCMC提供了条件平均估计，这是我们的主要目标。我们采用了所谓的单分量Metropolis–Hastings抽样方案，也称为Metropolis-with-in-Gibbs[11]。然而，MCMC方法因其速度慢而闻名，尤其是在具有多模态后验分布的高维问题中。对于层析成像问题，我们还讨论了高斯-牛顿型优化的最大后验估计，并将具有Cauchy先验的层析成像重建与具有高斯和全变差先验的滤波后投影估计和条件平均估计进行了比较。然而，我们想强调的是，柯西先验函数的应用可以包括任何具有尖锐特征的反问题。

本文的其余部分组织如下：在节中2首先介绍了贝叶斯统计反演框架，然后从Lévyα-稳定过程出发构造Cauchy先验。我们考虑了为什么Cauchy随机游动先验促进保边反演，最后我们将先验推广到二维。在节中三，我们回顾了单组件Metropolis–Hastings算法。第节考虑了数值一维反褶积和二维层析成像示例4。我们将所开发的层析成像算法与条件均值估计、高斯先验和TV先验以及滤波后的反向投影重建进行了比较。在节中5，我们对研究进行了总结，并对未来的研究做了一些注释。

贝叶斯反演的2个柯西先验

在贝叶斯统计反演中，我们的目标是估计未知连续参数对象的后验分布𝒳⁢(x个)，x个∈ℝd日，d日=1，2，…，给定噪声扰动的间接测量。用算符方程形式给出测量值

(2.1)米=𝒜⁢𝒳+e（电子），

哪里米∈ℝn个是已知的测量矢量，并且𝒜是从某些函数空间（如可分Banach空间）到有限维向量空间的已知线性算子。我们假设测量噪声为高斯分布e（电子）∼𝒩⁢(0，Σ)，其中Σ∈ℝn个×n个是已知协方差矩阵。我们进一步假设e（电子）在统计上独立于𝒳。

为了计算的目的，我们在方程中离散化连续观测模型(2.1)并将离散观测模型表示为

(2.2)米=A类⁢X（X）+e（电子），

哪里A类∈ℝn个×k个、和X（X）∈ℝk个。有限维贝叶斯统计逆问题的解是一个后验概率分布。通过贝叶斯公式，我们给出了后验概率密度

(2.3)D类(X（X）|米)=D类(X（X）)D类(米|X（X）)D类⁢(米)∝D类(X（X）)D类(米|X（X）)，

哪里D类(米|X（X）)是基于离散观测模型(2.2)和D类⁢(X（X）)是反映我们的先验概率密度在进行任何实际测量之前的未知信息。此外，D类⁢(米)只是一个标准化常数，因此我们可以放下它并使用非正常密度。在实践中，先验是估计算法中唯一可调的参数。

在下一小节中，我们将考虑Lévyα-稳定过程，从中可以获得Cauchy和Gaussian随机游动，如下所示特殊情况。通过证明柯西先验可以是单峰的，也可以是双峰的，后者促进了保边反演，我们研究了保边逆思想。我们还与其他常见的先验进行了比较，然后将先验推广到ℝ2。

2.1 Lévyα稳定、Cauchy和Gaussian随机游走先验

我们在这篇论文中的主要目的是做数值计算，但是在这里，我们将介绍我们感兴趣的过程的基本和更正式的概念，因为我们认为这是对读者有益的重要背景信息。但更重要的是，这为未来的研究指明了方向，在未来的研究中，将取得更多的技术成果，以证明后验分布的离散化方差。

让{𝒳⁢(吨)，吨∈[0，T型]}对某些人来说是一个随机过程0<T型<∞。此外，让S公司α⁢(τ，β，μ)是莱维参数为α，τ的α-稳定分布，β和μ，对于某些0<α≤2，-1≤β≤1。我们打电话给𝒳⁢(吨)从零开始的连续Lévyα-稳定过程，如果𝒳⁢(0)=0，𝒳具有独立增量和

𝒳⁢(吨)-𝒳⁢(秒)∼S公司α⁢((吨-秒)1α，β，0)

对于任何0≤秒<吨≤T型。有关Lévyα-稳定过程的更多详细信息，请参见[37]。为了说明参数的选择，如图1，我们绘制了α和β变化的α-稳定概率密度函数。这个案子α=2是高斯随机游动，这是一种特殊情况，也是稳定分布的高端。Cauchy walk由以下公式获得α=1。参数β对概率密度的偏度进行了建模，这在许多应用中都是一个有用的特征，但这里我们主要关注对称密度，因此β=0。

图1

具有不同α和β参数的概率密度函数。

例如，α-稳定过程可以通过使用独立分散的测度来构建。随机测量M（M）在[0，T型]映射（勒贝格）可测集A类⊂[0，T型]随机变量M（M）⁢(A类)以这样的方式M（M）⁢(⋃k个A类k个)=∑k个M（M）⁢(A类k个)几乎可以肯定的是，对于不相交的可测集A类k个⊂𝕀.随机测量M（M）称为独立分散，如果M（M）⁢(A类1)，…，M（M）⁢(A类n个)是当可测集为A类我⊂[0，T型]，我=1，…，n个，n个∈ℕ，两两不相交。我们要求

M（M）⁢(A类)∼S公司α⁢(米⁢(A类)1α，β，0)

斜度β为常数。在我们的案例中，衡量标准米⁢(⋅)在[0，T型]是勒贝格的尺度|⋅|乘以常数λ>0。那么从零开始的α-稳定Lévy运动是

𝒳⁢(吨)=M（M）⁢(χ⁢([0，吨]))，

哪里χA类是集合的特征函数A类。

对于Lévyα-稳定随机游动的连续极限，我们应用独立分散测度。让我们表示离散随机游动吨=j个⁢小时通过X（X）j个，其中j个∈ℤ+和小时>0是离散化步骤。我们获得了一个随机行走近似值（参见[37，第3.3节]）

(2.4)X（X）j个-X（X）j个-1∼S公司α⁢((λ⁢小时)1α，β，0)。

很容易看出这一点这种随机行走近似收敛于连续Lévyα稳定运动小时→0在Skorokhod空间中分布着右连续且有左极限的函数。

我们特别感兴趣的是柯西随机漫步。给定方程式(2.4)和选择α=1，β=0，我们可以写作为概率密度的柯西随机游动的联合分布

(2.5)D类⁢(X（X）)=C类⁢∏j个=1J型(λ⁢小时(λ⁢小时)2+(X（X）j个-X（X）j个-1)2)，

其中λ是正则常数C类是一个归一化常数。我们可以将这个概率密度直接应用为方程中贝叶斯公式中的先验概率密度(2.3).

在图中2，我们绘制了柯西和高斯随机游动的实现。后者是Lévyα稳定随机游动α=2。柯西随机游动实现要么有小扰动，要么有大跳跃。因此，直观地看，使用柯西分布进行差分会导致保边反演是合理的。高斯随机游走实现不喜欢跳跃，而喜欢连续路径，因此在贝叶斯反演中倾向于平滑。

图2

柯西和高斯随机游动的实现。

（a）

柯西随机游走。

（b）

高斯随机游走。

2.2保边反演

据我们所知，对于不同的先验值，还没有任何简单的标准可以表明先验值是否倾向于保边解而不是平滑解。直观的想法是构造一个先验，它促进不连续性，即跳跃。因此，通过否定，我们可以说，在保边反转中，先验函数并不能促进所有地方的平滑。

为了简单起见，我们首先考虑一维装箱取货小时=λ=1。我们将在中推广二维情况的构造接下来的部分。因此，让我们看看三个连续的坐标X（X）j个-1，X（X）j个，X（X）j个+1随机向量的X（X）。我们可以这么说X（X）更喜欢点的平滑度j个如果X（X）j个（很可能）在X（X）j个-1和X（X）j个+1并且在点上有边j个如果X（X）j个（很可能）在X（X）j个-1或在X（X）j个+1。上述想法意味着我们对以下条件先验分布感兴趣X（X）j个条件是X（X）j个-1和X（X）j个+1已知。在不失一般性的情况下，我们可以假设X（X）j个-1=-一和X（X）j个+1=一。Cauchy差分的条件密度X（X）j个然后可以写为概率密度

D类(X（X）j个|X（X）j个-1=-一，X（X）j个+1=一)∝11+(X（X）j个-一)211+(X（X）j个+一)2。

这只是两个柯西概率密度函数的乘积，每个柯西概率密度函数来自X（X）j个。为了了解这些函数的性质，考虑概率密度函数相对于X（X）j个零：

D类′′⁢(0)<0 什么时候⁢|一|<1:最大0，单峰，

D类′′⁢(0)=0 什么时候⁢|一|=1:尽可能平坦0，

D类′′⁢(0)>0 什么时候⁢|一|>1:最小值为0，双峰。

在图中三，我们绘制了这些概率密度函数。使用|一|<1我们有一个单峰函数，因此正则化可以在X（X）j个=0。使用|一|>1我们有一个双峰函数，因此正则化促进了X（X）j个=±一围绕这两个值平滑变化。这意味着这种先验会促进跳跃或小局部振荡。

图3

上部面板：Cauchy概率密度函数X（X）2，给定固定值X（X）1=-一和X（X）三=一.底部面板：相同的情况高斯差先验和总变差之前。

（a）

单峰柯西一<1。

（b）

扁平Cauchy|一|=1。

（c）

双峰柯西一>1。

（d）

高斯分布。

（e）

总变化。

我们可以为高斯差分先验和总变异先验绘制类似的图。高斯差分先验只是一个高斯形函数，其最大值为X（X）j个=0。因此，提供平滑特性。总变化先验是恒定的[-一，一]在该区域外呈指数衰减。因此，总变差既不会惩罚不连续性，也不会促进保边反演，因为它会平等地促进这两点之间的所有值。这与之前的Cauchy非常接近|一|=1即，在某种意义上，密度函数的形状类似于总变分先验密度函数。二者对重建的影响α→0和α=2-¦Β可以用类似的论据来理解。对于α→0，模态变为尖峰，而对于α=2-¦Β是双峰的，对于非常小的ϵ，它变成了三峰的，最后是单峰的。如果α=2，这可视为特殊情况。

2.3二维差异先验

让我们考虑一个二维晶格(小时⁢j个，小时'⁢j个')，其中(j个，j个')∈𝕀2⊂ℤ2以及每个坐标方向的离散化步骤小时，小时'>0。高斯差分先验可以通过具有i.i.d.增量的差分方程来构造

(2.6)X（X）j个，j个'-X（X）j个-1，j个'∼𝒩⁢(0，σ2⁢小时小时')，X（X）j个，j个'-X（X）j个，j个'-1∼𝒩⁢(0，σ2⁢小时'小时)。

术语σ2>0是正则化参数。然后通过条件正态分布的乘积构造先验

D类公共关系⁢(X（X）)⑪经验⁡(-(X（X）T型⁢我1T型⁢C类1-1⁢我1⁢X（X）+X（X）T型⁢我2T型⁢C类2-1⁢我2⁢X（X）))，

哪里我1和我2是两个不同坐标方向的差分矩阵，以及C类1=σ2⁢小时小时'⁢我和C类2=σ2⁢小时'小时⁢我，其中我是一个单位矩阵。这些结构相当有名，参见示例[23]。我们强调这个等式(2.6)不可显式求解，因为只存在逆协方差矩阵。对于最小二乘解的唯一性，我们需要施加边界条件。

让我们用一个简单的例子来考虑一下，为什么我们按小时小时'和小时'小时。考虑方程中的第一个差分方程就足够了(2.6)，因为另一个是通过相同的参数获得的。我们分两步考虑网格加密：（1）我们考虑更密集的离散化j个-方向，我们称之为正则化坐标方向，和（2）更密集的离散化j个'-方向。让X（X）j个，j个'-X（X）j个-1，j个'∼𝒩⁢(0，σ2)。如果我们沿着正则化坐标方向加密晶格j个，这将对应于一维高斯随机游走，因此我们需要在方差中添加离散化步长，因此X（X）j个，j个'-X（X）j个-1，j个'∼𝒩⁢(0，σ2⁢小时)。

第二个坐标方向有点复杂。让我们考虑一个简化的情况，如图所示4。因此，我们定义像素平均值

X（X）1=X（X）11+X（X）122 和 X（X）2=X（X）21+X（X）222。

缩放2的原因是因为这是标准的像素级算术平均值，即我们想要像素的自然缩放。现在增量

X（X）2-X（X）1=(X（X）21+X（X）22)-(X（X）11+X（X）12)2=(X（X）21-X（X）11)+(X（X）22-X（X）12)2。

如果我们选择X（X）2-X（X）1∼𝒩⁢(0，σ2)，那么我们必须选择X（X）21-X（X）11∼X（X）22-X（X）12∼𝒩⁢(0，2⁢σ2)。因此，根据以下比例1小时'在方程式中(2.6).

作为替代方案，使用随机偏微分方程可能很有吸引力Δ12⁢𝒳=𝒲，其中Δ是拉普拉斯算子。然而，在概率分布的意义上，你可以考虑一个方程组，写为

(2.7)∂1⁡𝒳∼𝒲1，∂2⁡𝒳∼𝒲2

其逆协方差算子在最小二乘意义下为拉普拉斯Δ。我们喜欢在方程式中使用公式(2.7)而不是Δ12⁢𝒳=𝒲，使用分数阶拉普拉斯算子Δ12在数值上相当繁琐。

图4

网格加密到第二个坐标方向。

同样，对于柯西差异，我们可以使用i.i.d.增量

(2.8)X（X）j个，j个'-X（X）j个-1，j个'∼柯西⁢(λ⁢小时)，X（X）j个，j个'-X（X）j个，j个'-1∼柯西⁢(λ⁢小时')。

结果是通过与高斯情况类似的论证得出的。如果我们使网格沿正则化方向更加密集，这同样对应于一维Cauchy随机行走，即我们只需使用离散化步长进行缩放小时。

在第二个坐标方向上，我们再次进行与图中相同的分割4。如果X（X）2-X（X）1∼柯西⁢(λ)，然后X（X）21-X（X）11∼X（X）22-X（X）12∼柯西⁢(λ)（注意Cauchy分布随机变量的独立线性组合）。

通过类似的论证，我们实际上可以看到，一般Lévyα-稳定差分先验的标度是以i.i.d.增量表示的

(2.9)X（X）j个，j个'-X（X）j个-1，j个'∼S公司α⁢(小时'⁣1-αα⁢小时1α，β，0)，X（X）j个，j个'-X（X）j个，j个'-1∼S公司α⁢(小时1-αα⁢小时'⁣1α，β，0)。

这自然意味着在更高维度上有类似的结构。

备注2.1。

方程的可解性(2.9)是保证，当(2.9)得到适当的调节。如果我们在定义方程之前用我⁢X（X）=Z轴，其中Z轴k个是独立的稳定随机变量，那么Z轴被限定在以下范围内我.可解方程为我X（X）=Z轴|对(我)，其中条件随机变量Z轴|对(我)∈对(我)概率为1。让我们将正交投影表示为对⁢(我)具有问.条件概率P（P）(Z轴∈A类|对(我))=∫A类c（c）D类Z轴(z（z）)d日z（z）用于Borel集合A类⊂对⁢(我)可以用通常的限制程序导出（在坐标变换后）

P（P）(Z轴∈A类|对(我))=林k个→∞P（P）(Z轴∈A类×B类k个)P（P）⁢(对⁢(我)×B类k个)

通过考虑类型集A类×B类k个在对⁢(我)×对⁢(我⊥)，其中递减超立方体B类k个正测度的具有交集0。我们之前的密度就是那个时候c（c）⁢D类Z轴⁢(问⁢我⁢x个)=c（c）⁢D类Z轴⁢(我⁢x个)，其中附加的标称常数c（c）提醒我们所执行的调节。附加的先验赋范常数在后验分布的表达式中抵消。先验密度D类公共关系⁢(x个)通过将所有概率密度函数

D类S公司α⁢(小时'(α-1)/α)⁢小时1/α)⁢(X（X）j个，j个'-X（X）j个-1，j个')⁢D类S公司α⁢(小时(α-1)/α)⁢小时'⁣1/α)⁢(X（X）j个，j个'-X（X）j个，j个'-1)

为所有人j个，j个'。

图5

二维Cauchy、Gaussian和TV先验的实现显示出块状和平滑结构。

（a）

二维柯西先验的实现。

（b）

二维高斯先验的实现。

（c）

二维电视的实现。

我们注意到，先验密度不是指数有界的。因此，无论何时，只要直接理论A类在里面(2.2)具有非平凡的null空间。

在图中5，我们绘制了二维高斯先验和柯西先验的实现，定义在方程中(2.6)和(2.8)分别是。我们还绘制了总变化先验的实现图。我们使用了零边界条件，但为了演示实现的块状或平滑结构，我们对图像进行了裁剪。柯西先验实现的一个显著特征是，它看起来块状，并且有明显的锐边。为了绘制实现，我们使用了单组件Metropolis–Hastings，这将在下一节中进行解释，即我们生成了Cauchy先验的实现链。

我们注意到，柯西先验不是旋转不变的，这可以从实现中清楚地看出。这可以解释为一个堆积问题，即为了具有旋转不变的柯西先验，边缘应该是球形的。但我们不能有效地包装球，因此实现了箱形。在当前构造中，盒形仅限于坐标方向。这可以通过进行坐标变换来改变，因此可以使用定向nabla算子，但这会给我们带来更复杂的公式。使用有向导数并让方向在超验函数中建模确实是一个有趣的想法。类似地，TV previor在其标准公式中也不是旋转不变的，但它可以通过重新公式的旋转不变来制作，如[15]. 从某种意义上说，各向异性TV先验与Cauchy先验一样具有坐标方面的边缘特征。

如果我们考虑二维Cauchy先验的保边特征，我们可以研究沿任一坐标方向的过程。这些一维过程看起来实际上类似于一维柯西随机游动。为了进行比较，我们还绘制了高斯差分先验和以类似方式获得的TV先验的实现。正如我们所见，高斯分布促进了平滑度，而不是跳跃。电视先验的实现不像柯西先验那样块状，但更像高斯先验。

3个单一组成的大都市——黑斯廷斯

为了获得贝叶斯反问题的实用解，我们需要根据后验分布计算点估计。最常用的点估计是最大后验（MAP）估计或条件平均（CM）估计。对于MAP估计，我们使用高斯-牛顿型优化方法。然而，我们的兴趣主要在于CM估计，即后验密度的平均值（见等式(2.3))。在高维问题中，CM的显式计算估计值（后验密度的积分）很少典型后向密度的可能数值积分导致不可行的方法。因此，CM估计值通常使用马尔可夫链蒙特卡罗（MCMC）计算采样方法，如Metropolis–Hastings（MH）算法或Gibbs取样器（参见示例[14，13]).MCMC方法的思想是生成样本{X（X）(秒)}从后验分布和近似CM估计使用系综平均值。同样的算法也可以用于从先前的分布中提取实现，如第2节所述。

在本文中，我们将所提出的Cauchy先验应用于一维反褶积问题X射线层析成像问题。这两个问题都有一个特定的属性，在评估后验密度，参数中单个分量的更新向量的计算成本非常低。这是因为前向理论矩阵是稀疏的，并且类似地，先验是沿分量变化的，即更新一个分量时，我们不需要计算完全后向，而只需要计算被改变的部分。利用该属性的常见算法有单组分吉布斯采样器和单组分大都会-黑斯廷斯（MH）算法；参见示例[14]。

这些算法的思想是从组件级采样条件分布。我们从一些初始值开始

X（X）(0)=(X（X）1(0)，X（X）2(0)，…，X（X）n个(0))

并设置秒=0。首先，我们考虑一维条件的分布X（X）1给定数据米以及所有其他参数：

(3.1)D类(X（X）1|米，X（X）2(0)，X（X）三(0)，…，X（X）n个(0))=C类1D类(X（X）1，X（X）2(0)，X（X）三(0)，…，X（X）n个(0)|米)，

哪里C类1是一个归一化常数。我们生成了一个样本X（X）1(1)来自该分布。然后继续到下一个参数并绘制示例X（X）2(1)来自分布

D类(X（X）2|米，X（X）1(1)，X（X）三(0)，…，X（X）n个(0))=C类2D类(X（X）1(1)，X（X）2，X（X）三(0)，…，X（X）n个(0)|米)，

哪里C类2是一个归一化常数。该过程一直持续到所有参数一次更新一个，我们已经X（X）(秒+1)。我们增加了秒←秒+1并重复样本生成N个秒获取一组样本的次数{X（X）(秒)}秒=1N个秒。

吉布斯采样器和单组分MH之间的差异取决于一维分布样本的方式绘制。在吉布斯采样器中，使用来自分配。然而，吉布斯采样器的缺点是一维条件分布很少以这种形式出现允许快速生成样本，并且是数值的（计算上低效）必须使用样本生成。参见示例[14]了解更多详细信息。

在单组分MH中，吉布斯采样器中的直接采样为替换为Metropolis–Hastings型采样。例如，要绘制样本X（X）1(1)从分布中(3.1)，我们先画一个候选人样品X（X）~1来自提案分发问1(⋅|X（X）1(0))，这是可以接受的

α1=最小值⁡(1，D类(X（X）~1|米，X（X）2(0)，…，X（X）n个(0))问1(X（X）1(0)|X（X）~1)D类(X（X）1(0)|米，X（X）2(0)，…，X（X）n个(0))问1(X（X）~1|X（X）1(0)))。

如果接受，我们设置X（X）1(1)=X（X）~1。否则我们设置X（X）1(1)=X（X）1(0)。对其他部件重复此步骤。

原则上，提案分配的选择是任意的，但这种选择对算法的收敛性和选择不当的提案分布可能需要很大的样本足够收敛的大小；参见示例[14，11]。我们注意到，通过使用基于无网格采样算法的建议分布，可以减轻这种影响[8]，则不同网格上的收敛性相似。本文中，建议分布被选择为高斯分布：问我(⋅|X（X）我(0))=𝒩(X（X）我(0)，σ2)。差异σ2选择的样品合格率为25%-50%。通常在开始时忽略一些样品（老化周期），因为收敛到平稳状态需要时间分布D类(X（X）|米)。本文选择老化期为样本大小。

4数值示例

为了证明具有Cauchy差分先验的Bayesian反演的保边和离散方差特性，我们考虑了一维反褶积问题和二维X射线层析成像的例子。

4.1一维反褶积

在图中6，我们有一个综合反褶积例子，在不同网格上具有Cauchy和总变分先验。原始未知量是一个分段常数函数，我们观察到卷积噪声测量。常数卷积核是余弦形的。在测量中，我们添加了带有标准偏差的白噪声σ=0.3，对应于∼10%相对比率。

图6

离散化水平对CM-估计的影响柯西和总变异先验：顶部面板：噪声去卷积观测值（实线）。中间面板：Cauchy之前对不同格子（实线），底部面板–与中间面板相同，但之前有总的变化。所有图形中的虚线对应于真实目标，灰色区域对应于2⁢σ可信度区间。

（a）

测量。

（b）

柯西N个=131。

（c）

柯西N个=261。

（d）

柯西N个=521。

（e）

总变化量N个=131。

（f）

总变化量N个=261。

（g）

总变化量N个=521。

我们已经用一些未知因素进行了重建N个=131，261，521。方程中的Cauchy正则化参数λ(2.5)选择的方式使重建看起来边缘保留。与使用任何高斯先验值类似，在选择λ时需要谨慎，以避免使用过强或过松的先验值-这需要进行一些调整，在未来的研究中，该参数应建模为类似于例如[11]，其中高斯先验参数在超先验中建模。在更稠密的网格上，我们简单地通过离散化步长缩放来改变正则化参数小时相应地。在不同网格上重建的Cauchy先验图明显提供了边缘保护功能。我们还绘制了2⁢σ可信度区间，这些看起来也很相似N个除了跳跃，在跳跃中，边不会收敛，而是增长。我们怀疑这是因为链条的长度太短。在测度弱收敛的意义上，后验是离散不变的。随机变量后验的一个简单例子X（X）2鉴于Y（Y）1=X（X）j个-1+ε和Y（Y）2=X（X）j个+1+ε，何时Y（Y）1=Y（Y）2=X（X）0=0，表明

D类邮递⁢(z（z）)=c（c）⁢经验⁡(-12⁢x个2-12⁢v（v）2)⁢1(1+x个2)⁢(1+(x个-z（z）)2)⁢(1+(z（z）-v（v）)2)。

因此X（X）j个是有界的，即使我们只有相邻点的观测值。同样，我们缩放TV先验正则化参数。我们注意到估计量的离散化方差行为不能可以从当我们增加未知数。关于前面示例中总变化的进一步讨论，请参见Lassas和Siltanen 2004[25]。

在单组件Metropolis–Hastings中，我们使用了1000000条长链以确保收敛。我们在单分量算法中需要长链，因为我们处理的是多峰分布，众所周知，单分量算法会陷入一种模式，要达到正确的模式需要很长时间。对于总变差和柯西情况，老化期均选择为链长的一半。MCMC链可以优化，但在本文中，我们的目标是简单地证明柯西先验，因此我们将MCMC优化研究推迟到后续论文中。

4.2二维X射线层析成像

现在我们将把所开发的方法应用于X射线层析成像。我们将谢普-洛根幻影视为最初的未知。测量值对应于二维扇形梁几何结构（见图7用于描述测量的几何形状）。源和探测器到旋转轴的距离为4和2个。检测器的宽度为3，检测器元素包含200像素。我们在41个不同的角度（从-10∘到190∘具有5∘步骤）。此外，我们添加了大约0.1%噪声级的白色测量噪声。

图7

合成二维层析成像示例的测量几何。

我们使用单组分Metropolis–Hastings计算CM房地产，如第节所述三。链是从反向投影估计开始的(A类T型⁢米适当缩放）。为了实现MCMC链的收敛，我们使用500万作为链长，其中250万个样品被忽略老化期。为了节省内存，我们只存储十分之一用于计算集合平均值的样本。这里，链长度是指单个组件更新的数量，因此在2D模拟中具有最高分辨率800×800×5×106=3.2×1012样品。然而，我们不需要评估全部可能性，因为我们可以简单地使用前向理论矩阵的稀疏性和差异先验的iid属性，这允许实际的计算时间。重建的计算时间200×200分辨率约为24小时，对于400×400分辨率，时间约为96小时。这表明算法的计算成本与预期的未知数（几乎）成正比。计算是通过带有双Xeon的PC工作站进行的E5-2680 CPU和128 GB内存。该算法的核心部分是使用MEX/C接口实现的。

为了进行比较，我们还使用计算成本较低的Cauchy先验来计算MAP估计。MAP估计是通过最小化损失函数来计算的

我μ(X（X）)=-日志D类(X（X）|米)-μ∑我日志(X（X）我+¦Β)，

其中添加额外的惩罚项以确保解的非负性，从而提高优化的稳定性。通过使用高斯-纽顿优化进行最小化，从而在迭代过程中减小μ和ϵ。我们注意到，作为非凸优化问题，高斯-牛顿迭代有收敛到（次优）局部极小值的趋势，这在这里也是如此（见下文）。

图8

晶格尺寸不同估计值的比较400×400像素。

（a）

未知：Shepp–Logan幻影。

（b）

过滤的反向投影估计。

（c）

Cauchy先验MAP估计。

（d）

柯西先验CM估计。

（e）

先用TV进行CM估算。

（f）

基于高斯先验的CM估计。

在图中8，我们用五种不同的方法绘制了原始未知和重建图。过滤后的反向投影估计恢复了大多数特征，但存在严重的伪影。具有Cauchy先验的CM估计很好地捕捉到了特征，而没有显著的伪影。带有TV先验的CM估计也很好地捕捉了特征，但不如Cauchy先验估计的尖锐。此外，当我们减少测量次数或增加噪声时，柯西估计似乎更稳健。高斯先验估计恢复了主结构，但它既有不良伪影，也有平滑结构。

图8还显示了使用Cauchy先验计算的MAP估计。该估计也很好地捕捉了特征，而CM估计在重建中具有较少的伪影。然而，如上所述，高斯-纽顿迭代有收敛到（次优）局部极小值的趋势，这里的情况似乎也是如此：从CM估计开始迭代可以得到成本更低的MAP估计我0（即较高的后验概率），且该估计值与CM估计值在视觉上不可区分。但当然，从构型管理开始并不是一种实际意义重大的方法。然而，寻找非凸成本函数的全局（最优）极小值通常被认为是一项困难的任务，因为我们的主要动机是MCMC估计，所以寻找最优MAP估计有待于未来的研究。

图9

三种不同格点上X射线层析成像问题的条件平均估计。

$（a） 200×200{200\乘以200}。$

（a）

200×200。

$（b） 400×400{400\乘以400}。$

（b）

400×400。

$（c） 800×800{800\乘以800}。$

（c）

800×800。

图10

计算域中间条件平均值估计值的横截面。

（a）

水平。

（b）

水平放大。

（c）

垂直。

（d）

垂直放大。

以数字表示9和10，我们对不同尺寸的格子进行了CM估计。如本节前面所述，已考虑离散化，因此重建看起来基本相同。

5结论与讨论

全变分先验作为非高斯先验在保边贝叶斯反演中是传统的常用选择。然而，如本文所示[25]，全变分先验不促进离散化-变分反演。此前，有人建议使用Besov-space优先解决此问题。Cauchy先验和更一般的Lévyα稳定先验可以为离散化保边贝叶斯反演提供一类非高斯先验。

柯西先验的吸引力来自这样一个事实，即我们可以将先验与柯西步行联系起来，柯西步行具有明确且公认的统计特性。我们在独立同分布增量的帮助下构建了Cauchy先验，这为我们通过条件分布理解保边性质提供了一种直观的方法。

为了进行数值验证，我们将开发的Cauchy先验函数应用于一维反褶积问题与二维X射线层析成像，其中我们计算了条件平均值估计使用单组件Metropolis–Hastings算法（最大a后验估计值也包括在X射线断层扫描中）。与具有总变差先验的贝叶斯反演相比提出的方法促进了边缘保持贝叶斯反演和当我们测量次数较少，测量噪声增加或当计算网格非常密集时。

在未来的研究中，我们应该考虑具有数学严密性的离散化方差。此外，还应该对Lévyα-稳定随机游动进行详尽的数值研究。而高斯随机游动α=2显然是某种特殊情况和稳定分布的高端，Cauchy随机变量α=1处于值的中间0<α≤2。对于0<α<1我们可能会在跳跃之间有比之前的柯西更稳定的行为。同样适用于案例1<α<2在柯西（Cauchy）和高斯（Gaussian）之间，我们的行为可能不太稳定，但平滑推理显然只适用于特殊情况α=2。此外，我们可以通过形式的随机偏微分来尝试更一般的先验模型𝒫⁢𝒳=𝒲，其中𝒫是一个线性微分算子，𝒳未知且𝒲低沉的噪音。这类似于Matérn场，即高斯马尔可夫随机场的构造，只是我们会有一般的Lévy噪声，而不是特殊情况α=2高斯噪声。此外，用Student的t分布代替Cauchy噪声将是一个有趣的发展，因为Cauchys是Student t的特例。MCMC采样方案的进一步开发也可以通过使用并行化和阻塞的Metropolis-with-in-Gibbs算法来完成。为了加快MCMC的收敛速度，还应该研究多模式提案分发。这些也将有助于增强我们的多峰分布单分量采样器，因为有时我们会陷入一种模式，无法摆脱它。

资金来源：工程和物理科学研究委员会

奖励标识/授予编号：EP/K034154/1号

资金来源：Luonontieteiden ja Tekniikan Tutkimuksen Toimikunta公司

奖励标识/授予编号：250215

资金报表：这项工作得到了英国工程和物理科学研究委员会（EPSRC参考：EP/K034154/1–大尺度逆问题的不确定性量化）和芬兰科学院（申请号250215、307741和313709）的资助。

工具书类

[1]J.M.Bardsley、D.Calvetti和E.Somersalo，PET图像保边重建的层次正则化，反问题26（2010），第3期，文章ID 035010。10.1088/0266-5611/26/3/035010在谷歌学者中搜索

[2]D.博林，非高斯噪声驱动的空间Matérn场，扫描。《J Stat.41》（2014），第3期，557–579。10.1111/约12046在谷歌学者中搜索

[3]D.Calvetti和E.Somersalo，用于恢复块状对象的高斯超模型，反问题23（2007），第2期，733–754。10.1088/0266-5611/23/2/016在谷歌学者中搜索

[4]D.Calvetti和E.Somersalo，贝叶斯科学计算导论，Surv公司。导师。申请。数学。科学。2,施普林格，纽约，2007年。在谷歌学者中搜索

[5]D.Calvetti和E.Somersalo，贝叶斯成像框架中的超模型，反问题24（2008），第3期，文章编号034013。10.1088/0266-5611/24/3/034013在谷歌学者中搜索

[6]A.M.Cormack，用线积分表示函数，有一些放射学应用。我，J.应用。物理学。34（1963），文章ID 2722。10.1063/1.1729798在谷歌学者中搜索

[7]A.M.Cormack，用线积分表示函数，有一些放射学应用。二、，J.应用。物理学。35 (1964), 195–207.10.1063/1.1713127在谷歌学者中搜索

[8]S.L.Cotter、G.O.Roberts、A.M.Stuart和D.White，函数的MCMC方法：修改旧算法使其更快，统计师。科学。28（2013），第3期，424–446。10.1214/13-STS421在谷歌学者中搜索

[9]A.C.Damianou和N.D.Lawrence，深高斯过程，第十六届国际人工智能与统计会议记录，PMLR，斯科茨代尔（2013），207-215。在谷歌学者中搜索

[10]M.E.Davison，有限角度层析成像问题的病态性质，SIAM J.应用。数学。43（1983），编号：2428-448。10.1137/0143028在谷歌学者中搜索

[11]M.M.Dunlop、M.A.Iglesias和A.M.Stuart，层次贝叶斯水平集反演，统计计算。27（2017），编号61555–1584。2007年10月17日/11222-016-9704-8在谷歌学者中搜索

[12]W.Feller，概率论及其应用导论。第二卷，约翰·威利父子公司，纽约，1966年。在谷歌学者中搜索

[13]A.Gelman、J.B.Carlin、H.S.Stern和D.B.Rubin，贝叶斯数据分析，第二版。，文本统计。科学。序列号。，查普曼和霍尔/CRC，博卡拉顿，2004年。2019年10月10日/9780429258480在谷歌学者中搜索

[14]W.R.Gilks、S.Richardson和D.J.Spiegelhalter，马尔可夫链蒙特卡罗实践，Interdiscip公司。统计人员。，查普曼和霍尔出版社，伦敦，1996年。10.1201/b14835在谷歌学者中搜索

[15]G.González、V.Kolehmainen和A.Seppänen，阻抗层析成像中的各向同性和各向异性总变差正则化，计算。数学。申请。74 (2017), 564–576.2016年10月10日/j.camwa.2017.05.004在谷歌学者中搜索

[16]T.Helin和M.Lassas，统计逆问题和Mumford–Shah泛函中的层次模型，反问题27（2011），第1期，文章编号015008。10.1088/0266-5611/27/1/015008在谷歌学者中搜索

[17]B.侯赛尼，具有无限可分和重尾先验测度的适定贝叶斯反问题，SIAM/ASA J.不确定性。数量。5（2017），第1期，1024–1060。10.1137/16M1096372在谷歌学者中搜索

[18]J.Kaipio和E.Somersalo，统计和计算反问题，申请。数学。科学。160,斯普林格，纽约，2005年。2007年10月1日/138659在谷歌学者中搜索

[19]O.卡伦伯格，现代概率基础，第2版。，普罗巴伯。申请。（纽约），施普林格，纽约，2002年。10.1007/978-1-4757-4015-8在谷歌学者中搜索

[20]E.Klann、R.Ramlau和W.Ring，用于SPECT/CT数据反演和分割的Mumford–Shah水平集方法，反向探测。《成像5》（2011），第1期，137-166。10.3934/ipi.2011.5.137在谷歌学者中搜索

[21]S.Lasanen，非高斯统计反问题。第一部分：后验分布，反向探测。《成像6》（2012年），第2期，215–266页。10.3934/ipi.2012.6.215在谷歌学者中搜索

[22]S.Lasanen，非高斯统计反问题。第二部分：近似未知数的后验收敛，反向探测。《成像6》（2012年），第2期，267–287页。10.3934/ipi.2012.6.267在谷歌学者中搜索

[23]S.Lasanen和L.Roininen，格林先验的统计反演，第五届工程反问题国际会议：理论与实践，Taylor&Francis，伦敦（2005），第3–32页。在谷歌学者中搜索

[24]M.Lassas、E.Saksman和S.Siltanen，离散变贝叶斯反演和贝索夫空间先验，反向探测。《成像3》（2009），第1期，第87–122页。10.3934/ipi.2009.3.87在谷歌学者中搜索

[25]M.Lassas和S.Siltanen，在保边贝叶斯反演中可以使用总变分先验吗？，反问题20（2004），第5期，1537-1563。10.1088/0266-5611/20/5/013在谷歌学者中搜索

[26]F.Lindgren、H.v.Rue和J.Lindström，高斯场和高斯-马尔可夫随机场之间的明确联系：随机偏微分方程方法，J.R.统计社会服务。B统计方法。73（2011），第4期，423–498。10.1111/j.1467-9868.2011.00777.x在谷歌学者中搜索

[27]A.K.Louis，x射线计算机断层扫描中的数据不完整问题。I.有限角度变换的奇异值分解，数字。数学。48（1986），第3期，251-262。2007年10月10日/BF01389474在谷歌学者中搜索

[28]D.芒福德和B.吉达斯，通用图像的随机模型，夸脱。申请。数学。59（2001），第1期，85–111。10.1090/qam/1811096在谷歌学者中搜索

[29]F.纳特勒，计算机断层成像的数学，B.G.Teubner，斯图加特，1986年。10.1007/978-3-663-01409-6在谷歌学者中搜索

[30]E.Niemi、M.Lassas、A.Kallonen、L.Harhanen、K.HäMäLäinen和S.Siltanen，使用时空水平集方法的动态多源X射线层析成像，J.计算。物理学。291 (2015), 218–237.2016年10月10日/j.jcp.2015.03.016在谷歌学者中搜索

[31]C.J.Paciorek和M.J.Schervish，使用一类新的非平稳协方差函数进行空间建模，Environmetrics 17（2006），第5期，483–506。10.1002/env.785在谷歌学者中搜索公共医学公共医学中心

[32]E.T.Quinto，X射线变换和有限数据层析成像的奇异性𝐑2和𝐑三，SIAM J.数学。分析。24（1993），第5期，1215–1225。10.1137/0524069在谷歌学者中搜索

[33]L.Roininen、M.Girolami、S.Lasanen和M.Markkanen，Matérn场的超验函数及其在贝叶斯反演中的应用，反向探测。《成像》第13期（2019年），第1期，第1-29页。10.3934/ipi.2019001在谷歌学者中搜索

[34]L.Roininen、J.M.J.Huttunen和S.Lasanen，贝叶斯统计反演的Whittle-Matérn先验及其在电阻抗断层成像中的应用，反向探测。《成像》第8期（2014年），第2期，561-586页。10.3934/ipi.2014.8.561在谷歌学者中搜索

[35]L.Roininen、M.S.Lehtinen、S.Lasanen、M.Orispää和M.Markkanen，相关性先验，反向探测。《成像5》（2011年），第1期，167–184页。10.3934/ipi.2011.5.167在谷歌学者中搜索

[36]L.Roininen、P.Piiroinen和M.Lehtinen，构造连续平稳协变量作为二阶随机差分方程的极限，反向探测。《成像》第7期（2013年），第2期，611-647页。10.3934/ipi.2013.7.611在谷歌学者中搜索

[37]G.Samorodnitsky和M.S.Taqqu，稳定的非高斯随机过程，斯托克。型号。，查普曼和霍尔出版社，纽约，1994年。在谷歌学者中搜索

[38]L.A.Shepp和J.B.Kruskal，计算机断层扫描：新的医用X射线技术，阿默尔。数学。《月刊》第85期（1978年），第6期，第420–439页。10.1080/00029890.1978.11994611在谷歌学者中搜索

[39]S.Siltanen、V.Kolehmainen、S.Järvenpää、J.P.Kaipio、P.Koistinen、M.Lassas、J.Pirttilä和E.Somersalo，少数射线医学x射线断层摄影的统计反演：I.一般理论，物理学。医学生物学。48 (2003), 1437–1463.10.1088/0031-9155/48/10/314在谷歌学者中搜索公共医学

[40]T·J·沙利文，适定贝叶斯反问题和重尾稳定拟巴拿赫空间先验，反向探测。《成像11》（2017），第5期，857–874。10.3934/ipi.2017040在谷歌学者中搜索

[41]S.Vänskä、M.Lassas和S.Siltanen，使用经验Besov先验的统计X射线层析成像，国际层析成像杂志。《统计》第11卷（2009年），编号S09，3–32。在谷歌学者中搜索

收到：2017-05-05

修订过的：2018-05-19

认可的：2018-12-25

在线发布：2019-01-30

印刷出版：2019-04-01

本作品根据Creative Commons Attribution-NonCommercial-NoDerivatives 3.0许可证授权。