文献翻译__人工智能时代医学图像重建中的凸优化算法(前3章)

文章下载–我的Gitee

Convex optimization algorithms in medical image reconstruction—in the age of AI

人工智能时代医学图像重建中的凸优化算法

徐景颜一和弗雷德里克·诺2

2022年3月23日发布• 2022医学物理与工程研究所
医学和生物学中的物理学, 第67卷,7号引用徐景颜和弗雷德里克·诺2022物理医学。生物。 6707TR01

https://iopscience.iop.org/article/10.1088/1361-6560/ac3842#)

摘要

在过去的十年里，基于模型的图像重建(MBIR)算法得到了迅速的发展，它们往往是最优化领域的凸优化算法的应用或改编。我们回顾了在医学图像重建中广泛流行的一些最新算法，强调了不同算法之间的已知联系，并讨论了诸如计算和存储开销等实际问题。最近，深度学习进入了医学成像领域，最新的发展试图利用深度学习和MBIR之间的协同作用来提升MBIR的性能。我们介绍了现有的方法和动态增强的MBIR方法的新趋势，特别关注了凸性和凸性算法在网络体系结构中的潜在作用。我们还讨论了如何利用凸性来提高DL网络的泛化能力和表示能力。

1.介绍

在过去的十年里，基于模型的图像重建(MBIR)方法在CT、MR、PET和SPECT中得到了广泛的研究。许多出版物记录了这些MBIR方法的好处，从减少图像伪影和总体上提高图像质量，到减少CT应用中的辐射剂量。MBIR问题通常被描述为优化问题，其中由数据拟合项和正则化项组成的标量目标函数关于未知图像被最小化。在如此大规模和数据密集型应用的推动下，同一时期数学界也对开发凸优化算法进行了密集的研究。凸优化概念的引入引发了许多新的研究方向，如具有快速收敛特性的MBIR算法，以及更好地捕捉先验图像信息的新颖的正则化设计。

最近，**深度学习方法(DL)**在许多复杂的现实世界任务中取得了超乎人类的表现。它们在解决医学成像问题方面的快速采用和适应也取得了有效成果。关于反问题的动态链接法的出版物数量呈爆炸式增长。作为这种快速发展的证据，一些特刊(Greenspan等人2016年、Wang等人2018年、Duncan等人2019年)和评论文章(McCann等人2017年、Lucaset al 2018年、Willemink和Noël 2019年、Lell和Kcherie 2020年)已经出版，以总结当前的最先进技术。

许多文章总体上讨论了人工智能和深度学习的优势和挑战，还有一些文章讨论了它们在医学成像中的作用和未来。一个值得警惕的观点是，DL的力量应该得到承认，但它并不是解决所有问题的灵丹妙药。似是而非的是，DL可以与传统方法协同工作，例如凸优化：传统方法的优势可能是DL的不足之处。例如，数字语言经常因为可解释性低而受到批评。另一方面，凸优化以其丰富的结构而闻名，它可以用来编码结构信息，并在与DL网络结合时提高可解释性。DL也需要数据(Marcus 2018)；它需要大量具有已知基本事实的培训数据来进行培训或评估。DL可用于增强传统MBIR方法的性能，进而为DL训练产生高质量的地面真实标签。

以此为背景，本文回顾了凸优化的基本概念，讨论了在MBIR问题中得到广泛应用的流行的一阶算法，并用文献中的实例应用来展示凸性在人工智能和数字图书馆时代的相关性。以下是本文主要内容的提纲。

部分2:凸优化中的元素

部分3:凸优化的确定性一阶算法

部分4:凸优化的随机一阶算法

部分5:非凸优化中的凸性

部分6:凸性、图像重建和深度学习的协同集成

部分7:结论

部分8:附录-其他主题，如Bregman距离，泊松似然的相对光滑性，以及一些计算实例。

2.凸优化中的元素

我们首先介绍在整篇论文中使用的通用符号。只与特定部分相关的符号将在本地引入。然后，我们解释凸分析的基本概念和结果，这有助于理解论文的内容，尤其是章节3,4，以及5。

2.1.注释

我们用 $I_C$ 表示集合C的示性函数，即如果 $x\in C$，则$I_C(x)=0$，否则，$+\infty$。一个集合 $C\subset R^n$ 是凸集当且仅当(当且仅当)对所有$x_1，x_2\in C$，$\alpha x_1+(1-\alpha )x_2\in C$。函数$f：R^n\to R\cup +\infty$的定义域定义为Dom $f=x|f(x)<\infty$；如果函数$f$的定义域不为空，则$f$是真的。一个函数$f$是闭的，如果它的上图epi $ f={(x,t) \in R^{n+1} | f(x)<=t,x \in dom f }$ 是闭的。一个函数f是下半连续的，如果它的上图是闭的(Bauschke等人，2011)，引理1.24。函数 $f：C\subset R^n\to R\cup +\infty$ 是凸的，如果 $C\subset R^n$ 是凸集，且对于$α \in $ [0，1]和$x_1，x_2\in C$，$f(\alpha x_1+(1-\alpha )x_2)<=\alpha f(x_1)+(1-\alpha )f(x_2)$。我们用缩写CCP来表示一个函数 $f$ 是凸的、闭的和真的。为方便起见，我们可以简单地把这类函数称为凸函数。
我们用$\langle .,.\rangle$表示 $a，b\in R^n$ 的两个向量的内积，即$\langle a,b\rangle ={\sum }_i{a}_i{b}_i$。内积诱导范数表示为 $\Vert \cdot {\Vert }_2$ 或简称为 $\Vert \cdot \Vert$，即 $∥x∥= \sqrt{〈a,b〉} $ 。如果没有特别说明，本文使用的范数是2-范数。

2.2.基本定义和属性

一阶算法根据它们所设计的目标函数的类型进行分类。在不同的类型中，光滑的目标函数是最常见的假设，也可能是最容易处理的。让 $Q\subseteq R^d$ 。如果一个凸函数f是可微的，并且它的梯度 $\nabla f$是Lipschitz连续的，即存在一个常数 $L>0$ 使得：
$$\Vert ▽f(x)-▽f(y)\Vert <=L\Vert x-y\Vert ,\forall x,y\in Q\text{\hspace{1em}(2.1)}$$
则 $f$ 在Q上是L-光滑的。从Nesterov等人2018年的定理2.1.5开始，这类函数可以等价地刻画为：
$$0<=f(y)-[f(x)+<▽f(x),y-x>]<=\frac{L}{2}\Vert y-x{\Vert }^2，x,y\in Q\text{\hspace{1em}(2.2)}$$
该关系表明L-光滑函数允许任何 $x，y\in Q$ 的二次函数。（2）中的常数L是梯度Lipschitz常数。

函数$f：Q\subseteq R^d\to R\cup +\infty$是 $\sigma -$强凸的，则
$$f(\alpha x_1+(1-\alpha )x_2)<=\alpha f(x_1)+(1-\alpha )f(x_2)-\frac{1}{2}\alpha (1-\alpha )\sigma \Vert x_1-x_2{\Vert }^2\text{\hspace{1em}(2.3)}$$
对于$ α\in$[0，1]，对所有$x_1，x_2 \in Q$，当函数f可微时，$σ$强凸函数可由下式给出：
$$f(x)+<▽f(x),y-x>+\frac{1}{2}\sigma \Vert x_1-x_2{\Vert }^2<=f(y)\text{\hspace{1em}(2.4)}$$
让$f:Q\subseteq R^d\to R\cup +\infty$是CCP，和$x\in Q$，f 在 x的次微分$\partial f(x)$，定义为：

∂ f ( x ) : = { u ∈ Q ∣ f ( y ) > = f ( x ) + < u , x − y > , f o r − a l l − y ∈ Q } (2.5) ∂f(x):=\{u \in Q|f(y) >=f(x)+<u,x-y>,for-all- y \in Q\}\tag{2.5} ∂f(x):={u∈Q∣f(y)>=f(x)+<u,x−y>,for−all−y∈Q}(2.5)
集合 $\partial f(x)$ 的元素称为x处的次梯度。真凸f的次微分 $\partial f(x)$ ，对于$x\in ri$ dom $f$ 是非空的(Bauschke等人，2011年，第228页)。CCP $f$ 的极小值可以用费马法则来刻画，该规则指出 $x$ 是 $f $ 的极小值当且仅当 $0\in \partial f(x)$ (Rockafella and Wets，2009年，第422页)。

共轭函数 $f^{\ast }$ 关于 $f$ ：$R^d\to R\cup +\infty$ 定义为
$$f^{\ast }(t)=su{p}_s<s,t>-f(s)\text{\hspace{1em}(2.6)}$$
当 f${}^{\ast }$ 可视为由 s 参数化的 t 的线性函数的逐点上确界，(2.6)中的 $f^{\ast }$ 总是所有 f 的凸函数。 $f^{\ast }$ 的共轭函数定义了双共轭：
$$f^{\ast \ast }(p)=su{p}_t<p.t>-f^{\ast }(t)\text{\hspace{1em}()}$$
此外，无论f如何，$f^{\ast \ast }$都是凸的。此外，还可以证明，如果f是CCP，则$f^{\ast \ast }=f$(Bauschke等人，2011，第13章)；其他的$f^{\ast \ast }<=f$，且对于任何凸函数$g<=f$，则$g<=f^{\ast \ast }$。也就是说，双共轭$f^{\ast \ast }$是f的最紧的凸下界，也称为凸包络f。下面的对偶关系将f的次微分和它的共轭$f^{\ast \ast }$联系在一起(Rockafella和Wets 2009，命题11.3)。对于任何CCP f，有 $\partial f^{\ast }=(\partial f{)}^{-1}$，并且$\partial f=(\partial f^{\ast }{)}^{-1}$；更具体地，
$$v\in \partial f(x)<==>x\in \partial f^{\ast }(v)<==>f(x)+f^{\ast }(v)=<v,x>\text{\hspace{1em}()}$$
总的来说，$f(x)+f^{\ast }(v)=<v,x>$，对所有的x,v。从上面来看，
a r g max ⁡ v { < v , x > − f ∗ ( v ) } = ( ∂ f ∗ ) − 1 ( x ) = ∂ f ( x ) (2.7) arg\max_v \{<v,x>-f^*(v)\}=(∂f^*)_{-1}(x)=∂f(x)\tag{2.7} argvmax​{<v,x>−f∗(v)}=(∂f∗)−1​(x)=∂f(x)(2.7)
类似地，
a r g max ⁡ x { < v , x > − f ( x ) } = ∂ f ∗ ( v ) (2.8) arg\max_x \{<v,x>-f(x)\}=∂f^*(v)\tag{2.8} argxmax​{<v,x>−f(x)}=∂f∗(v)(2.8)
作为基本的例子，当$f(x)=\Vert x\Vert ,x\in R$时，则$f^{\ast }(y)=I_{[-1,1]}$；二次函数$\frac{1}{2}\Vert .{\Vert }^2$是自共轭的。其他凸共轭对可以在(Bauschke等人2011年，第13章)、(Boyd等人2004年，第3章)和(Beck 2017，附录B)中找到。

如果f是CCP且μ-强凸的，则其共轭f是1/μ-光滑的(Bauschke等人，2011年，命题14.2)。反之，如果f是CCP且L-光滑的，则其共轭f是1/L-强凸的。出于这个原因，有时L-光滑CCP函数也被称为L-强光滑(Ryu和Boyd 2016)。

对于CCP $f:R^d\to R\cup +\infty$和参数$\mu >0$。近似映射和Moreau包络(或莫罗-约西达正则化)有下式定义
p r o x u f ( x ) : = a r g min ⁡ y { f ( y ) + 1 2 u ∥ x − y ∥ 2 } e u f ( x ) : = min ⁡ y { f ( y ) + 1 2 u ∥ x − y ∥ 2 } (2.9) prox_{uf}(x):=arg\min_y\{f(y)+\frac1{2u}\|x-y\|^2\}\\\tag{2.9} e_uf(x):=\min_y\{f(y)+\frac1{2u}\|x-y\|^2\}\\ proxuf​(x):=argymin​{f(y)+2u1​∥x−y∥2}eu​f(x):=ymin​{f(y)+2u1​∥x−y∥2}(2.9)
当f(y)是凸的，目标函数(2.9)或者(2.10)是强凸的，因此近邻映射$pro{x}_{\mu f}$总是单值的。当$f(y)=I_C(y)$时，则$pro{x}_f:=x_{\ast }$是与$x_{\ast }$最接近的点，使得$x_{\ast }\in C$，即投影运算。在这个意义上，近似映射(11)是凸集上投影的推广，其中f不限于指示函数。简单函数的近邻映射计算的例子可以在(Combettes和Pesquite 2011，Parikh和Boyd 2014，Beck 2017)中找到，无论是使用闭合形式的解还是使用有效的数值算法。接下来，某些函数可能被称为简单，这是以相同的方式解释的，即它们的近邻映射很容易计算或以封闭的形式存在。

如果f是CCP，则Moreau包络(2.10)是1/μ-光滑的；它的梯度∇eμf，由下式给出：
$$▽e_{u}f(x)=\frac{1}{u}(x-y_{\ast }),where{y}_{\ast }=pro{x}_{uf}(x)\text{\hspace{1em}(2.11)}$$
是1/μ利普希茨连续的(包施克等人，2011年)。从这个角度来看，Moreau包络(12.121)提供了一种从下面用光滑函数逼近潜在非光滑函数f的一般方法。更确切地说，在(Rockafella And Wets2009)中，定理1.25中$e_{\mu }f<\infty$,并且$e_{\mu }f(x)$是μ和x的连续函数，使得对所有x，$e_{\mu }f(x)\nearrow f(x)$，如$\mu \searrow 0$。${}^3$众所周知的f和$e_{\mu }f$对是：$f(y)=I(y)$，$e_{\mu }f$是障碍函数的二次型；(2)$f(y)=|y|，y\in R$，$e_{\mu }f(x)$是Huber函数。

3这句话对非凸函数f也是有效的，只要asf是从下面有界的。然而，对于非凸函数，不能保证μf是光滑的。

莫罗恒等式描述了函数f的近邻映射与其共轭函数f*之间的关系：
$$x=pro{x}_{\tau f}(x)+\tau pro{x}_{f^{\ast }/\tau }(\frac{x}{\tau })\text{\hspace{1em}(2.12)}$$
继续近似映射是广义投影概念的类比，当专门用于正交投影时，莫罗恒等式(13)可被解释为向量通过其在线性子空间L上的投影及其正交补$L^{\perp }=y|\langle y，x\rangle =0，\forall x\in L$的分解(Parikh And Boyd2014)。

可以通过将(2.9)中的二次距离替换为Bregman距离来推广最近映射(2.9)。设h是一个可微的强凸函数${}^4$，考虑由h参数表示的以下‘距离’：
$$D_h(y;x)=h(y)-[h(x)+<▽h(x),y-x>]\text{\hspace{1em}(2.13)}$$
这是由Bregman(Bregman 1967)首先研究的，14年后由 Censor 和Lent( Censor 和Lentt，1981)跟进，随后进行了更多的工作( Censor 和 Zenios 1992，Bauschke和Borwein 1997)。${}^5$ h的凸性意味着对任何x，y，$D_h(y;x)>=0$；h的强凸性意味着当y=x时$D_h$达到唯一极小值.当$h=\frac{1}{2}\Vert .{\Vert }_2^2$时，则定义(2.13)得到$D_h(y;x)=\frac{1}{2}\Vert y-x{\Vert }_2^2$。从这个意义上说，$D_h(y;x)$确实是二次距离函数的推广。作为另一个例子，如果h是加权平方2-范数，即$h=\frac{1}{2}\Vert .{\Vert }_2^2$，其中$M>0$是正定对称矩阵，则$D_h(y;x)=\frac{1}{2}\Vert y-x{\Vert }_M^2$。一般来说，与距离函数不同，$D_h(y;x)$在y和x之间不是对称的；换句话说，$D_h(y;x)!=D_h(x;y)$是可能的。

通过在(2.13)中插入Bregman距离(2.9)来定义Bregman近端映射，即，
y + ( x ) = a r g min ⁡ y { f ( y ) + 1 u D h ( y ; x ) } ( ) y_+(x)=arg \min_y\{f(y)+\frac1{u}D_h(y;x)\}\tag{} y+​(x)=argymin​{f(y)+u1​Dh​(y;x)}()
通过选择适合问题几何意义的h函数，Bregman距离$D_h(y;x)$可以用来简化计算。例如，当f(y)是$R^d$中的单位单纯形时，即$f(y)=I_C(y)$，其中$C=y|{\sum }_i{y}_i=1，y_i\in [0，1]$，近邻映射(在单纯形上的投影)没有闭合形式的解；但是选择$h(x)={\sum }_i{x}_{i}log{x}_i$，Bregman近邻映射可以用闭合形式计算(Tseng 2008)。为方便起见，我们可以简单地用$D_h(y;x)$来表示Bregman距离，而不显式地指定h函数。

Moreau包络(2.10)是两个CCP函数小卷积的特殊情况，定义为:
f ( x ) = inf ⁡ y { f 1 ( y ) + f 2 ( x − y ) } (2.14) f(x)=\inf_y\{f_1(y)+f_2(x-y)\}\tag{2.14} f(x)=yinf​{f1​(y)+f2​(x−y)}(2.14)
由于映射$(x,y)\to f_1(y)+f_2(x-y)$在x和y上是联合凸的，并且部分最小化保持了凸性，所以f的下卷积是一个凸函数。如果f1和f2都是CCP，此外，如果f1是强制性的，f2从下有界，则(16)中得到下界，可以用min代替(Bauschke et al 2011，命题12.14)。对于CT应用，已使用下值卷积(16)来组合具有互补性质的正则化(Chambolle and Lions 1997, Bredieset al 2010, Xu and Noo 2020)。粗略地说，(16)中的‘inf’运算可以‘算出’f1和f2之间的哪个分量导致较低的成本f(x)，因此更适合本地图像内容。

4这里，强凸性如(3)中所定义，但关于一般范数，不一定是由内积诱导的2-范数。有关更多详细信息，请参阅附录A.1。

5感兴趣的读者可以在(Facchinei和Pang 2003，第1232页)中找到简短的书目评论。

3.凸优化的确定性一阶算法

我们介首先介绍了一阶算法及其加速算法，然后讨论了它们在反问题求解中的应用。就内容而言，这一部分与同一主题的几篇评论论文(Cevheret al 2014，Komodakis and Pesquite 2015)、书籍或专著(Bubeck 2015，Chambolle and Pock 2016，Beck 2017)有部分重叠。感兴趣的读者应该查阅这些出版物，以获得我们没有涵盖的材料。我们的讨论集中在各种算法之间的相互关系，以及相关的内存和计算问题，当它们应用于典型的图像重建问题时。另一个目的是准备第6节，其中来自凸优化和深度学习的元素交织在一起，以利用它们之间的协同作用。

3.1.凸优化中的一阶算法

在优化领域已经开发了许多一阶算法。这些算法只利用了函数值及其梯度的信息，即使对于像图像重建这样的大规模问题也很容易计算。不同算法之间的差异通常在于它们对问题模型/结构的假设。

本节包含三个小节。在前两节中，我们讨论了原始-对偶混合梯度(PDHG)算法和(预处理)ADMM算法。这两种算法在成像应用中受到了极大的欢迎。在最后一节中，我们讨论了最小化三个函数之和的最新进展，其中一个是与线性算子复合的非光滑函数；相关的三块算法可以比传统的两块类型的前两个算法更有效地存储。

3.1.1.非光滑凸优化的原始对偶算法

考虑以下优化模型:
$$\underset{x}{\min }g(x)+h(Kx)\text{\hspace{1em}(3.1)}$$
其中g，h都是CCP，$x\in R^d$和$K：R^d\to R^q$是一个线性算子，其范数为$\Vert K\Vert$，已知。由于通常很难按原样处理复合形式$h(K\cdot )$，原始对偶算法将目标函数(17)重新表示为最小-最大凸凹问题。我们首先使用$h(\cdot )$的(双)共轭函数重写$h(\cdot )$:
$$h(Kx)=\underset{z}{\max }<Kx,z>-h^{\ast }(z)\text{\hspace{1em}(3.2)}$$
然后得到(3.1)的原始-对偶改写为:
$$\underset{x}{\min }\underset{z}{\max }g(x)+<Kx,z>-h^{\ast }(z)\text{\hspace{1em}(3.3)}$$
双重目标函数由下式给出${}^6$
max ⁡ z { − sup ⁡ x { < − K t z , x > − g ( x ) } − h ∗ ( z ) } = max ⁡ z − g ∗ ( − K t z ) + h ∗ ( z ) (3.4) \max_z\{-\sup_x\{<-K^tz,x>-g(x)\}-h^*(z)\}=\max_z-{g^*(-K^tz)+h^*(z)}\tag{3.4} zmax​{−xsup​{<−Ktz,x>−g(x)}−h∗(z)}=zmax​−g∗(−Ktz)+h∗(z)(3.4)
原始-对偶混合梯度(PDHG)算法在原始下降步和对偶上升步之间交替。一个简单的变种(Chambolle And Pock 2011)如下:
z k + 1 : = a r g min ⁡ z { < K x ~ k , z > − h ∗ ( z ) − 1 2 σ ∣ ∣ z − z k ∣ ∣ 2 } = p r o x σ h ∗ ( z k + K x ~ k ) x k + 1 : = a r g min ⁡ x { g ( x ) + < K x , z k + 1 > + 1 2 τ ∣ ∣ x − x k ∣ ∣ 2 } = p r o x τ g ( x k − τ K ∗ z k + 1 ) x ~ k + 1 = x k + 1 + θ ( x k + 1 − x k ) (3.5) z_{k+1}:=arg\min_z \{<K\widetilde{x}_k,z>-h^*(z)-\frac1{2\sigma}||z-z_k||^2\}=prox_{\sigma h^*}(z_k+K\widetilde{x}_k)\\\tag{3.5} x_{k+1}:=arg\min_x \{g(x)+<Kx,z_{k+1}>+\frac1{2\tau}||x-x_k||^2\}=prox_{\tau g}(x_k-\tau K^*z_{k+1})\\ \widetilde{x}_{k+1}=x_{k+1}+\theta(x_{k+1}-x_k) zk+1​:=argzmin​{<Kx k​,z>−h∗(z)−2σ1​∣∣z−zk​∣∣2}=proxσh∗​(zk​+Kx k​)xk+1​:=argxmin​{g(x)+<Kx,zk+1​>+2τ1​∣∣x−xk​∣∣2}=proxτg​(xk​−τK∗zk+1​)x k+1​=xk+1​+θ(xk+1​−xk​)(3.5)
当θ=1，且(21)中的步长$\sigma ，\tau$满足$\tau \sigma \Vert K{\Vert }^2<1$时，在(Chambolle And Pock2011)中显示了该算法以O(1/k)的遍历率${}^7$以部分原始-对偶间隙的形式收敛。

6我们用p和d分别表示(17)和(18)中的原始目标值和对偶目标值。一般而言，弱对偶性成立，即$p>=d$。两者的等价性(强对偶)可以在g，h和线性映射K上的较温和的条件下建立，这是ʼ对偶定理的推广。有关更多细节，请参见(Rockafella 2015，第31节)。

7这些速率是根据迭代的加权平均值来衡量的，而不是迭代本身。对于(21)，证明了O(1/k)对于$x^k=({\sum }_i^k{x}_i)/k$，其中$x_i$来自(21b)。

3.1.2.非光滑凸优化的ADMM

ADMM考虑了以下约束问题:
$$\begin{gathered} \underset{x,z}{\min }\varphi (x,z)=g(x)+h(z) \\ s.t.v=Kx+Bz\text{\hspace{1em}(3.6)} \end{gathered}$$
式中$g：R^d\to R\cup +\infty 、h：R^m\to R\cup +\infty$均为CCP。问题数据由$(v，K，B)$，$K：R^d\to R^q$和$B：R^q\to R^p$组成，是线性映射，而$v\in R^q$是给定的向量。目标函数$\varphi$在满足(22)中耦合约束的未知数$x，z，x\in R^d，z\in R^p$中是可分的。我们引入了约束的拉格朗日乘子$\in R^q$，形成了增广拉格朗日函数:

$$\begin{gathered} L(x,z;\lambda )=g(x)+h(z)+<\lambda ,Kv+Bz-v>+\frac{u}{2}||Kx+Bz-v|{|}^2 \\ =g(x)+h(z)++\frac{u}{2}||Kx+Bz-v+\frac{\lambda }{u}|{|}^2-\frac{1}{2u}||\lambda |{|}^2\text{\hspace{1em}(3.7)} \end{gathered}$$
其中，μ>0是恒定步长参数。ADMM算法的基本版本以交替的方式与以下更新公式一起更新(23)中的原始变量x、z和拉格朗日乘数λ:
z k + 1 ∈ a r g min ⁡ z { h ( z ) + u 2 ∣ ∣ K x k + B z − v + λ k u ∣ ∣ 2 } x k + 1 ∈ a r g min ⁡ x { g ( x ) + u 2 ∣ ∣ K x k + B z k + 1 − v + λ k u ∣ ∣ 2 } λ k + 1 : = λ + u ( K x k + 1 + B z k + 1 − v ) (3.8) z_{k+1}\in arg\min_z\{h(z)+ \frac {u}{2}||Kx_k+Bz-v+\frac {\lambda_k}{u}||^2 \}\\\tag{3.8} x_{k+1}\in arg\min_x\{g(x)+\frac {u}{2}||Kx_k+Bz_{k+1}-v+\frac {\lambda_k}{u}||^2\}\\ \lambda_{k+1}:=\lambda+u(Kx_{k+1}+Bz_{k+1}-v) zk+1​∈argzmin​{h(z)+2u​∣∣Kxk​+Bz−v+uλk​​∣∣2}xk+1​∈argxmin​{g(x)+2u​∣∣Kxk​+Bzk+1​−v+uλk​​∣∣2}λk+1​:=λ+u(Kxk+1​+Bzk+1​−v)(3.8)
对偶序列{${\lambda }_k$}和原目标{$g(x_k)+h(z_k)$}的收敛可以建立当子问题(24a、24b)都有解时，即迭代继续。保证子问题解存在的较弱条件和反例可以在(Chen等人2017)中找到。

在应用中的一种常见情况是两个线性映射K，B中的一个是简单的。${}^8$假设B很简单，即$B=I$或$B^{t}B=I$，则(24a)中的更新允许以$pro{x}_{h/\mu }(\cdot )$的形式的解。
如果没有进一步的假设K，更新$x_{k+1}$可能不接受直接解决方案。已经提出了带有预条件或线性化的ADMM的变体，以使子问题(24b)更容易。算法3.1是ADMM(Beck 2017)的这样一种变体，在x更新时具有预条件矩阵M。

算法3.1 问题(22)的预条件ADMM算法。
Input: Choose $x_0,{\lambda }_0$, let $u>0$.
Output: $x_k,z_k,{\lambda }_k$
1 for $iter=0,...,K-1$ do
2 $z_{k+1}:= arg\min_z{h(z)+ \frac {u}{2}
3 $x_{k+1}\in arg\min_x{g(x)+\frac {u}{2}
4 ${\lambda }_{k+1}:=\lambda +u(K{x}_{k+1}+B{z}_{k+1}-v)$ /* dual ascent */

如果选择M为
$$M=\frac{1}{{\tau }^{\prime }}I-u{K}^{t}K\text{\hspace{1em}(3.9)}$$
如果$0<{\tau }^{\prime }<u||K|{|}^2$，则M是正定矩阵；算法3.1的$x_{k+1}$更新中的极小值问题允许以邻近算子$pro{x}_{g{\tau }^{\prime }}(.)$的形式存在唯一解，从而简化了问题.。在(Beck 2017)中可以找到推广版本的算法3.1(也带有z更新的预条件矩阵)的收敛分析，其中建立了关于原始目标和约束满足的O(1/k)遍历率。

算法3.1中的预处理器$\frac{1}{2}||.|{|}_M^2$可以用多种方式解释。对于(25)中M的选择，其结果与在(24b)中为二次项$\frac{u}{2}||Kx+B{z}_{k+1}-v+\frac{{\lambda }_k}{u}|{|}^2$找到一个优选项的替代项一致。或者，预条件矩阵M通过在原始问题(22)中引入$\tilde{x}=M^{1/2}$，ADMM形式的冗余约束并应用原始对偶算法来求解(Nien和Fessler 2014)。

在 (Chambolle and Pock 2011) 中指出，为了最小化相同的问题模型 $g(x)+h(Kx)$，算法 3.1 的序列 $x_k$，当$ μ = σ，\tau’ = \tau$，M 在 ( 25) 与 (21) 一致。换句话说，原始对偶算法（21）可以作为算法3.1的特例获得。此外，表明（OĆonnor 和 Vandenberghe 2020）ADMM（26）和 PDHG（21）都可以作为 Douglas-Rachford 分裂（DRS）的特殊实例获得。然后来自 DRS 的收敛和收敛速率导致 ADMM 和 PDHG 的相应收敛语句。

8如果我们使用PDHG的相同问题模型(17)，则只有一个线性映射。

3.1.3.三个凸函数和的优化算法

(3.1)或(3.6)中的问题模型，具有两个凸函数和一个线性算子的和，在某种意义上可以非常限制反问题，因为我们经常需要通过将项分组并在更高维空间中定义新函数来适当地重新表示我们的目标函数(Sidky等人，2012)以符合(3.1)或(3.6)。这种重新表述通常涉及引入额外的对偶变量，这会增加内存和计算。

对于三个凸函数之和的问题，已经提出了许多算法。具体地说，它们解决了以下最小化问题：
$$\underset{x}{\min }\varphi (x)=g(x)+h(Kx)+f(x)\text{\hspace{1em}(3.10)}$$
如前所述，g和h是CCP，K是线性算子；g和h都可以是非光滑但简单的。新的分量f是CCP和$L_f$-光滑的。当f不存在时，(3.10)与(3.1)相同，并可重新表示为(3.6)中的约束形式。

在(2-块)PDHG的推导中，我们使用其共轭函数重写(3.10)中的复合形式h(K·)，然后获得(3.10)的原始对偶公式如下：
$$\underset{x}{\min }\underset{z}{\max }\tilde{\varphi }(x,z)=g(x)+<Kx,z>-h^{\ast }(z)+f(x)\text{\hspace{1em}(3.11)}$$
在(Condat 2013，Vũ2013，Chambolle and Pock 2016)中提出了求解(3.11)的扩展(3.5b)，它简单地将(3.5b)替换为以下公式：
x k + 1 : = a r g min ⁡ x { f ( x k ) + < ▽ f ( x k ) , x − x k > + g ( x ) + < K x , z k > + 1 2 τ ∣ ∣ x − x k ∣ ∣ 2 } (3.12) x_{k+1}:=arg\min_x \{f(x_k)+<\triangledown f(x_k),x-x_{k}>+g(x)+<Kx,z_k>+\frac1{2\tau}||x-x_k||^2\}\tag{3.12} xk+1​:=argxmin​{f(xk​)+<▽f(xk​),x−xk​>+g(x)+<Kx,zk​>+2τ1​∣∣x−xk​∣∣2}(3.12)
与(3.5b)相比，(3.12)中的目标函数以(2.2)的形式增加了新分量f(x)的二次上界。用新的步长建立了O(1/k)的遍历收敛速度，类似于当f=0时：
$$1/\tau -L_f>\sigma ||K|{|}^2\text{\hspace{1em}(3.13)}$$
当$L_f$=0时，即当f不存在时，它也减少到(3.5)。

算法3.2 问题(3.11)的一个原始对偶算法(严2018)
Input: Choose $x_0,z_0$, set ${\tilde{x}}_0=x_0$,set $\tau ,\sigma >0$.
Ouput: $x_K,z_K$
1 for $iter=0,...,K-1$ do
2 $z_{k+1}:=arg\min_z {<K\widetilde{x}_k,z>-h^*(z)-\frac1{2\sigma}
3 $x_{k+1}:=arg\min_x {f(x_k)+<\triangledown f(x_k),x-x_{k}>+g(x)+<Kx,z_k>+\frac1{2\tau}
4 ${\tilde{x}}_{k+1}=x_{k+1}+(x_{k+1}-x_k)+\tau (▽f(x_k)-▽f(x_{k+1}))$ /* extrapolation*/

其他直接与三个函数之和一起工作的算法可以在(Chen等人2016，Latafat和Patrinos 2017，Yan 2018)中找到。其中，(严2018)中的工作以其步长参数范围较大和每次迭代计算代价较小而值得注意。9该算法被称为算法3.2，当参数为：
$$\tau \sigma ||K|{|}^2<1,\tau L_f<2,\text{\hspace{1em}(3.14)}$$
与(3.13)相比，步长规则(3.14)消除了$\Vert K\Vert$和$L_f$对参数$\tau ，\sigma$的影响，有效地扩大了保证收敛的步长取值范围。步长取值范围的扩大是以保持$\nabla f(x)$的两个梯度向量的存储量增加为代价的，这两个梯度向量分别以$k$和$k+1$连续迭代。与基于(3.12)的3-块扩展类似，该算法在原始-对偶间隙内具有O(1/k)遍历收敛速度。当一个分量函数不存在时，算法3.2专门用于其他众所周知的两块算法，例如当f不存在时的2块PDHG(3.5)，以及当g不存在时的最近交替预测-校正器(PAPC)算法(Loris和Verhoven 2011，Chen等人2013，Droriet al 2015)。

最近，在(Davis And Yen 2017)中提出了一种三算子分裂格式${}^{10}$，作为DRS的扩展。DRS 在二算子拆分方面非常出色：它可以用来推导PDHG算法(OĆonnor和Vandenberghe2020)；当应用于约束2块问题的对偶(3.6时)，结果立即是ADMM3.8.。以类似的方式，当应用于下面的3块约束最小化问题的对偶问题时，可以使用三个算子DRS(Davis And Yen 2017)来导出如(OĆonnor和Vandenberghe2020)中所示的3块PD算法3.2。

9在梯度评价的数量方面。一些3块扩展每次迭代需要两次渐变求值，而(Yan2018)中的一次只需要一次。

10有时称为前向Douglas-Rachford分裂，因为与DRS相比，它包括额外的协迫算子(前向算子)。

$$\begin{gathered} \underset{x,y,z}{\min }{f}_1(x)+f_2(y)+f_3(z) \\ s.t.Ax+By+Cz=b\text{\hspace{1em}(3.15)} \end{gathered}$$

结果是3块ADMM，如算法3.3所示。

算法3.3 ADMM(Davis And Yen 2017)，问题(3.15a)。
Input: Choose $x_0,z_0$, set ${\tilde{x}}_0=x_0$,set $u<2\sigma/
Ouput: $x_K,z_K$
1 for $iter=0,...,K-1$ do
2 x k + 1 : = a r g min ⁡ x { f 1 ( x ) + < λ k , A x > } x_{k+1}:=arg\min_x \{f_1(x)+<\lambda_k,Ax>\} xk+1​:=argminx​{f1​(x)+<λk​,Ax>} /$f_1$ $\sigma -$strongly convex/
3 $y_{k+1}\in arg\min_y{f_2(y)+\frac u{2}
4 $z_{k+1}\in arg\min_z {f_3(z)+\frac u{2}
5 $ \lambda_{k+1}:=\lambda_k+u(Ax_{k+1}+By_{k+1}+Cz_{k+1}-b)$

算法的收敛要求$f_1(x)$是σ-强凸的，并且收敛速度继承于三个运算分裂(Davis And Yen 2017)的收敛速度O(1/k)。在实际应用中，ADMM有时以3块或多块的形式应用，在更新拉格朗日乘子之前更新三个或更多原始变量的序列。如Chen等人(2016)所示，将2块ADMM天真地扩展到3块ADMM并不一定收敛。算法3.3仅与步骤2中的这种朴素扩展不同，在步骤2中，目标函数不是增广拉格朗日，而是拉格朗日本身。

3.2.(非)光滑凸优化的加速一阶算法

在最后一节中，一个明显的疏忽遗漏是光滑最小化的经典梯度下降算法。这种疏忽遗漏是由于MBIR中广泛使用的非光滑、稀疏诱导正则化所推动的原始-对偶算法的巨大流行。然而，由于加速梯度算法(Beck和Teboulle，2009)的(重新)发现，梯度下降算法仍然是重要的，并进一步获得了动力，因为它们的收敛速度符合复杂性理论的下界(Nomerovski和Yudin 1983)，因此是最优的。这些加速梯度法反过来又促进了加速原对偶方法的发展。这些加速方法，包括原始对偶型和原始(唯一)型，将是本节的主题。

3.2.1.加速一阶原始对偶算法

随着对问题结构的更多假设，3.1节中的许多原始对偶类型算法都可以加速。例如，通过采用依赖于迭代的步长参数${\tau }_k、{\sigma }_k、{\alpha }_k$，可以如算法3.4所示加速PDHG算法(3.5)。此外，它将布雷格曼距离(钱伯勒和Pock2016)并入对偶更新方程中。${}^{11}$

算法3.4 问题(3.3)的原始对偶算法。
Input: Choose $x_0,z_0$, let ${\tilde{x}}_0=x_0,{\tau }_k>0,{\sigma }_k>0$,s.t. $\sigma_0 \tau_0
Output:$x_K,z_K$
1 for $iter=0,...,K-1$ do
2 $z_{k+1}:={\arg max}_z\{⟨K{\tilde{x}}_k,z⟩-h^{\ast }(z)-\frac{1}{{\sigma }_k}{D}_2(z;z_k)\}$ /* dual ascent*/
3 x k + 1 : = arg ⁡ m i n x { g ( x ) + ⟨ K x , z k + 1 ⟩ + 1 2 τ k ∥ x − x k ∥ 2 } {x}_{k+1}:= {\arg min}_{x}\{g(x)+ \langle {Kx},{z}_{k+1} \rangle +{\textstyle \tfrac{1}{2{{\tau }}_{k}}}\parallel x-{x}_{k}{\parallel }^{2}\} xk+1​:=argminx​{g(x)+⟨Kx,zk+1​⟩+2τk​1​∥x−xk​∥2} /* primal descent*/
4 ${\tilde{x}}_{k+1}:=x_{k+1}+{\alpha }_k(x_{k+1}-x_k)$ /* extrapolation

在(Chambolle And Pock2016)中证明了，如果g是$\gamma$-强凸的，通过设置参数${\sigma }_k+1=\frac{{\sigma }_k}{{\alpha }_k},{\tau }_{k+1}={\alpha }_k{\tau }_k,{\alpha }_k=\frac{1}{\sqrt{1+2\gamma {\tau }_k}}$，可以将算法的收敛速度提高到O($1/k^2$)，其中$\gamma$是g的强凸性参数。

实现加速的另一种方法是利用不同算法之间的联系，而不是从头开始重新推导。如3.1节中所讨论的，DRS可用于导出PDHG算法(OĆonnor和Vandenberghe 2020)；这种关联可用于从加速的DRS导出加速的PDHG算法(Davis And Yen 2017)。同样，由于预条件ADMM(算法3.1)等价于应用于对偶问题的PDHG，因此可以从加速的PDHG(算法3.4)得到预条件ADMM的加速版本。

11这个版本的算法(Chambolle And Pock 2016)比(Chambolle And Pock 2011)中提出的算法稍微更通用一些

同样的策略也适用于3块算法。如(OĆonnor and Vandenberghe2020)所示，3-算子分裂DRS和3-块原始-对偶算法3.2之间的等价性，意味着算法3.2的加速版本可以从已经完成的加速3-算子分裂(Davis And Yen 2017)得到(Condatet Al 2020)。

在这些原始-对偶加速算法中，一个共同的假设是目标函数要么是强凸的，要么是L-光滑的，以实现从O(1/k)到O(1/$k^2$)的加速。如果目标函数同时包含线性分量中的光滑分量(具有Lipschitz连续梯度)和非光滑分量，则这些算法的收敛速度将由非光滑分量控制，最多为O(1/k)。

这种情况并不令人满意，确实可以改善。如(Nester Ov 2005)所证明的，有可能获得“模块化”的最优收敛速度，其对目标函数的光滑分量具有O(1/$k^2$)依赖关系，对(结构化)非光滑分量具有O(1/k)依赖关系。虽然这种算法的整体收敛速度仍然由O(1/k)控制，但它能更好地处理问题模型中的大梯度Lipschitz常数，这可能是许多成像反问题的情况。加速的原始对偶算法(Chen等人，2014)和加速的ADMM算法(欧阳等人，2015)也达到了这种“最优”的收敛速度。

3.2.2.加速(近似)梯度下降(AGD)算法

许多关于加速一阶方法的工作都是受Nester ov 1983年的开创性论文(Nester Ov 1983)的启发，这篇论文的最简单形式是考虑最小化f(x)的问题，其中f(x)是$L_f$光滑的。对于这类问题，著名的标准梯度下降算法$x_{k+1}=x_k-1/L_f\nabla f(x_k)$以$O(L_f/k)$的速度收敛于目标值，即$f(x_k)-f(x_{\ast })<=O(L_f/k)$，其中$x\in arg{\min }_{x}f(x)$是假定存在的。Nesterov法展示了以下两步序列：
x ~ k + 1 = a r g m a x y { < ▽ f ( y k ) , y − y k > + L f 2 ∣ ∣ y − y k ∣ ∣ 2 } = y k − 1 L f ▽ f ( y k ) y k + 1 = x ~ k + 1 + θ k + 1 θ k ( 1 − θ k ) ( x ~ k + 1 − x ~ k ) , k = 0 , 1 , . . . (3.16) \tilde{x}_{k+1}= arg max_{y}\{ <\triangledown f(y_k),y-y_k>+\frac{L_f}{2}||y-y_k||^2\}=y_k-\frac1{L_f}\triangledown f(y_k)\tag{3.16}\\ y_{k+1}=\tilde{x}_{k+1}+\frac{\theta_{k+1}}{\theta_k}(1-\theta_k)(\tilde{x}_{k+1}-\tilde{x}_{k}),k=0,1,... x~k+1​=argmaxy​{<▽f(yk​),y−yk​>+2Lf​​∣∣y−yk​∣∣2}=yk​−Lf​1​▽f(yk​)yk+1​=x~k+1​+θk​θk+1​​(1−θk​)(x~k+1​−x~k​),k=0,1,...(3.16)
以及一个复杂的内插参数序列。${}^{12}$
$$\frac{1-{\theta }_{k+1}}{{\theta }_{k+1}^2}=\frac{1}{{\theta }_k^2},{\theta }_0=1\text{\hspace{1em}(3.17)}$$
导致$O(L_f/k^2)$的$f({\tilde{x}}_k)$的收敛速度加快。就其对k和$L_f$的依赖性而言，该速率是最优的，即不可侵犯的，因为它匹配用于仅使用一阶信息最小化平滑函数的较低复杂度界限。

Nesterov的论文(NesterOv 1983)也考虑了约束极小化问题${\min }_{x\in C}f(x)$，其中C是闭凸集。通过将(3.16a)替换为梯度投影步长，即${\tilde{x}}_{k+1}=pro{j}_C(y_k-\frac{1}{L_f}▽f(x_k)$，即可得到解，其中$pro{j}_C$是凸集C上的正交投影。(3.16)的这个受限版本可以被视为著名的Fista(Beck和Teboulle 2009)的前身。

在过去的十年左右的时间里，人们对Nesterov的加速算法进行了广泛的分析，提出了许多变体。一个这样的变体，算法3.5，例如参见(Auslender and Teboulle 2006，Tseng 2008)，考虑最小化复合目标函数$f(x)=f(x)+g(x)$，其中f如前所述是$L_f$光滑的，g是简单的，并且假设存在$x_{\ast }=argmin\varphi (x)$。

算法3.5 min f+g，f是$L_f$光滑的，g是单的。
Input: Choose ${\tilde{x}}_0=x_0$,and let ${\theta }_k$ follow(3.17)
Output:$x_K$
1 for $iter=0,...,K-1$ do
2 $ y_{k}=(1-\theta_k)\tilde{x}{k}-\theta_k{x}{k}$
3 x k + 1 = a r g max ⁡ x { g ( x ) + f ( y k ) + < ▽ f ( y k ) , x − y k > + θ k L f D ( x ; x k ) } x_{k+1}= arg\max_x \{g(x)+f(y_k)+<\triangledown f(y_k),x-y_k>+\theta_k L_fD(x;x_k)\} xk+1​=argmaxx​{g(x)+f(yk​)+<▽f(yk​),x−yk​>+θk​Lf​D(x;xk​)}
4 ${\tilde{x}}_{k+1}=(1-{\theta }_k){\tilde{x}}_k-{\theta }_k{x}_{k+1}$

12 '=‘可通过（3.17）可以用’<='代替，参见例如（Tseng 2008）。例如，θk=2/(k+2)满足（Nesterov 2005）中使用的不等式。通过这个选择，外推步（3.16b）被简化为$y_{k+1}={\tilde{x}}_{k+1}+\frac{k}{k+3}({\tilde{x}}_{k+1}-{\tilde{x}}_k)$。

请注意，算法3.5维护三个序列，${\tilde{x}}_k、x_k$和$y_k$，这比两个序列更新方程(3.16)更复杂。然而，与仅限于二次距离的(3.16a)不同，增加的复杂性是由于梯度下降步骤(3)包含Bregman距离的灵活性而得到的。当g(x)不存在，且$D(x;y_k)=1/2\Vert x-y_k{\Vert }^2$时，可以证明算法3.5的序列${\tilde{x}}_k、y_k$与(3.16)重合。类似于(3.16)，算法3.5中${\tilde{x}}_k$的收敛速度满足$\varphi ({\bar{x}}_k)-\varphi (x_{\ast })⩽\mathcal{O}(L_f/k^2)$。

算法3.4和算法3.5之间的一种有趣的等价关系在(Lan和周2018)中被发现，使用算法3.4的对偶上升步中的Bregman距离的特殊化。${}^{13}$设算法3.4的对偶上升步骤中的$D_2(x;z_k)$是由h*自身产生的Bregman距离，即，
$$D_2(z;z_k)=h^{\ast }(z)-[h^{\ast }(z_k)+<▽h^{\ast }(z_k),z-z_k>]\text{\hspace{1em}(3.18)}$$
然后双上升步骤变成：
z k + 1 = a r g max ⁡ z { < K x ˉ k , z > − h ∗ ( z ) − 1 σ k D 2 ( z ; z k ) } ( 3.18 ) —— > z k + 1 = = a r g max ⁡ z { < K x ˉ k + 1 σ k h ∗ ( z k ) , z > − ( 1 + 1 σ k ) h ∗ ( z ) } ( a ) —— > = a r g max ⁡ z { ( w k + 1 , z ) − h ∗ ( z ) } ( 2.7 ) —— > = ▽ h ( w k + 1 ) (3.19) z_{k+1}=arg\max_z\{<K\bar x_k,z>-h^*(z)-\frac1{\sigma_k}D_2(z;z_k)\}\\ (3.18)——>z_{k+1}==arg\max_z\{<K\bar x_k+\frac1{\sigma_k}h^*(z_k),z>-(1+\frac1{\sigma_k})h^*(z)\}\\ (a)——>=arg\max_z\{(w_{k+1},z)-h^*(z)\}\\ (2.7)——>=\triangledown h(w_{k+1})\\ \tag{3.19} zk+1​=argzmax​{<Kxˉk​,z>−h∗(z)−σk​1​D2​(z;zk​)}(3.18)——>zk+1​==argzmax​{<Kxˉk​+σk​1​h∗(zk​),z>−(1+σk​1​)h∗(z)}(a)——>=argzmax​{(wk+1​,z)−h∗(z)}(2.7)——>=▽h(wk+1​)(3.19)
其中，在a中，我们将$w_{k+1}$定义为下划线术语的缩放版本：
$$w_{k+1}:=\frac{K{\bar{x}}_k+{\sigma }_k^{-1}▽h^{\ast }(z)}{1+{\sigma }_k^{-1}}-(3.19)->=\frac{K{\bar{x}}_k+{\sigma }_k^{-1}{w}_k}{1+{\sigma }_k^{-1}}\text{\hspace{1em}(3.20)}$$
将(3.19)、(3.20)与算法3.4相结合，然后由下式给出专门的原始-对偶更新步骤：
w k + 1 : = K x ˉ k + σ k − 1 w k 1 + σ k − 1 x k + 1 = a r g m a x { g ( x ) + < K x , ▽ h ( w k + 1 ) > + 1 τ k D 1 ( x ; x k ) } x ~ k + 1 = x k − α k ( x k + 1 − x k ) (3.21) w_{k+1}:=\frac{K\bar x_k+\sigma^{-1}_kw_k}{1+\sigma^{-1}_k}\\ x_{k+1}=argmax\{g(x)+<K x,\triangledown h(w_{k+1})>+\frac1{\tau_k}D_1(x;x_k)\}\\ \tilde x_{k+1}=x_k-\alpha_k(x_{k+1}-x_k)\tag{3.21} wk+1​:=1+σk−1​Kxˉk​+σk−1​wk​​xk+1​=argmax{g(x)+<Kx,▽h(wk+1​)>+τk​1​D1​(x;xk​)}x~k+1​=xk​−αk​(xk+1​−xk​)(3.21)
将算法3.5的f(x)识别为求解h(Kx)+g(x)的PDHG算法(算法3.4)中的h(Kx)，附录A.3中的进一步处理表明两个算法的参数可以匹配，使得(3.21b)中的序列$x_k$与来自算法3.5的序列$x_k$一致。从算法3.5的第3行开始，${\bar{x}}_k$和$x_k$之间的关系是${\bar{x}}_k$是$x_k$的加权平均值。算法3.5中${\bar{x}}_k$以O($1/k^2$)的速率收敛，然后转化为(加权的)$x_k$在相同速率下的遍历收敛，这与算法3.4中的结论相同。

13严格地说，(Lan和周2018)中建立的关系是关于算法3.4的一个变体，该算法允许Bregman距离出现在原始和对偶更新方程中。有关更多详细信息，请参阅(兰和周2018)。

3.3.一阶算法在成像问题中的应用

在这一节中，我们将讨论如何使用前几节的算法来解决反问题。我们首先定义了一个通常用于CT重建的原型问题。然后，我们将一些有代表性的算法应用于原型问题。通常需要将我们的问题重新表述为模型形式((3.1)、(3.6)或(3.11))。我们探索了不同的选择，并讨论了相关的内存和计算成本。

3.3.1.问题定义

CT重建通常可以表示为以下最小化问题：${}^{14}$

$$\underset{x}{\min }\varphi (x),\varphi (x)=\frac{1}{2}||y-Ax|{|}_w^2+H(w)+G(x)\text{\hspace{1em}(3.22)}$$
其中$y\in R^m$是测量的投影数据，$A\in R^{m\times d}$是系统矩阵或前向投影算子，$0<w<R^m$是与投影数据y相关的统计权重，$x\in R^d$是要重建的未知图像。设$$x_{\ast }\in arg\underset{x}{\min }\varphi (x)$$，我们总是假设$x^{\ast }$存在。

14CT中常用的是二次数据拟合模型。对于PET和SPECT重建，数据拟合项通常是负Poisson对数似然，其梯度不是(全局)Lipschitz连续的。更多详情见附录A.2。

在不失一般性的情况下，我们假设统计权重被按缩放，使得对于$0<w_j<=1,j=1,2,...,m$。${}^{15}$比例因子可以被吸收到正则化的定义H(x)和G(x)中，它们编码了我们关于x的先验知识。在这里，我们假设G是一个简单的函数，而H不是。H(X)在压缩感知中的一个流行的例子是TV正则化，由下式给出
$$H(x)=\sum _i\tilde{H}(K_{i}x)\text{\hspace{1em}(3.23)}$$
其中$K_{i}x=[x_i-x_{i,i_1},x_i-x_{i,i_2},x_i-x_{i,i_3}]$，$i=1,...,d$是有限差分算子，$x_{i,i_2}，x_{i,i_2}，x_{i,i_3}$表示$x_i$的三维邻域。如果$\tilde{H}(z)={\sum }{j}| {z}{j}| $，则H(x)是各向异性TV；如果$\tilde{H}(z)=\parallel z\parallel $，则H(x)是各向同性TV。

通过将$K_i$指定为通用线性算子，例如(学习的)卷积滤波，并且通过将H(x)指定为可以是(非)光滑或(非)凸的通用势函数，(3.23)中H(x)的简单表达式确实可以包含多种正则化。(3.22)中的最后一项G(x)编码了对未知x的简单(稀疏)约束。例如，有时将x限制在凸集C中是有物理意义的，当x表示人体的线性衰减系数时，则C是非负的正交体。在这种情况下，$Gx)=I_C(x)$。为方便起见，我们还使用$F(x)=\parallel y-Ax{\parallel }_w^2/2$来表示(3.22)中的数据拟合项。

3.3.2.使用两块PDHG算法(3.5)

在CT重建的背景下，正则化H(x)可以是(非)光滑的，并且可能经常涉及线性算子，例如有限差分算子。因此，我们很自然的将原型问题转换成问题(3.1)：

$$\begin{gathered} F(x)+G(x)<-->g(x) \\ H(x)=\sum _i\tilde{H}(K_{i}x)<-->h(Kx)\text{\hspace{1em}(3.24)} \end{gathered}$$
根据双共轭关系(3.2)，我们可以写：
$$\sum _i\tilde{H}(K_{i}x)=\sum _i(\underset{z_i}{\max }<K_{i}x,z_i>-{\tilde{H}}^{\ast }(z_i))\text{\hspace{1em}()}$$
其中对偶变量$z_i,i=1,...,n$是可分的。这种重新表述会导致以下更新公式。对应于(3.5a)和(3.5b)：

1.双重更新:
$$z_{k+1}=arg\underset{z}{\max }\sum _i(<K_{i}x,z_i>-{\tilde{H}}^{\ast }(z_i)-\frac{1}{2\sigma }||z_i-z_{i,k}|{|}^2)\text{\hspace{1em}(3.25)}$$
请注意，最大化问题在$z_i$中是可分离的，因此可以并行完成。这个更新本质上需要计算$pro{x}_{\sigma {\tilde{H}}^{\ast }}$，这很容易用莫罗恒等式(2.12)和我们假设$\tilde{H}$是简单的来计算。

2.原始更新:
x k + 1 = a r g min ⁡ x { F ( x ) + G ( x ) + ∑ i < K i x , z i , k + 1 > + 1 2 τ ∣ ∣ x − x k ∣ ∣ 2 } (3.26) x_{k+1}=arg\min_x\{F(x)+G(x)+\sum_i<K_ix,z_{i,k+1}>+\frac{1}{2\tau}||x-x_{k}||^2\}\tag{3.26} xk+1​=argxmin​{F(x)+G(x)+i∑​<Ki​x,zi,k+1​>+2τ1​∣∣x−xk​∣∣2}(3.26)
同样，此更新需要计算F(x)+G(x)的近端映射。由于F(x)是数据拟合项，无论G(x)是简单的，该更新可能不能以封闭形式计算或以其他方式有效地获得。作为一种实用的选择，$x_{k+1}$通常通过运行(近端)梯度下降算法的几次迭代来近似。在绝对可求和误差的条件下，尽管更新的近似性质，理论收敛结果仍然可以建立。${}^{16}$

或者，我们可以使用加权二次差分来应用一般的近端映射步骤，${}^{17}$类似于我们在预条件ADMM中所做的那样(参见(3.9))，即，

15在第4.5节中需要进行此缩放，其中权重以Bregman距离表示。

16详情见第3.4节讨论。

17提出了只使用二次距离的两块PDHG算法；在算法的非加速版本中，PDHG的三块扩展包含了原始更新和对偶更新的Bregman距离。

x k + 1 = a r g min ⁡ x { F ( x ) + G ( x ) + ∑ i < K i x , z i , k + 1 > − 1 2 ∣ ∣ x − x k ∣ ∣ M 2 } (3.27) x_{k+1}=arg\min_x\{F(x)+G(x)+\sum_i<K_ix,z_{i,k+1}>-\frac{1}{2}||x-x_{k}||^2_M\}\tag{3.27} xk+1​=argxmin​{F(x)+G(x)+i∑​<Ki​x,zi,k+1​>−21​∣∣x−xk​∣∣M2​}(3.27)

由于$F(x)=\frac{1}{2}||y-Ax|{|}_w^2$，如果我们选择M为：
M = 1 τ I − A t d i a g { w } A (3.28) M=\frac{1}{\tau}I-A^tdiag\{w\}A\tag{3.28} M=τ1​I−Atdiag{w}A(3.28)
并且$$\tau$$，满足 1 τ > ∣ ∣ A t d i a g { w } A ∣ ∣ + σ ∑ i ∣ ∣ K i ∣ ∣ 2 \frac{1}{\tau}>||A^tdiag\{w\}A||+\sigma\sum_i||K_i||^2 τ1​>∣∣Atdiag{w}A∣∣+σ∑i​∣∣Ki​∣∣2(参见(3.13))，然后插入F(x)和M到(3.27)：
$$\begin{gathered} F(x)+\frac{1}{2}||x-x_k|{|}_M^2=\frac{1}{2}||y-Ax|{|}_w^2+\frac{1}{2}(x-x_k{)}^t(\frac{1}{\tau }I-A^{t}wA)(x-x_k) \\ =-<y-A{x}_k,wA(x-x_k)>+\frac{1}{2\tau }||x-x_k|{|}^2+constant \\ =<▽F(x_K),(x-x_K)>+\frac{1}{2\tau }||x-x_k|{|}^2+constant \\ \text{\hspace{1em}()} \end{gathered}$$
则(3.27)中的$x_{k+1}$允许闭合形式的解：
x k + 1 = a r g m i n { G ( x ) + ∑ i < K i t z i , k + 1 + ▽ F ( x k ) , x − x k > + 1 2 τ ∣ ∣ x − x k ∣ ∣ 2 } = p r o x τ G x K − τ [ ∑ i < K i t z i , k + 1 + ▽ F ( x k ) ] (3.29) x_{k+1}=argmin\{G(x)+\sum_i<K_i^tz_{i,k+1}+\triangledown F(x_k),x-x_k>+\frac{1}{2\tau}||x-x_k||^2\}\\\tag{3.29} =prox_{\tau G}{x_K-\tau[\sum_i<K_i^tz_{i,k+1}+\triangledown F(x_k)]} xk+1​=argmin{G(x)+i∑​<Kit​zi,k+1​+▽F(xk​),x−xk​>+2τ1​∣∣x−xk​∣∣2}=proxτG​xK​−τ[i∑​<Kit​zi,k+1​+▽F(xk​)](3.29)
综上所述，我们选择了一个特殊的预条件矩阵M，它消除了数据拟合函数F(x)中的二次项，得到了闭合形式的原始更新$x_{k+1}$。

3.3.3.使用三块PD算法3.2

由于算法3.2直接处理三个函数的和(3.10)，因此我们的原型问题(3.22)和(3.10)之间的自然对应关系如下：
$$
F(x)<–>f(x)\\tag{}

G(x)<–>g(x)\

H(x)=\sum_i\tilde H(K_ix)<–>h(kx)\
$$
算法通过依次计算F(x)的梯度以及G和${\tilde{H}}^{\ast }\ast$的近邻映射来进行，这些都是容易计算的。更新方程类似于(3.25)和(3.29)，并且具有不同的外推步骤(算法3.2的第4行)，其中应用了梯度校正。收敛的步长要求是：$\sigma \tau {\sum }_i\Vert K_i{\Vert }^2$等于1，且$\tau \Vert A^{t}wA{\Vert }^2<2$。

3.4.讨论

我们讨论了获得最优收敛速度的一阶算法的加速变体，例如，对于光滑优化，改进为O(1/k)到O($1/k^2$)。除了这些技术外，还经常通过过度松弛来经验地观察到加速度。给定形式为$x_{k+1}=T{x}_k$的定点迭代，过度松弛指的是通过以下方式更新$x_k$：
$$x_{k+1}=x_k+\rho (T{x}_k-x_k)\text{\hspace{1em}(3.30)}$$
其中${\rho }_k$是(迭代相关的)超松弛参数。超松弛不动点迭代(3.3)收敛于α平均算子，其形式为$T$=$T_{\alpha }$=$(1-\alpha )Id+\alpha N$，其中Id是单位映射，N是非扩张映射，0<α<1。如果算子T是1/2平均的，即α=1/2，则松弛参数${\rho }_k$≡ρ可以逼近2，而不动点迭代(3.30)仍然是平均算子，因此仍然保证(3.3)的收敛。

我们讨论的许多迭代算法都是α平均算子。L-光滑函数f的简单梯度下降算法$Tx=x-\frac{1}{L}\nabla f(x)$平均是1/2；(2-块)PDHG算法(具有θ=1)和ADMM算法是近似点算法的实例，其平均是1/2；用于最小化三个函数和的Yen的算法(Yan 2018)和Davis-Yenʼ的三个算子分裂(Davis和Yen 2017)也是平均算子。如果选择适当的超松弛参数${\rho }_k={\rho }_k(\alpha )$，所有这些算法都可以有像(3.30)这样的超松弛形式，并保证收敛。过度松弛的理论证明确实表明收敛界可以减少ρ>1，从O(1/k)到O(1/ρk)，例如见(Chambolle And Pock2016)，定理2。

正如我们在第3.3节中遇到的那样，有时很难准确地计算$T{x}_k$，例如，什么时候T是复函数的近似映射。不精确Krasnoselskii-Mann(Km)定理考虑以下形式的不精确更新：$x_{k+1}=T_{{\alpha }_k}{x}_k+{\alpha }_k{ϵ}_k$其中$T_{{\alpha }_k}={\alpha }_{k}N+(1-{\alpha }_k)\mathrm{Id}$是${\alpha }_k$-平均算子，而${\alpha }_k{ϵ}_k$量化更新$x_{k+1}$中的错误。如果错误满足${\sum }_k{\alpha }_k\Vert {ϵ}_k\Vert <\infty$和${\sum }_k{\alpha }_k(1-{\alpha }_k)\to \infty$，则迭代$x_k$仍然收敛到N的不动点(梁等人，2016)。对于过于宽松的版本(3.30)，适当选择松弛参数${\rho }_k$，不动点迭代(3.30)保持平均，并且不精确的KM定理仍然适用。

上一节中的例子展示了将一阶算法应用于CT图像重建所涉及的典型步骤：问题的重构和子问题的求解通常都需要特定问题的工程工作。此外，开发这样的算法还需要大量的研究人员时间。从实践者的角度来看，解决定义明确的优化问题的理论保证应该与这些努力背后的开发时间进行权衡。如果人们愿意放弃算法的精确度，那么可以通过优化法获得启发式解(Herman等人2012年，Censoret等人2017年)。

超优化适用于复合最小化问题，其中扰动弹性算法被引导到减少正则化函数，同时保持与数据保真度诱导约束的兼容性。超优化可以成为一个自动过程，将一个算法转化为它的上级版本，从而可以最大限度地减少算法开发和实现的研究时间。与我们在本章中讨论的精确算法不同，优化法是启发式的，因为不能保证结果接近目标函数的最小值。关于这一方法的更多信息可以从最初的提出者之一(Censor 2021)维护的书目网站上找到。