数学工具：图像重建技术

此处提供本书中的算法所涉及的一些数学工具，主要涉及矩阵论、凸分析中的一些基本概念与结论，最后介绍两个在成像问题中被广泛应用的最优化算法——ADMM［241］和FIS⁃TA［200］。

A1　矩阵论

A1.1　矩阵子空间

行空间：矩阵A的行空间是由A的行向量张成的子空间。

列空间：矩阵A的列空间是由A的列向量张成的子空间。

零空间：矩阵A的零空间是所有满足Ax＝0的x的集合。

秩：rank A是线性独立的列或行向量的数量。矩阵A∈Rm×n是满秩的，如果rank A＝min｛m，n｝。

正交子空间：两个子空间S1，S2∈Rm×n是正交的，如果，则∀s1∈S1，s2∈S2。

A1.2　矩阵分解

特征分解：实对称矩阵A∈Rn×n可以分解为

此处，Q为单位正交矩阵；QTQ＝I；Λ＝diag（λ1，λ2，…，λn），特征值λi按降序排列。此外：

•A当且仅当其所有特征值非零时可逆，A－1＝QΛ－1QT；

•A是半正定的，如果其所有特征值非负。

奇异值分解：任何rank r的矩阵A∈Rm×n可以分解为

此处，U∈Rm×r满足UTU＝I，V∈Rn×r满足VTV＝I。Σ＝diag（σ1，σ2，…，σn），奇异值σi降序排列。此外：

•ATA的特征分解，矩阵W使得［VW］为正交阵；

•A的条件数cond A＝σ1／σr；

•A的伪逆A†＝VΣ－1UT。

A1.3　范数

•‖x‖0，x中的非零元素个数；尽管这不是一个严格意义上的范数，但是在作为一个稀疏规整化项时，有着广泛应用。

矩阵的Lp范数：

谱范数／算子范数：‖X‖2＝σ1（X），X的最大奇异值。

迹：X的所有奇异值之和。

A1.4　矩阵微分

最常用的矩阵微分有：

A1.5　Weyl定理

Weyl定理表述的是一个厄米特矩阵被添加一个加性扰动时，其特征值的变化情况［315］。令｛A′，A，ΔA｝为三个N×N的厄米特矩阵，其特征值｛λm（A′），λm（A），λm（ΔA）｝降序排列，即

同时，｛A′，ΔA｝的特征值也如此排列，下标1代表最大特征值，而N代表最小特征值。Weyl定理表明，如果矩阵A被扰动，即

则对所有1≤n≤N，新矩阵A′的特征值有以下限制。

特别地，最大特征值按如下方式扰动。

如果扰动矩阵为半正定的，即ΔA≥0，则由（A－3）可知，对所有1≤n≤N，有

A1.6　盖尔圆定理（Gershgorin's Theorem）

盖尔圆定理给出了一个矩阵的特征值所在的圆形区域范围。考虑一个N×N的矩阵A，在复平面上，其对角线元素aii与一个中心位于aii的圆相关联，其半径为

也就是说，ri是与aii同一行上所有非对角线元素的模之和。将这个圆记作Di，它由所有满足以下条件的点构成。

该定理表明，矩阵A的谱［即A的所有特征值的集合，记作λ（A）］包含在所有N个盖尔圆的并集之中。

当一些盖尔圆彼此分离时，还有一个关于盖尔圆定理的更强的论述：如果有L个盖尔圆的并集与另外N－L个盖尔圆的并集不重叠，那么该定理保证A的L个特征值就位于这L个盖尔圆的并集之中，而剩余的N－L个特征值位于另外的N－L个盖尔圆的并集之中。

A1.7　加权二次函数的代数求解

考虑如下加权二次函数的最优化问题。

式中，f是二次函数。

而A∈Rp×n，S∈Sn＋为对称正定矩阵，x∈Rn，v∈Rp。假设P＋ρATA是可逆的，则x＋为v的仿射函数。

因此，x＋的计算即求解一个线性系统，其系统矩阵为一个正定矩阵P＋ρATA，而常向量为ρATv－q。进一步利用矩阵结构，可以显著提升计算速度。这一求解过程称为反向求解（back－solve）。该方法中矩阵分解与反向运算的计算复杂度取决于矩阵F的稀疏性和其他特性。

A1.7.1　直接求解

直接求解一个线性系统Fx＝g基于以下方式：首先进行矩阵分解，将矩阵分解为一系列简单矩阵乘积，F＝F1F2…Fk，然后通过一系列简单问题的求解最终得到x＝zk。

计算F＝P＋ρATA的代价为O（pn2）；之后进行Cholesky分解的代价是O（n3）；反向求解的代价是O（n2）；计算g＝ρATv－q的代价为O（np），相对前面的代价可忽略不计。当p和n是同一量级或比n大时，整体的计算复杂度是O（pn2）；当p的量级小于n时，其后的矩阵逆引理的复杂度可降为O（p2n）。

A1.7.2　利用稀疏性

当A和P的结构使F是稀疏的（F中大部分元素为零）时，还存在更高效的矩阵分解与反向求解方法可以利用。一个极端的例子是，A和P是n×n对角矩阵，此时矩阵逆的计算复杂度是O（n），整体的求解复杂度由计算g＝ρATv－q的代价O（np）占主导。如果A和P都是带状的，则F也是带状的。记F的带状宽度为k，则矩阵分解代价是O（nk2），而反向求解的代价是O（nk）。这样，式（A－8）的求解代价为O（nk2）再加上计算F的代价O（pk2）。当F＝P＋ρATA具有更一般的稀疏分布特点时，可以利用置换Cholesky分解，而置换操作的目的就是使分解后的矩阵尽量稀疏。

A1.7.3　缓存矩阵分解

现在考虑求解多个线性系统Fx（i）＝g（i），i＝1，2，…，N，它们具有相同的系数矩阵F但是不同的右侧向量g（i）。由于系数矩阵F是固定的，故其分解可以只计算一次并存储起来，从而使反向求解每次只针对不同的g（i）进行。如果将矩阵分解的代价记作t，而反向求解的代价记作s，则总代价是t＋Ns，而不是每次都进行矩阵分解的代价N（t＋s）。根据A1.7.1节中的讨论，使得在每次反向求解x（i）时的代价为O（n2），比直接求解O（n3）效率提高了n倍。

如果进一步考虑A和P的特殊结构，t和s的比例一般会小于n，但是仍然会比较可观，从而带来显著计算效率的提升。但是，当变换参数ρ时，需要权衡变化ρ带来的好处和不重新分解F＝P＋ρATA二者之间的收益。通常，在实现过程中，可以在需要改变ρ时重新分解F，而在其他时候存储分解结果用于反向求解。

A1.7.4　矩阵逆引理

在递推最小二乘法估计问题中，因为每次推导运算时必须计算矩阵和的逆，这样做工作量非常大，为了简化，通常使用矩阵求逆引理来简化计算量。矩阵求逆引理要解决的问题是：

已知一个高维矩阵A的逆矩阵，当A矩阵产生了一个非常小的变化（秩低于rank A）时，能不能根据已知的A的逆矩阵，求发生微小变化后的矩阵的逆。矩阵逆引理如下：

若矩阵A，C∈CN×N均为非奇异矩阵，矩阵B∈CN×M，D∈CM×N，则矩阵A＋BCD具有逆矩阵：

可以利用上面的矩阵逆引理来处理式（A－8）中的（P＋ρATA）－1，有

这意味着如果矩阵P－1容易计算，而且p比n小或者至少在数量级上不大于n，式（A－8）的求解可以更有效地通过A1.7.3节中的矩阵分解与缓存方法来实现：首先分解I＋ρAP－1AT，注意它是p×p的，而此处假定p≤n，然后缓存矩阵分解并将其用于反向求解。

这里讨论一个特殊情况：矩阵P为对角的，而且p≤n。按照A1.7.3节中的方法，直接计算矩阵分解的代价是O（n3），而存储P＋ρATA的分解矩阵后的反向求解代价为O（n2）。使用矩阵逆引理的话，矩阵分解的代价为O（np2）——这是因为当P为对角矩阵时，矩阵逆引理中的主要代价是计算AP－1AT、O（np2），而对I＋ρAP－1AT进行分解的代价是O（p3）＜O（np2）。可见此种情况下，利用矩阵逆引理进行矩阵分解，可以获得比A1.7.3节中的方法高（n／p）2倍的增益；在存储矩阵分解结果后，反向求解代价是O（np），可以获得比A1.7.3节中的方法高（n／p）倍的增益。

利用矩阵逆引理，在变化参数ρ时x＋的计算代价更低。例如，当矩阵P为对角的，可以只计算一次AP－1AT，代价为O（np2）；然后在第k步迭代时计算和分解I＋ρAP－1AT的代价为O（p3）。当p2在量级上小于等于n时，计算一次AP－1AT的代价O（np2）小于计算和分解P＋ρATA的代价O（n3）；而每一步反向求解时，A1.7.3节中的方法代价为O（n2），这使得O（p3）＋O（np）＜O（n2）。也就是说，相对于直接分解存储法，矩阵逆引理在此条件下，即便是使用变化的参数ρ，仍然不增加计算量。

A2　凸分析基础

A2.1　凸函数

一个拓展实值函数称之为凸的，如果它的上图（epigraph）：

是凸集，即如果λ≥f（x），μ≥f（y），t∈［0，1］，则有tλ＋（1－t）μ≥f（tx＋（1－t）y）。

特展（proper）函数：f不取值－，而且不取值为＋。

对于特展函数f，其为凸函数当且仅当对所有x，y∈X，t∈［0，1］，有

严格凸（strictly convex）：上面的不等式（A－10）对任何x≠y，0＜t＜1都严格成立（不取等号）。

下准连续（lower semi－continuous，l.s.c.）：对所有x∈X，xn→x，有

一个重要的例子是集合C的示性函数（indicator function）。

当C为闭合非空凸集，其示性函数为特展l.s.c.凸函数。

A2.2　次梯度（subgradient）

对于一个实值特展l.s.c.凸函数f：X→［－，＋］，其在点x∈X的次梯度定义为

次梯度可以用来判定非平滑凸函数的最优化条件：

x∈X是f的全局最小值点⇔0∈∂f（x）

强凸性（strongly convex，μ－convex）：对于x，y∈X，p∈∂f（x），有

一个强凸函数显然也是严格凸的，因为对任何x，y∈X，t∈［0，1］，有

如果f是强凸的，并且x是最小值点（0∈∂f（x）），则对所有y∈X，有

定义域（domain）：f的定义域是集合dom f＝｛x∈X：f（x）＜＋｝，而∂f的定义域是集合dom∂f＝｛x∈X：∂f（x）≠∅｝。

A2.3　Lengendre－Fenchel共轭（conjugate）

对任意函数f：X→［－，＋］，其Lengendre－Fenchel共轭（凸共轭）为

f∗（y）是线性连续函数族〈y，x〉－f（x）的上确界，因此是l.s.c.凸函数。

双重共轭f∗∗是在f之下的最大的l.s.c.凸函数，f∗∗≤f。当f为凸函数时，f∗∗＝f。由此梯度定义式（A－13）可知，当且仅当y∈∂f（x）时，x达到式（A－16）中的上确界，从而有

f（x）＋f∗（y）＝〈y，x〉

此情况下，f∗∗（x）＝f（x）＝〈y，x〉－f∗（x），也就是y使得f∗∗（x）达到Lengendre－Fenchel共轭定义中的上确界，而且x∈∂f∗（y）。这即引出著名的Lengendre－Fenchel等式：

由定义可知，一个凸函数f的次梯度∂f是单调算子。

若f是强凸的，则∂f是强单调的，即

一个l.s.c.凸函数f是μ－强凸的，当且仅当它的共轭∂f∗是1阶连续且有（1／μ）－Lipschitz梯度（gradient）。

则称f（x）为梯度μ－Lipschitz连续的。

极大单调算子（maximalmonotone operator）：对于多值算子T：X→P（X），满足

且其图（graph）｛（x，p）：x∈Tp｝⊂X×X在所有满足式（A－21）的算子的图中是极大的。可见，次梯度算子是极大单调算子的一个特例。而且，根据定义，任何极大单调算子T都是可逆的，定义为x∈T－1p⇔p∈Tx。相应的，在此意义下∂f和∂f∗是互逆的。

A2.4　迫近图（Proimalmap）与预解算子（Resolvent）

在最优化理论中，一个具有重要作用的概念是迫近图或迫近算子（proximity operator）。如果f是特展l.s.c.凸函数，则对任何x，都存在如下强凸问题的唯一最优解y＾。

同时由于是强凸的［利用式（A－15）］，也满足对任何y，有

将式（A－22）记作≜proxτf（x），即函数f的迫近算子。proxτf（x）是1－Lipschitz的单调算子。通过次梯度微分运算，有换言之，迫近算子proxτf（x）由极大单调算子τ∂f在x处的预解算子（I＋τ∂f）－1x给出。通常，算子T是极大单调的，当且仅当其预解算子（I＋τT）－1是良定义且单值的。预解算子是弱收缩的，或者称为“1／2－平均算子”［316］。

将预解算子与式（A－17）相结合，可以导出Moreau identity等式：

事实上，对任何极大单调算子T和T－1，式（A－25）都成立。通过式（A－25），如果知道如何计算proxτf（x），就能计算proxf∗／τ。有时也会出现在某个测度下计算迫近算子的情况，记作

A2.5　Fenchel－Rockafellar对偶（duality）

有时对偶性可以用于将一个凸问题转换成另一个更容易求解的问题。考虑如下问题：

在有限维情况下，如果存在一个点x存在于dom g内部，而且Kx在dom f的内部，则可以交换式（A－28）中min和sup的次序［317］，即

最后一个问题就是对偶问题。在上述假设下，式（A－29）至少存在一个解y∗。记拉格朗日函数（Lagrangian）为

如果x∗是原始问题式（A－27）的任何一个解，则（x∗，y∗）是原始—对偶问题的一个鞍点：对于任何（x，y）∈X×Y，有

由最优化条件，有

此时，原始—对偶隙（primal－dual gap）为

原始—对偶隙始终为非负值，当且仅当（x，y）为鞍点时原始—对偶隙为零。

最后，最优化条件式（A－32）可以通过一个极大单调算子T来表达：其中

A3　ADMM算法

A3.1　ADMM算法简介

ADMM（Alternating Direction Method of Multipliers）算法融合了原始—对偶算法的可分解性与乘子法优良的收敛性两重优势。ADMM可以用于求解以下形式的问题。

此处，变量x∈Rn，z∈Rm，矩阵A∈Rp×n，B∈Rp×m，c∈Rp；函数f和g为凸函数。变量x分裂为x和y两部分，从而使得目标函数也变成分别依赖x和y可分离的。式（A－35）的最优解为

基于乘子法，可以构造增广拉格朗日函数为

ADMM由以下迭代构成：(www.xing528.com)

式中，惩罚参数ρ＞0。

ADMM算法包含一个关于x的最小化步骤［式（A－38）］、一个关于z的最小化步骤［式（A－39）］，以及一个对偶变量更新［式（A－40）］。与乘子法相同，对偶变量更新所用步长等于增广拉格朗日乘子参数ρ。

式（A－35）的乘子法形式为

即增广拉格朗日函数是同时相对于两个原始变量联合优化的。而ADMM是按照交替方式轮流优化x和z的，因此可以看作Gauss－Seidel的变量更新方式。把关于x和z的最小化分成两步，恰好使f或g可分离时，将原始复杂的目标函数分解为简单的子问题。

A3.1.1　Scaled ADMM

ADMM可以写为另一种更方便的方式，即将增广拉格朗日中的线性和二次项结合起来并且缩放对偶变量。定义残差为r＝Ax＋Bz－c，则有

此处，u＝y／ρ称为比例化对偶变量（scaled dual variable）。使用比例对偶变量，ADMM可以写作

定义第k步的残差为rk＝Axk＋Bzk－c，则可以看出

是残差项的累积和。

比例形式的ADMM［式（A－41）～式（A－43）］相对于非比例形式的ADMM［式（A－38）～式（A－40）］更为简洁方便，因此后续介绍中将使用比例形式的ADMM。

A3.2　ADMM收敛性

ADMM在非常一般的条件下，即可保证该收敛性：

条件1：扩展实值函数f：Rn→R∪｛＋｝和g：Rn→R∪｛＋｝是闭合、特展、凸函数。

该条件意味着x－更新［式（A－38）］和z－更新［式（A－39）］是有解的，即存在x和z（不一定唯一）可以最小化增广拉格朗日函数。条件1重要的一点是，允许函数f和g是不可微的，并且允许函数取值为＋。这使得可以取f为一个闭合非空凸集C的示性函数，即f（x）＝0如果x∈C，否则f（x）＝＋。在这种情况下，最小化x步骤将涉及求解在凸集上的约束二次规划问题，即f的实际定义域。

条件2：增广拉格朗日函数L0（x，z，y）有一个鞍点。

也就是说，存在（x∗，z∗，y∗），并不一定唯一，使

L0（x∗，z∗，y）≤L0（x∗，z∗，y∗）≤L0（x，z，y∗）

对所有x、y、z都成立。

由条件1可知，对任何鞍点（x∗，z∗，y∗），L0（x∗，z∗，y∗）都是有限的。这意味着（x∗，z∗）是原始问题式（A－35）的一个解，所以Ax∗＋Bz∗＝c，并且f（x∗）＜＋，g（z∗）＜＋。这也意味着这组解是对偶最优的，并且原始与对偶问题的最优值相等，即达到强对偶性。注意，此处并未对A、B和c做任何要求（当然条件2有一些隐含约束），所以并不要求A、B是满秩的。

A3.2.1　ADMM收敛性讨论

满足条件1和2，ADMM的迭代将有如下特性。

•残差收敛，当k→＋时，rk→0，即迭代逼近可行解；

•目标函数值收敛，当k→＋时，f（xk）＋g（zk）→p∗，即目标函数值逼近最优；

•对偶变量收敛，当k→＋时，yk→y∗，即对偶变量逼近最优对偶点y∗。

该收敛性证明过程可参见文献［241］。

A3.2.2　实际问题中的收敛性

有实验表明，在要求达到极高的精度（如10－10）时，ADMM收敛速度会比较慢。但是对于大多数实际应用问题，中等精度即可满足需要，此时ADMM通常只需要几十次迭代即可达到要求。通常，对于大规模问题，如成像中的逆问题，当一个像素值的估计达到极高精度，并不会对最终重建结果带来实质性的区别。因此，ADMM对众多仅需中等精度的问题来说是一个很好的选择。

A3.3　最优条件和终止准则

ADMM问题式（A－35）的充分必要最优条件是原始可行，即

以及对偶可行，即

此处，如果f和g是可微函数，则可以用梯度∇f、∇g代替次梯度∂f、∂g，同时用＝代替∈。可见，ADMM的对偶可行条件式（A－45）、式（A－46）与式（A－34）中的极大单调算子最优化条件是一致的。

由于zk＋1要最小化Lρ（xk＋1，z，yk），故有

这意味着zk＋1和yk＋1总是满足条件式（A－46）。这样，最优化条件简化为满足式（A－44）和式（A－45）。下面继续验证这两个条件。

由于xk＋1要最小化Lρ（x，zk，yk），故有

这等价于：

也就是说，可以定义：

将其看作对偶可行条件式（A－45）在第k＋1次迭代时的对偶残差。同时，称rk＋1＝Axk＋1＋Bzk＋1－c为第k＋1次迭代时的原始残差。

总体来说，ADMM问题的最优化条件由式（A－44）～（A－46）构成。最后的式（A－46）对于每次迭代（xk＋1，zk＋1，yk＋1）都成立；而式（A－44）和（A－45）的残差项分别为原始残差rk＋1与对偶残差sk＋1。这两项残差在ADMM迭代过程中收敛为零。

A3.3.1　终止准则

可以将最优化条件的残差项rk和sk用作当前的迭代点的优化程度的上界，即衡量目标函数值与最优值之间的距离f（xk）＋g（zk）－p∗。可以证明［241］：

这表明，当残差项很小时，目标函数的次优性也一定很小。由于不知道真实的x∗，不等式（A－47）并不能被直接用于迭代的终止准则。但是如果能够估计或猜测‖xk－x∗‖≤d，则可以得到

上面不等式中的第2或第3项都可以用于衡量目标函数的次优性。这意味着可以设置一个合理的终止准则，使得原始和对偶残差都很小，即

此处εpri＞0和εdual＞0分别为原始可行条件式（A－44）、对偶可行条件式（A－45）的可行容差。这些容差可以通过使用绝对和相对准则来设置，例如：

式中，εabs＞0是绝对容差；εrel＞0是相对容差。而因子和是表征l2范数，分别定义在Rp和Rn空间上。根据应用需要，一个合理的相对终止容差可以是εrel＝10－3或10－4，而绝对容差εabs则根据典型变量值的量级设置。

A3.4　ADMM拓展

目前存在很多标准ADMM算法的变体，它们可以显著提升算法的收敛性能，下面予以简单介绍。

A3.4.1　变惩罚参数

与使用固定的惩罚参数ρ不同，可以在每个ADMM迭代步骤中使用不同的惩罚参数ρk以提高收敛速度，并减少对惩罚参数初始设置的依赖。尽管很难证明当ρ变化时的ADMM的收敛性，固定ρ的收敛证明仍然能保证当经过有限次变化之后ρ保持固定值时，该种方法仍然具有收敛性保证。

一个简单有效的方案［318，319］是：

式中，μ＞1，τincr＞1，τdecr＞1是参数，需根据具体问题设置。典型选择是μ＝10，τincr＝τdecr＝2。这种处理方式的想法是让原始和对偶残差范数在收敛到0的过程中彼此保持在μ倍的比例之内。

ADMM更新方程表明，大的ρ值会对原始可行偏离给予大的惩罚，从而获得小的原始残差。反过来，对偶残差sk的定义又表明，小的ρ值会减小对偶残差，但是会降低对原始可行偏离的惩罚从而引起大的原始残差。式（A－49）中的调整方案会在原始残差相对于主残差过大时将ρ增大τincr倍，而当原始残差相对于主残差过小时将ρ降低为原来的1／τdecr。这个方案也可以进一步考虑调整εpri和εdual的相对大小来获得提升。

需要指出的是，当在比例形式的ADMM中调整惩罚参数时，比例对偶变量uk＝yk／ρ也必须相应地更新，即如果ρ降低了一半，uk应该增加为原来的2倍。

A3.4.2　更一般的增广项

除上述以外，还可以对每个等式约束变量使用不同的惩罚参数，比如，将二次项（ρ／2）‖r‖22替换为（1／2）rTPr，P是一个对称正定矩阵。当P保持恒定时，这种方式可以看作将标准ADMM的等式约束r＝0替换成了Fr＝0，而FTF＝P（注意任何正定矩阵都有此分解形式）。

A3.4.3　过松弛

在z更新步骤式（A－39）和y更新步骤式（A－40）中，可以将Axk＋1替换为

αkAxk＋1－（1－αk）（Bzk－c）

式中，αk∈（0，2）称为松弛参数。当αk＞1时，该技术称为过松弛；当αk＜1时，称为欠松弛。文献［320，321］中的实验表明，当取过松弛αk∈［1.5，1.8］时可以提升性能。

A4　FISTA算法

FISTA（Fast Iterative Shrinkage－Thresholding Algorithm）［200］是一个在求解信号和图像重建问题中常用的简单高效算法，能够适用于涉及稠密矩阵的大规模线性问题求逆。FISTA在ISTA（Iterative Shrinkage－Thresholding Algorithm）算法［322，323］的基础上，将ISTA算法中每次迭代步骤仅依赖于前一次迭代结果xk－1，变为依赖于前两次结果｛xk－1，xk－2｝的一个巧妙的线性组合，从而将ISTA算法的最坏收敛速率O（1／k）提升到O（1／k2），而且通过实验验证了FISTA的收敛速度比ISTA算法提高了几个数量级。

A4.1　ISTA算法原理

ISTA可以用于求解如下形式的凸优化问题：

其中：

•f：Rn→R是C1，1型平滑凸函数，即连续可微，具有Lipschitz连续梯度：

‖·‖为向量2－范数；L（f）＞0为梯度∇f的Lipschitz常量；

•g：Rn→R是连续凸函数但可能不平滑；

•问题式（A－50）可解，即解集非空，X∗：＝argmin F≠∅，而且对于x∗∈X∗，取函数值F∗＝F（x∗）。

在式（A－50）的模型中，如果令g（x）≡0，则问题变为普通的无约束平滑凸优化问题，经典梯度下降法即可获得最优解。如果f（x）＝‖Ax－b‖2，g（x）≡‖x‖1，则问题变为经典的LASSO（Least Absolute Shrinkage and Selection Operator）问题［324］，此种情况下，梯度∇f的最小Lipschitz常量L（f）＝2λmax（ATA）。下面的讨论给出了两种情况之间的联系。

对于无约束连续可微函数f：Rn→R的最小化问题：

其梯度下降解法将得到点序列｛xk｝：

式中，tk为合适的步长参数。该解法可以看作是f在xk－1处的线性逼近的迫近规整化（proxi⁃mal regularizaition）：

令f（x）＝‖Ax－b‖2，将此方法应用到非平滑的l1规整化问题上：

可得到迭代方案：

此处，Tλtk：Rn→Rn为一个收缩算子（或称软阈值算子［325］），定义为

式（A－55）和式（A－56）即ISTA算法的求解方法，是一种一阶优化方法，在梯度计算上仅涉及矩阵—向量乘法，计算复杂度较低。但是ISTA算法的一个严重缺陷是其仅具有次线性全局收敛率（sublinear global rate of convergence）：

F（xk）－F（x∗）＝O（1／k）

式中，k为迭代次数。FISAT算法将收敛速率提高到了O（1／k2），下面予以介绍。

A4.2　FISTA算法步骤

FISTA算法的简单高效在于，其每一次迭代仅涉及目标函数F（x）中的平滑项f（x）的一次梯度计算，而梯度计算的主要计算代价为两次涉及矩阵A和AT的矩阵—向量乘法，尤其是当A具有特殊结构，如为卷积或稀疏矩阵时，计算代价还可进一步降低。对于梯度下降所需步长，可以使用∇f的Lipschitz常量L（f）＝2λmax（ATA），也可以使用变步长策略以避免大规模矩阵的谱半径估计。下面分别介绍这两种FISTA算法的实现流程。

使用固定步长的FISTA算法流程如下：

输入：L＝L（f），L（f）为∇f的Lipschitz常量。

Step 0：初始化y1＝x0∈Rn，t1＝1。

Step k：对于k≥1，依次计算：

比较式（A－55）中的ISTA迭代与式（A－57）中的FISTA迭代，可以看出，ISTA的软阈值算子只作用在前一次迭代结果xk－1上，而FISTA则是作用在前两次迭代结果｛xk－1，xk－2｝的线性组合yk上。正是式（A－58）中的迭代参数tk的递归设计策略，使得FISTA获得了O（1／k2）的收敛率。

由于L（f）＝2λmax（ATA）需要计算ATA的最大特征值，这在A的规模较大时比较困难，此时可以使用反向跟踪法来确定步长参数。

使用反向追踪步长的FISTA算法流程如下：

Step 0：取L0＞0，参数η＞1，初始化y1＝x0∈Rn，t1＝1。

Step k：对于k≥1，找到最小的非负整数ik，使得能够满足如下条件：

取Lk＝ηikLk－1，依次计算：

按照反向追踪法选择出来的步长Lk，满足如下条件。

式中，α＝β＝1对应固定步长设置，对应反向追踪步长设置。由于选择步长时满足条件由A4.3节中的引理1可知，。由于ik是最小非负整数满足，所以将小于等于L（f），即Lk≤L（f）。

A4.3　FISTA算法原理

对于任何L＞0，考虑F（x）≜f（x）＋g（x）在给定点y的二次逼近：

其具有唯一最小值点：

可见，pL（y）即式（A－22）中的迫近算子取为

下面首先给出三个引理，它们将用于FISTA的收敛性证明与步长选择。

引理1　令f：Rn→R是连续可微凸函数，具有Lipschitz连续梯度常量L（f），则对任意L≥L（f），有

引理2　对任何y∈Rn，记z＝pL（y），如果存在γ（y）∈g（z），则有

取式（A－60）关于x的最优化条件：

显然，当x取值z＝pL（y）时，式（A－63）成立。

引理3　令y∈Rn，L＞0，使得

则对任何x∈Rn，有

注意，式（A－64）将成为可行步长L选择的判断标准。根据引理1，只要L≥L（f），式（A－64）总是成立的。下面的定理给出了FISTA的非渐进收敛率。

定理FISTA收敛率　令｛xk｝，｛yk｝为FISTA算法生成的序列，则对任何k≥1，有

式中，α＝1对应常步长设置，而α＝η为反向追踪步长设置。

该定理表明，FISTA算法具有1／k2的收敛率。为了达到指定的精度ε，即迭代结果满足F（x～）－F∗≤ε，FISTA所需的迭代次数最多为