1引言

本综述关注在反问题背景下优化问题的有效求解。当数据集规模庞大或前向算子计算代价高昂时，这一问题尤为重要。为铺垫背景，我们旨在从间接含噪测量 \(v \in \mathcal{V}\) 中恢复一个对象 \(u^{\dagger} \in \mathcal{U}\)。我们假设 \(u\) 与 \(v\) 之间的关系可由前向模型近似

其中 \(K: \mathcal{U} \to \mathcal{V}\) 是 Hilbert 空间 \(\mathcal{U}\) 与 \(\mathcal{V}\) 之间的映射。许多任务可通过这一表达式来描述，例如图像处理中的去模糊、成像中的层析重建，或统计学中的回归，仅举数例。层析成像中的具体应用包括计算机断层扫描（CT）、正电子发射断层扫描（PET）和磁共振成像（MRI）。

在本综述中，我们主要考虑线性反问题，其中映射 \(K\) 是线性算子，并用 \(Ku\) 表示它对元素 \(u\) 的作用。尽管非线性引入了额外的数学和计算挑战（如非凸性），但本文介绍的许多技术可以 readily 推广到非线性问题。

在实际应用中，解 \(u\) 和数据 \(v\) 通常都是离散化的，我们寻求从 \(\mathbb{R}^s\) 中的测量恢复 \(\mathbb{R}^d\) 中的一个元素。复数可以通过将 \(\mathbb{C}\) 与 \(\mathbb{R}^2\) 等同而方便地纳入。分量 \(v_i\) 代表第 \(i\) 个测量值（例如来自数字探测器），而 \(u_i\) 对应于解的第 \(i\) 个参数（例如某一图像体素的值）。

许多现代应用的一个共同点是：\(d\) 和 \(s\) 都很大，或者前向算子 \(K\) 的计算代价高昂。因此，计算反问题 (1) 的近似解在计算上是具有挑战性的。这两个问题通常高度相关，例如矩阵-向量乘积的代价随输入和输出维度线性增长。为了说明数据的高维度，考虑 PET-MR 扫描仪 Siemens Biograph mMR，其未压缩数据（所谓"span 1"）有超过 3.5 亿个元素，而图像通常建模为约 2000 万个体素。需要强调的是，"大"和"昂贵"是相对的，必须与特定应用中可用资源相比较。

与其直接求解 (1)，不如通过优化问题来近似解：

其中 \(D: \mathcal{V} \times \mathcal{V} \to \mathbb{R}_{\infty} := \mathbb{R} \cup \{\infty\}\) 衡量 \(Ku\) 近似 \(v\) 的程度，通常被称为数据保真项（data-fidelity）或数据拟合项（data-fit term）。数据保真项 \(D\) 通常是非负函数，在参数相同时最小化，即 \(\inf_{v'} D(v',v) = D(v,v) = 0\)。允许 \(D\) 取值为 \(\infty\) 有时需要将 \(D\) 的定义扩展到向量空间（例如见下面的 Kullback-Leibler 散度）或建模约束。

另一种观点是给定 \(Ku\) 建模并最大化观察到 \(v\) 的似然。数据保真项可以定义为负对数似然，所寻求的估计量即为最大化似然的那个。

最常用的数据保真项是可分的。在有限维数据设定中，\(\mathcal{V} = \mathbb{R}^s\)，这意味着它们可以写成

其中函数 \(d: \mathbb{R} \times \mathbb{R} \to \mathbb{R}_{\infty}\) 可以使用数据的统计性质或测量噪声来确定。例如，使用高斯分布、拉普拉斯分布和泊松分布分别导出最小二乘、绝对值差异和 Kullback-Leibler 散度，定义如下：

• 最小二乘：\(d(v'_i, v_i) = \dfrac{1}{2\sigma_i^2}(v'_i - v_i)^2\)
• 绝对值差异：\(d(v'_i, v_i) = |v'_i - v_i|\)
• Kullback-Leibler 散度： \[ d(v'_i, v_i) = \KL(v'_i \mid v_i) = \begin{cases} v'_i \log\dfrac{v'_i}{v_i} - v'_i + v_i, & \text{若 } v'_i > 0 \text{ 且 } v_i > 0, \\ v_i, & \text{若 } v'_i = 0 \text{ 且 } v_i \geq 0, \\ \infty, & \text{其他情况}. \end{cases} \]

虽然高斯假设应用最广泛，但泊松分布数据也获得了大量关注。

1.1变分正则化

反问题 (1) 常常是不适定的（ill-posed），在这个意义上它可能没有解，有无穷多解，或者即使解存在且唯一，它对测量中的噪声也高度敏感。求解优化问题 (2)——它最小化重建对测量的保真度——可以帮助克服解的不存在性。然而，它继承了其他问题，如解的非唯一性和噪声放大。例如，如果 \(K\) 是病态的，观测中的微小变化（例如由于噪声）可能导致重建解的巨大变化。因此，仅最小化数据保真度将导致对噪声过拟合的重建。解决这个问题的常用方法是向优化目标中添加正则化项 \(R: \mathcal{U} \to \mathbb{R}_{\infty}\) 并改为求解

其中 \(\lambda > 0\) 被称为正则化参数。这一概念被称为变分正则化，在过去几十年中得到了广泛研究。多年来已使用了许多正则化项 \(R\)。我们不一一全面介绍，而是讨论具有与反问题相关性质的代表性正则化项。

若 \(\mathcal{U} = \mathbb{R}^d\)，最简单的正则化项是平方 2-范数：\(R(u) = \frac{1}{2}\|u\|_2^2 = \frac{1}{2}\sum_i u_i^2\)。这通常被称为 Tikhonov 或 Miller 正则化，可用于稳定重建并使模型解唯一。或者，使用 1-范数 \(\|u\|_1 = \sum_i |u_i|\) 可以促进解中的一种特定性质：稀疏性，即解中只有少数系数非零。

用 1-范数和平方 2-范数进行正则化不能很好地捕捉成像应用中通常遇到的期望解的结构。它们促进具有较小 entries 的重建，这对图像来说是不可取的。此外，它没有考虑图像中的局部邻域或平滑性特征。这可以通过不直接惩罚 \(u\) 而是惩罚其梯度来实现：\(R(u) = \frac{1}{2}\|\nabla u\|_2^2\)。然而，严格来说，这只对光滑函数 \(u\) 有良好定义。另一种观点是将 \(\nabla\) 解释为计算离散化向量 \(u\) 的（空间）有限差分的线性算子。

使用 \(\nabla u\) 的平方 2-范数会促进相邻像素之间的相似性，并可能导致物体边缘的模糊。为了保留边缘，通常用 1-范数替代梯度的平方 2-范数，即

这定义了全变分（Total Variation, TV），如果期望解可以被很好地描述为分段常数函数——自然图像的一个常见模型——那么它是一个特别好的模型。这个公式对光滑函数 \(u\) 有良好定义，但可以推广到具有较低正则性的其他函数。

已知当图像具有分段线性结构时，TV 会引入阶梯伪影。为了规避这个问题，已经提出了 TV 的几种高阶变体。也许最受欢迎的选择是广义全变分（Total Generalised Variation, TGV），定义为

其中 \(\beta > 0\)。这里 \(w\) 是与 \(u\) 的分段线性分量相关的向量场，而 \(\nabla u - w\) 是 \(u\) 的分段常数分量的向量场。TV 和 TGV 的比较示于图 1。

TV 的一种非光滑替代是所谓的 Huber TV，它将 1-范数替换为光滑近似。在离散化设定中可以写成 \[ R(u) = \sum_i h_{\gamma}(\|(\nabla u)_i\|_2), \] 其中求和遍历梯度场 \(\nabla u\) 的所有分量，且 \[ h_{\gamma}(t) = \begin{cases} t & \text{若 } |t| > \gamma, \\ \dfrac{1}{2\gamma}t^2 + \dfrac{\gamma}{2} & \text{若 } |t| \leq \gamma. \end{cases} \]

Huber TV 是可微的，具有 \(\|\nabla\|_2/\gamma\)-Lipschitz 连续梯度，也可以帮助克服 TV 引入的阶梯伪影。

约束提供了另一种正则化来源。与之前讨论的模型（鼓励解具有某些性质）不同，约束可用于强制期望的性质。成像中最流行的约束是 box 约束，我们只寻求 \(d\) 维 box 中的解 \[ C = \{ x \in \mathbb{R}^d : \ell_i \leq x_i \leq r_i, \, i = 1,2,\ldots,d \}. \] 这包括非负性约束，其中取 \(\ell_i = 0\) 和 \(r_i = \infty\)。

在本综述中，我们关注优化视角，即关注高效计算 MAP 估计量。虽然随机近似的其他统计方面对某些读者也可能感兴趣，例如最小均方误差估计，但这超出了本优化聚焦综述的范围。

1.2优化模板

如前所述，在变分正则化的背景下，一旦选择了数据拟合项、正则化项和正则化参数，反问题的解就可以通过优化 (4) 来计算。上述所有优化问题都是如下优化模板的特例

其各分量满足以下假设，我们将在整篇文档中默认使用：

注意，我们不想也不能计算 \(f \circ A\) 的近端算子，而是希望将函数 \(f\) 与算子 \(A\) 解耦，并分别处理它们。实现这一点的技术在 第 2 节中讨论。

从给定的变分正则化最小化问题到优化模板的转换不是唯一的，如下面的示例所示。

虽然上述示例主要涉及 TV 正则化，但它们可以 readily 推广到其他类型的正则化，例如 Huber TV 或小波正则化（通过将梯度算子 \(\nabla\) 替换为小波变换 \(W\)）。

1.3延伸阅读

高级建模。本综述聚焦于变分正则化最常用的数学模型。然而，这只是冰山一角。例如，变分正则化的一个缺点是正则化项引入的偏差，特别是在数据信息不足时。偏差可以通过 [BBRS17] 的两步技术或迭代正则化 [BB18] 来减少。另一个有趣的方向是多测量的耦合，可以对应于不同模态或能量通道，或者将问题视为随时间变化的动态反问题。此类前向模型需要定制的数学建模，可以在变分正则化的背景下完成。

机器学习。近年来，出现了针对反问题的学习正则化项，例如字典学习、fields-of-expert 模型、梯度步去噪器、regularisation-by-denoising、deep image prior with TV、生成正则化项、input-convex 神经网络（ICNN）或网络 Tikhonov 正则化。其中，重点在于学习信息丰富且凸的正则化项，它们利用特殊参数化的神经网络，如 ICNN 和凸 ridge 正则化项。

1.4其他综述

据我们所知，这是第一篇聚焦于反问题随机优化算法的综述。然而，已有其他关于反问题高效计算的综述，其中一些也关注优化方法。[CFKS13] 和 [CKN15] 讨论了大规模反问题的高效求解，聚焦于来自数值线性代数的确定性算法。[BSS16, CP16] 和 [Val16] 关注反问题和成像的高效一阶优化算法。虽然我们也详细讨论一阶算法，但我们的特别关注点是克服问题大规模特性的随机方法。[BBM23] 综述了用于图像分类的神经网络训练的随机方法，但不考虑反问题。[XN22] 讨论了凸优化与机器学习和医学成像 AI 的联系。他们确实触及了随机优化，但没有深入对反问题感兴趣的细节或 specifics。

2变分正则化的确定性算法

基于 第 1.2 节的讨论，我们回忆优化模板 (7) 的一般形式 \[ \min_{x \in \mathcal{X}} \{ \Phi(x) = f(Ax) + g(x) + h(x) \}, \] 在假设 (A.1)–(A.3) 下。为这类优化问题开发高效求解器一直是一个活跃且富有成效的研究领域，特别是自新千年以来。在现有方法论中，一阶方法——仅利用目标函数的梯度信息——一直是最受欢迎的选择。此外，正如我们将在 第 3 节中看到的，一阶方法是大多数随机优化算法的骨干。

在本节中，我们呈现一阶方法的概述，从基础原理到最新算法进展。在 第 2.1 节中，我们讨论求解一般变分优化问题的困难，以引入一阶方法的基本构件和设计原理。在 第 2.2 节中，我们呈现用于最小化 (7) 的不同变体的代表性一阶算法。我们在 第 2.3 节中以进一步主题的讨论作为结尾，包括高阶方案、非欧几里得方法、非凸优化和非线性问题。

我们强调，一阶优化是一个丰富且快速发展的研究领域。因此，我们的介绍并不全面。我们参考专著 [BC11, Bec17] 和综述 [CKCH23, BB09, BB18, BBES21] 以获取更多细节。

2.1一阶算法的基本构件

求解 (7) 的挑战来自以下方面：非光滑性、函数-算子复合以及多项之和。一阶优化算法的发展反映了反问题、计算机科学和工程中的挑战日益复杂，以及对计算效率的日益增长的需求。

一阶优化算法依赖梯度信息来指导最小化目标函数的迭代过程。这一方法家族在过去几十年中显著演化，始于简单但基础的梯度下降算法，扩展到更先进的技术，如动量加速、近端算法和现代算子分裂技术。其演化中的每个阶段都是由解决早期方法的局限性、适应新问题结构和提高计算效率的需求所驱动的。

在一阶优化算法的历史发展中，分而治之（divide and conquer）的思想一直是基石。它作为解决上述挑战的指导原则，包括将原问题分解为更简单的子问题，分别处理每个子问题，并将所有组件组装成一个可证明且可实现的数值方案。

2.1.1基本要素

经过数十年的发展，分解优化问题（如 (7)）的技术已经变得成熟。一般来说，可以应用以下原则：直接分解目标、基于最优解的最优性条件分解、或就优化变量进行分解。

无论采用何种分解方法，"征服"步骤本质上归结为管理 \(h\)、\(g\) 和 \(f \circ A\)。在本节中，基于 第 1.2 节 中概述的项的性质 (A.1)–(A.3)，我们讨论分别处理这些项的基本方法。

梯度下降法（Gradient Descent, GD）。追溯到至少 19 世纪中叶，GD 是无约束优化最古老的迭代方法之一。当最小化一个凸的、光滑的、可微的目标 \(h\) 时，它迭代地沿当前迭代点的负梯度方向执行一步

其中 \(\tau_k > 0\) 是步长。正如其名称所示，GD 是一种下降方法，在步长 \(\tau_k\) 选择适当时每次迭代减少函数值 \(h(x^{(k+1)}) \leq h(x^{(k)})\)。当函数 \(h\) 是 \(L\)-光滑时，即 \(\nabla h\) 是 \(L\)-Lipschitz 连续的，GD 在适当的常数步长下达到 \(O(1/k)\) 的目标函数值收敛速率 [Bec17]。

近端点算法（Proximal Point Algorithm, PPA）。当最小化非光滑函数时，仍然可以通过用次梯度替换梯度来应用 GD（定义见 (A.3)）。然而，所得次梯度方法不再是下降算法，且不一定在常数步长下收敛。相反，通常需要递减的步长来保证收敛。此外，次梯度下降只获得 \(O(1/\sqrt{k})\) 的收敛速率，比 GD 差。一种更流行的最小化非光滑目标的方法是使用 \(g\) 的近端算子（定义见 (A.7)），这产生了 PPA [Roc76]。

对于适当的、闭的、凸的函数 \(g\)，相应的 PPA（或近端最小化算法）取形式

其中 \(\tau_k > 0\) 是步长。在步长不太快趋于零时保证收敛到最小化点，例如 \(\sum_k \tau_k = \infty\)。PPA 也是一种下降方法 [Bec17]，在常数步长下它承认 \(O(1/k)\) 的目标函数值收敛速率。

注意，每次迭代需要求解一个优化问题来计算近端算子。然而，许多广泛使用的（正则化）函数的近端算子是已知闭式的 [CP11b, CCCP]。

加速。加速是任何有效一阶方法的基本组成部分。最成功的加速技术是惯性加速（inertial acceleration）或动量（momentum），追溯到重球法 [Pol64]，它可以显著加速 GD，但缺乏全局加速保证。Nesterov 加速梯度（NAG）方法是第一个惯性加速方法，它可证明地将光滑凸目标的 GD 目标函数收敛速率从 \(O(1/k)\) 提高到 \(O(1/k^2)\) [Nes83]。快速迭代收缩阈值算法（FISTA）[BT09] 将此加速推广到非光滑凸情形。我们参考 [ABC23] 以获取过去十年惯性加速发展的更全面的处理。其他加速技术包括预处理 [PC11, LXY21] 和向量外推 [And65, SdB16]，仅举数例。我们参考专著 [dST21] 以获取更多关于加速的细节。

函数-算子分解。复合项 \(f \circ A\) 的近端算子通常难以计算，我们参考附录以了解一些例外。这通常通过两种流行方法之一来解决：对偶性或多重乘子法。前者使用 Fenchel-Moreau 定理，它指出一个适当的、闭的凸函数的二次共轭（定义见 (A.4) 和 (A.5)）等于函数本身，这给出

这允许将底层优化问题重构为鞍点问题，然后可以用原始-对偶算法 PDHG（见算法 3）来求解。另一种选择是通过辅助变量将 \(f \circ A\) 分解为约束问题 \[ \min_{x \in \mathcal{X}, z \in \mathcal{Z}} \{ f(z), \text{ s.t. } z = Ax \}. \] 通过将约束加回目标函数并通过 Lagrange 乘子，该问题可以通过对偶上升、多重乘子法 [BPC+11] 或交替方向乘子法（ADMM）（见算法 4）来求解。PDHG 和 ADMM 有着深刻的联系，在某些条件下彼此等价，见 [CP11a, CP16, OV20] 的讨论。

坐标下降法（Coordinate Descent, CD）。不是在每次迭代中更新优化变量 \(x\) 的所有 entries，CD 方法只更新一个或一小部分坐标

这里 \(\mathcal{B}_k\) 是 \(x\) 坐标的子集，\(e_i\) 是第 \(i\) 个标准基向量，\(\nabla_i h(x^{(k)}) = \partial h(x^{(k)})/\partial x_i\) 表示沿第 \(i\) 个坐标的方向导数。当 \(\mathcal{B}_k\) 是单元素集时，(11) 是标准 CD 方案；否则它被称为块 CD。CD 方法通常用于具有可分结构或"坐标友好"的问题 [STXY16]。坐标更新顺序的流行确定性选择包括循环和固定排列方案。随机化顺序在 第 3.3 节 中讨论。CD 在优化中的引入归功于 [BT15]。收敛结果最初针对可微目标 [LT92] 展示，后来扩展到非可微目标 [Tse01]。我们参考 [Nes12, BT13, RT14] 以获取 CD 算法的理论进展。

2.1.2设计一阶算法

有了上述基本要素，设计一阶方法包括以下步骤：

• 问题分解：根据各函数的性质将原优化问题分解为更简单的子问题。
• 子问题处理与组装：为每个子问题分配一种处理方法（例如使用表 1 作为指南）。然后将所得组件组装成迭代方案。原则上，最优性条件可以作为算法组装的参考。
• 精化：迭代方案可能需要进一步改进以解决现有问题。首先，其一个或多个组件可能需要求解额外的优化问题，例如计算近端算子。其次，理论收敛保证可能缺乏或过于限制，因此需要修改方法。第三，理论或实际收敛可能较慢，因此需要加速或其他技术来提高性能。

需要强调的是，上述步骤不一定唯一。在下文中，我们呈现可应用于求解 (7) 的实例的经典一阶算法。我们从较简单的情况开始，其中 (7) 中只有两项非平凡，然后转向更一般的三项情形。

2.2高级一阶算法

在本节中，我们呈现用于最小化两项和三项优化问题的代表性一阶算法。

2.2.1两项和的方法

函数两项和的最小化在文献中得到了广泛研究，有几种流行的方法。此外，两项问题的方法通常用作更高级方法的构建模块。

近端梯度下降法（PGD）。当 \(f = 0\) 时，(7) 简化为两项和

其中 \(h\) 是光滑的，\(g\) 是近端友好的。求解 (12) 最流行的方案是 PGD（见算法 1），也被称为前向-后向分裂方法 [CW05]。PGD 对光滑项 \(h\) 执行 GD 步，然后计算 \(g\) 的近端算子。

PGD 是 GD 到非光滑情形的简单扩展，它共享相似的收敛性质：在 \(h\) 是 \(L\)-光滑且步长满足 \(\tau < 2/L\) 的条件下，PGD 以 \(O(1/k)\) 速率在目标函数中收敛 [BT09]。

加速 PGD。尽管 PGD 简单，但在实践中受到其缓慢的 \(O(1/k)\) 收敛速率的阻碍。快速迭代收缩阈值算法（FISTA）[BT09] 将 NAG 的加速方法论扩展到非光滑凸情形。

在 \(h\) 是凸的且 \(L\)-光滑且 \(\tau \leq 1/L\) 的条件下，FISTA 将目标函数值中的收敛速率提高到 \(O(1/k^2)\)。可以证明 FISTA 在修改的步长规则 \(a_k = (k-1)/(k+d), d > 3\) [CD15] 下收敛到全局最小化点。此外，作者还表明目标函数收敛速度是 \(o(1/k^2)\)。

尽管具有理论优势，但在实践中 FISTA 不一定比 PGD 快，特别是在寻求高精度解时。原因是 FISTA 不是下降方法，表现出（周期性）振荡行为，通常在问题（局部）强凸时。为了规避这个问题，重启方案可以显著改善 FISTA 的性能 [OC15, SBC16]。

原始-对偶分裂方法。当 \(h = 0\) 时，(7) 简化为

其中 \(f\) 和 \(g\) 都是非光滑但近端友好的。如前所述，\(f \circ A\) 的近端算子通常没有闭式表达，需要替代方法来处理 \(f \circ A\)。

使用对偶性，将 (10) 代入 (13) 得到鞍点问题

Chambolle 和 Pock [CP11a] 提出通过一种在原始变量中进行下降、在对偶变量中进行上升的算法来求解 (14)，应用于原始变量的外推，详见算法 3。

该算法被称为原始-对偶混合梯度（Primal-Dual Hybrid Gradient, PDHG）或 Chambolle-Pock 算法。PDHG 可以等价地写成 \(\mathcal{X} \times \mathcal{Y}\) 上的 PPA [HY12]。因此，可以获得 \((x^{(k)}, y^{(k)})\) 收敛到鞍点 \((x^*, y^*)\)，其中 \(x^*\) 是 (12) 的解。最著名的 PDHG 收敛结果是：在 \(\sigma\tau\|A\|^2 < 1\) [CP11a] 的条件下，原始-对偶间隙以 \(O(1/k)\) 收敛。在最近的工作 [MCJH23] 中，步长条件被改进为 \(\sigma\tau\|A\|^2 < 4/3\)。我们还参考 [GLY+13, GLY15] 以获取替代的自适应步长方案。

ADMM（交替方向乘子法）。通过多重乘子法处理 \(f \circ A\) 将 (13) 转化为约束问题

目标 \(\Phi\) 在 \(x\) 和 \(z\) 中是可分的，约束强制它们的耦合。关联的增广 Lagrangian 为

其中 \(\tau > 0\) 且 \(y\) 是 Lagrange 乘子。为了找到 \(\mathcal{L}(x,z;y)\) 的鞍点，ADMM（算法 4）一次更新 \(z\)、\(x\) 和 \(y\) 一个。

我们参考 [GM75, GM76] 以获取 ADMM 的早期发展和 [BPC+11] 以获取更现代的综述。ADMM 等价于将 Douglas-Rachford 分裂方法应用于 (13) 的对偶，见 [EB92, CP07]。

注意，辅助变量的更新可以写成 \(z^{(k+1)} = \prox_{f/\tau}(Ax^{(k+1)} - z^{(k)} + \tfrac{1}{\tau}y^{(k)})\)，因此是唯一确定的。另一方面，\(\mathcal{L}(x, z^{(k)}, y^{(k)})\) 的最小化点，即迭代 \(x^{(k+1)}\)，不需要唯一，除非 \(A\) 具有满列秩。

ADMM 的收敛分析通常比较复杂。然而，当 \(f\) 和 \(g\) 是闭的、凸的，且 \(\mathcal{L}\) 允许鞍点时，可以显示目标函数和约束中的收敛 [BPC+11]。原始序列 \(x^{(k)}\) 和 \(z^{(k)}\) 的收敛更难以建立，但当 \(A\) 具有满列秩时，可以显示 \((x^{(k)}, z^{(k)})\) 的收敛 [EB92]。

2.2.2三项和的方法

在上一节中，我们讨论了目标分解为两项之和的标准方法。那里介绍的大多数方法可以扩展到三项优化模板 (7)。然而，有时需要更大的灵活性。

Condat-Vũ 算法。Condat [Con13] 和 Vũ [Vũ13] 独立开发了一种用于 (7) 的原始-对偶前向-后向算法，作为包含光滑项梯度步的 PDHG 的推广，见算法 5。

Condat-Vũ 算法的迭代在 \(\tau(\sigma\|A\|^2 + L/2) < 1\) 的条件下收敛，其中 \(L\) 是 \(\nabla h\) 的 Lipschitz 常数 [Con13, Vũ13]。注意，每一步需要一次 \(A\)、\(A^*\)、\(\prox_{\tau g}\)、\(\prox_{\sigma f^*}\) 的求值和梯度 \(\nabla h\)。

PD3O。原始-对偶三算子算法（PD3O）[Yan18] 是 PDHG 和其他两项及三项算法的另一种推广。PD3O 与 Condat-Vũ 的区别在于 \(y^{(k+1)}\) 更新中的最后一项，其中 PD3O 添加了梯度动量项 \(\tau(\nabla h(x^{(k)}) - \nabla h(x^{(k+1)}))\)。

PD3O 比 Condat-Vũ 需要多一次梯度 \(\nabla h\) 的求值，但这些代价可以通过存储 \(\nabla h(x^{(k+1)})\) 并在下一次迭代中使用来摊销。已证明它在 \(\sigma\tau\|A\|^2 < 1\) 和 \(\tau L < 2\) [Yan18] 的条件下收敛。重要的是，收敛条件将来自 PDHG 和 GD 的条件解耦。这与 Condat-Vũ 算法形成对比，后者的收敛标准将参数 \(\tau\) 和 \(\sigma\) 与 \(\|A\|^2\) 和 \(L\) 耦合。

进一步的分裂方法。前面各节介绍的方法只是求解优化模板 (7) 的各种实例的一阶方案中的一小部分，列出所有现有算法超出了本文的范围。两项和的其他重要算法包括 Douglas-Rachford 分裂算法 [LM79, CP07] 和 Loris-Verhoeven 算法 [LV11]（也称为近端交替预测校正算法）。三项和的算法包括原始-对偶不动点算法（PDFP）[CHZ16]、原始-对偶 Davis-Yin 算法（PDDY）[DY17] 和非对称前向-后向伴随分裂（AFBA）[LP17]。现有的两项和三项算法通常密切相关。例如，在特定参数和函数选择下，PD3O 是 Loris-Verhoeven、Davis-Yin、PDHG 和 Douglas-Rachford 算法的推广。我们参考 [JV23, CKCH23] 以获取这些算法之间的进一步联系。ADMM 也可以扩展为求解三项和问题，然而，如 [CHYY16] 所示，所得算法不一定收敛，因此需要修改。

多于三项的优化问题（有时称为多块问题）的数值方法也已在文献中研究。[RFP13] 中作者提出了一种所谓的"广义前向-后向分裂"算法，用于最小化光滑函数和多个近端友好项。有几种原始-对偶分裂方法用于求解多块问题，例如 [CP12, Vũ13, CV14]。一种广泛采用的策略是块更新，它与 CD 有着相同的精神。

2.3延伸阅读

在这里，我们讨论与所介绍方法相关的高级主题，包括：高阶方法、Bregman 方法、非线性问题、学习优化和增量算法。

二阶方法。一阶方法有几个缺点。例如，它们通常表现出缓慢的收敛速度，不是尺度不变的，且通常需要繁琐的参数调优。高阶方法如牛顿法通过使用高阶导数来估计目标函数，从而解决了一些问题 [NW06]。然而，在许多成像应用中，计算 Hessian 的代价高得令人望而却步。拟牛顿方法尝试仅使用梯度信息来近似 Hessian 以克服这一点。也许最著名的拟牛顿方法是 Broyden-Fletcher-Goldfarb-Shanno（BFGS）算法 [Bro70, Fle70, Gol70, Sha70]。BFGS 在适当的凸性和光滑性条件下达到超线性收敛，也已应用于非光滑问题 [LO13, YVGS08]。然而，其内存需求可能很高，这导致了有限内存 BFGS（L-BFGS）[LN89] 的发展，它已应用于多个反问题 [Vog00, TBE+18]。

牛顿方法的超线性收敛可以推广到半光滑函数类的非光滑优化问题 [QS93]。半光滑牛顿方法已应用于 Lasso 问题 [LST18]，并已扩展到 PGD 和 Douglas-Rachford 分裂类型算法 [XLWZ18, HTPW22, LST23]、原始-对偶算法 [WFO23]，仅举数例，应用包括 TV 最小化 [NQYH07] 和地震层析成像 [BU15]。

Bregman 近端方法。本综述的大部分内容限制在 Hilbert 空间（或欧几里得空间）设定。推广是一个活跃的研究领域，通常依赖 Bregman 距离 [Bre67]。设 \(\varphi: \mathcal{X} \to \mathbb{R}\) 是严格凸且连续可微的函数。与 \(\varphi\) 关联的 Bregman 距离定义为 \(D_\varphi(x, y) = \varphi(x) - \varphi(y) - \langle \nabla \varphi(y), x - y \rangle\)。Bregman 距离通常不对称，不满足三角不等式。然而，它是非负的，且当且仅当 \(x = y\) 时等于 0。L-光滑性的概念可以推广为关于相应 Bregman 距离的"相对光滑性" [BBT17, Teb18]。在此框架内可以设计各种 Bregman 一阶方法，如 Bregman 近端算法 [BBT17]、Bregman ADMM [WB14]、Bregman PGD [MO19, HRX21] 和 Bregman 原始-对偶分裂 [CC22, JV22]，仅举数例。

非线性反问题。虽然许多应用可以用线性前向模型很好地建模，但也有许多应用线性模型是不够的，例如电阻抗断层成像（EIT）、相位恢复、光声断层成像、全波形反演（在地球物理学中有应用）、一般光学成像和特定的扩散光学层析成像，仅举数例。从优化的角度来看，这是具有挑战性的，因为所得最小化问题通常是非凸的，其正则性性质是异质的，大多数算法只保证收敛到局部最优解。一些凸优化算法已扩展到非凸优化，并已成功应用于反问题 [FG98, CGOT00, Val14, Val21, CZX+21]。近年来，通过利用 Polyak-或 Kurdyka-Łojasiewicz 条件等性质，在非凸问题的优化算法理论基础上取得了很大进展 [ABRS10, BST14, KNS16, FEHK22]。确定哪些非线性反问题满足此条件仍然是一个开放的研究问题。

学习优化（Learning-to-optimise, L2O）。随着深度学习的出现，学习优化正在成为优化中一个充满希望的方向 [LM16, CCC+22]。这一工作线在其早期旨在加速机器学习模型的训练，主要聚焦于使用强化学习改进网络训练 [ADG+16]。然而，大多数此类方法缺乏收敛保证，这激发了可证明收敛的 L2O 方案的研究，旨在学习算法的部分，例如 PDHG 的步长 [BRA+20] 和镜像下降算法的镜像映射 [TMTS23, TMTS24a]，应用包括机器学习和反问题成像。在 [EFG24] 中，作者开发了一个贪婪 L2O 框架，保证在训练数据上收敛。L2O 方案的理论泛化保证最近得到了发展 [SFO24]，引领了该领域的新研究方向。

迭代正则化。到目前为止，我们讨论了通过变分正则化 (4) 用优化方法求解反问题。另一种选择是不直接包含正则化项，而是提前停止迭代。例如，Landweber 迭代将 GD 直接应用于数据保真度。在步长、迭代次数和初始化选择适当时，Landweber 方法可以证明收敛到 (1) 的最小范数解（如果存在）[Lan51]。选择好的停止准则是重要且依赖于任务的问题。Landweber 迭代对非线性前向算子 \(K\) 的分析见 [HNS95]。将优化方法转化为正则化方法通常被称为迭代正则化，已用于线性化 Bregman 迭代 [YOGD08, Yin10] 及其非凸变体 [BBES21]、对偶算法 [GRV18]、原始-对偶算法 [MMRV21]，并在 Banach 空间中分析 [SLS06]，仅举数例。

3变分正则化的随机算法

第 2 节 中讨论的一阶数值格式在过去几十年中一直是求解反问题的主要方法，由于它们的效率、简单性和相对较快的收敛速度。然而，由于近年来问题复杂性和数据规模的快速增长，与一阶方法相关的计算成本正变得高得令人望而却步。这是开发具有低每次迭代计算成本的算法的主要驱动力，这些算法通常本质上是随机的。在本节中，我们呈现对反问题特别相关的现有随机优化方法的概述。本节布局与 第 2 节 相似。在 第 3.1 节 中，我们讨论一阶优化算法随机化的基本要素，在 第 3.2 和 第 3.3 节 中，我们介绍几种代表性算法。我们在 第 3.4 节 中以进一步主题的讨论作为结尾。

3.1设计一阶随机算法

在设计随机优化算法时，我们的目标是降低每次迭代的计算成本，同时不过度降低性能，从而创建一种整体上更快、更高效的方法。为了讨论随机算法的设计原理，我们首先重写问题 (7)，揭示有限和结构

对所有 \(i\)，算子 \(A_i: \mathcal{X} \to \mathcal{Y}_i\) 是线性有界的，\(f_i: \mathcal{Y}_i \to \mathbb{R}_{\infty}\) 和 \(g_i: \mathcal{X} \to \mathbb{R}_{\infty}\) 是适当的、闭的、凸的、近端友好的，\(h_i\) 是凸的且 \(L_i\)-光滑的。此外，我们记 \(\mathcal{Y} = \mathcal{Y}_1 \times \cdots \times \mathcal{Y}_{\ell}\)。

类似于 第 2.1 节，我们首先讨论解决这些挑战和设计随机算法的基本方法。

3.1.1基本要素

随机梯度近似。随机梯度方法的基本思想是用随机估计量来近似光滑项的全梯度。典型示例是随机梯度下降（SGD）[RM51]，它在每一步通过单个函数 \(h_i\) 的梯度 \(\widetilde{\nabla} h(x^{(k)}) = n \nabla h_{i_k}(x^{(k)})\) 来近似全梯度 \(\nabla h(x^{(k)}) = \sum_{i=1}^{n} \nabla h_i(x^{(k)})\)，使用随机选择的索引 \(i_k\)，有放回均匀抽样。因此，每次迭代的计算工作量减少了 \(1/n\) 的因子，从而 GD 的一次迭代对应于执行 \(n\) 次 SGD 迭代。这被称为一个 epoch 或数据遍历。SGD 的缺点是随机过程的方差引入的误差不会沿迭代减少。多年来已经开发了大量算法来解决这个问题，它们与 SGD 一起将在 第 3.2 节 中更详细地讨论。注意，索引的其他抽样形式也已被研究，例如无放回均匀抽样 [Sha16] 和不同形式的重要性抽样 [ZZ15, GDE22]。

随机 PPA。随机 PPA 是非光滑问题的 PPA 的随机对应物，在每次迭代中只计算一个随机选择的近端算子，即 \(x^{(k+1)} = \prox_{\tau_k g_{i_k}}(x^{(k)})\) [Bia16]。类似于光滑设定，可以推导方差缩减变体 [TASV23]。

随机化对偶 CD。与确定性算法的情况一样，随机算法也可以应用于问题的对偶形式。即，\(\sum_{i=1}^{\ell} f_i(A_i x)\) 的对偶问题计算为 \(\min_y \sum_{i=1}^{\ell} f_i^*(y_i)\) 约束为 \(\sum_{i=1}^{\ell} A_i^* y_i = 0\) [SSZ13]。注意，原始问题中的有限和结构导致对偶问题在 \(y\) 中是可分的，因此自然地适合对偶 CD 策略 [RT14]。这一想法已在 (17) 的原始-对偶算法中得到利用，如 QUARTZ [QRZ15] 和随机原始-对偶混合梯度（SPDHG）[CERS18]。

3.1.2设计随机算法

一旦确定了确定性优化算法的计算昂贵组件，它们就可以用随机替代方案替换。然而，仍有几个微妙的任务。例如，将 (7) 分解为有限和形式 (17) 不是唯一的，选择 \(\ell, m, n\) 和划分 \(A\) 以优化算法性能通常不是直接的任务。

在本节的其余部分，我们介绍几种代表性随机优化算法并讨论相关主题。

3.2随机梯度方法

在本节中，我们详细回顾基于梯度随机近似的算法，这一概念可以追溯到 1950 年代。在他们开创性的工作 [RM51] 中，Robbins 和 Monro 提出了一种随机近似程序，用于寻找可表示为期望值的函数的根。一个重要的洞察是许多大规模现实世界应用中的数据表现出大量真实或近似冗余 [BCN18]。因此，在每次计算中使用所有可用数据是低效的。此外，GD 已被证明相当鲁棒，即即使下降方向中有小误差，也能保证向最小值收敛，只要误差在长期内减少。这具有重大的实际好处，允许高效求解大规模优化问题。

自 2010 年代以来，由于机器学习复兴，随机梯度方案再次受到关注，现有方法库已得到显著丰富。新的算法修改包括自适应步长、方差缩减方案等。

3.2.1朴素 SGD

当我们考虑单个近端友好函数而没有复合项时，即 \(\ell = 0, m = 1\)，模板 (17) 简化为

SGD [RM51] 是一种经典的随机优化方法，用于高效求解这类问题，见算法 7。在每次迭代 \(k\) 中，它沿随机抽样索引 \(i_k \in \{1, \ldots, n\}\) 的（负）方向 \(\nabla h_{i_k}(x)\) 使用预定义步长计划（学习率）\(\tau_k > 0\) 执行一步。

与其在每次迭代中使用单个函数 \(h_i\)，在机器学习中通常使用几个随机选择的函数，产生所谓的 mini-batch SGD。相比之下，对于反问题成像，通过划分前向算子 \(K\) 来预分批数据并从划分定义函数 \(h_i\) 更为常见。在这种背景下，批次通常被称为子集，因此批次数量和批次大小变为子集数量和子集大小。我们将在 第 4 节 中更详细地讨论子集选择。SGD 已用于图像重建，例如 PET [TAK+23]、CT [GB18]、光学层析成像 [CLL18] 和冷冻电镜 [SS20, PRFB17]。

SGD 中的下降方向是全梯度的无偏估计量，即 \(\mathbb{E}[n \nabla h_i(x)] = \nabla h(x)\)。这种随机近似允许我们将 SGD 重写成以下形式

其中 \(e^{(k)} = \widetilde{\nabla} h(x^{(k)}) - \nabla h(x^{(k)})\) 是近似（或扰动）误差。从扰动角度看，为了证明收敛，我们需要控制近似误差。如果误差 \(e^{(k)}\) 是确定性的，根据 [SRB11] 的研究，在梯度误差下 PGD 的收敛性，(20) 在以下条件下收敛 \[ \sum_{i=0}^{k} \tau_i \|e^{(i)}\| = o(\sqrt{k}). \]

SGD 收敛分析通常施加有界方差假设：\(\mathbb{E}[\|n \nabla h_i(x) - \nabla h(x)\|^2] \leq C\)，\(C > 0\)；或施加所谓的"期望光滑性"假设 [KR20]，它提供了有界方差的较弱替代。

虽然上述假设只确保误差 \(e^{(k)}\) 是可控制的，但要实现收敛需要递减的步长，因为 SGD 对于常数步长不收敛。典型的递减步长采取以下形式

例如，\(\tau_k = k^{-\delta}\)，\(\delta \in (1/2, 1]\) 满足上述条件。然而，递减步长制度最终会导致非常慢的收敛速度。为了规避收敛速度的减慢，通常有两种类型的修改：方差缩减方案和自适应步长。

3.2.2方差缩减方案

随机方差缩减算法旨在通过历史梯度信息来减少 SGD 中存在的方差。在传统统计文献中，这种量有时被称为控制变量（control variates）。这些方法将 SGD 的低每次迭代复杂度与常数步长（以及更快的收敛速率）相结合，以增加存储需求为代价。方差缩减算法的研究始于随机平均梯度（SAG）[RSB12]、改进的随机平均梯度（SAGA）[DBLJ14] 和随机方差缩减梯度（SVRG）[JZ13]，但此后已显著增长。我们参考 [DLS22, GSBR20] 以获取概述。

方差缩减算法已用于许多图像重建问题，如 PET [TAK+23, KTA+21]、CT [KW17]、光学层析成像 [MAP20]、地震层析成像 [HCHY23] 和盲去模糊 [DTL+21]。在下文中，我们描述 SAGA 和 SVRG，两种 SGD 的代表性方差缩减变体。

SAGA 和 SAG。为了减少方差，SAG 和 SAGA 维护一个存储每个子函数 \(h_i\) 最近计算梯度的表。整个表在每次迭代中用于近似全梯度，利用所有梯度 \(\nabla h_i\) 的信息，从而降低方差。然而，表仅在最近使用（即随机抽样）的函数索引位置更新，这减少了计算工作量。

当 \(h_i\) 是凸的时，在 \(\tau = 1/(3nL_{\max})\) [Def15] 的条件下，SAGA 的平均迭代以 \(O(1/k)\) 的次线性速率收敛到目标的最小化点。注意，SAGA 的收敛是以平均迭代给出的，而不是最终迭代，这与 第 2 节 中的确定性算法形成对比。这是随机优化算法的一个相当普遍的问题，由证明技术引起。然而，在实践中通常使用最终迭代而不是平均迭代，无论缺失的理论保证如何。

SAG 遵循几乎相同的迭代，除了 \(\nabla h_{i_k}(x^{(k)}) - \gamma_{i_k}^{(k)}\) 前面没有 \(n\) 因子，这降低了方差并使 SAG 有偏。SAG 的原始收敛分析需要计算机辅助方法来完整证明，尽管收敛最近已无需计算机辅助而得到证明 [DLS22]。

当子集数量很大或数据维度很高时，存储历史梯度表可能变得不可行。然而，当应用于线性模型（在 [GSBR20] 的意义上）时，在某些应用中存储需求可以显著改善，见 第 4 节 以获取细节。

SVRG。SVRG 不存储梯度表，而是存储一个参考迭代点（通常称为锚点或快照），以及在参考迭代点处求值的全梯度。参考迭代点和全梯度都在迭代期间的预定义间隔处更新。为此，SVRG 需要两个循环：在外循环中，我们更新参考迭代点并计算全梯度；在内循环中，我们使用当前梯度和外循环中全梯度的组合来执行随机更新。内循环的长度按惯例设置为凸问题的 \(2n\) 和非凸问题的 \(5n\)。

SVRG 最初仅为光滑目标设计 [JZ13]。复合目标的近端版本在 [XZ14] 中开发和分析。他们表明，对于强凸问题，在 \(\tau < 1/(4nL_{\max})\) 和足够大的 \(t\) 的条件下，SVRG 的平均迭代在目标函数值中以线性速率收敛 [XZ14]。

SVRG 相对于 SAG(A) 的优势在于它只需要存储两个数据块，\(\tilde{x}^{(j)}\) 和 \(\gamma^{(j)}\)。因此，其存储需求与 \(n\) 无关。然而，其缺点是需要两个循环，这引入了额外的调参参数，以及每次更新参考迭代点时需要求值全梯度（SAGA 不需要）。

自适应子集大小。光滑项的数量 \(n\)，或等价地批量大小，对 SGD 的方差有直接影响，且可以为具有有界方差的无偏随机估计量的一般情况推导 [BCN18]。具体而言，随机估计量的方差随子集数量的增加而以大致线性的速率增加。因此，方差可以通过自适应或动态选择批量或子集大小来减少，而不是通过减小步长。这一主题主要从深度学习的角度进行了研究。然而，需要注意的是，增加 \(n\) 也会增加每次迭代的计算成本，因此可能遇到与 PGD 相同的可扩展性问题。因此，及时的批量增加在实践中至关重要，通常设计特定标准来解决这个问题，例如 [BCNW12]。虽然这一方法论的不同变体已在某些反问题中使用，如扩散光学层析成像 [MAP20] 和 PET [TAH+19]，但该主题尚未得到广泛采用。

3.2.3收敛模式与收敛速率

收敛模式和收敛速率的研究是随机迭代方法的关键重要主题，已在广泛假设下进行了研究，尽管这主要由机器学习社区针对相关优化任务完成。由于迭代的随机性质，随机迭代方法的收敛从概率角度进行研究。即，迭代（或其函数）被视为随机变量，可以考虑几种收敛模式之一：均值收敛（即在期望中）、几乎必然收敛（概率 1）或概率收敛（作为高概率保证）。对于随机方法，这些研究针对最后迭代或平均迭代进行。

收敛速率（如果存在）通常表示为迭代次数的函数（通常是幂），并严重依赖于目标的性质。历史上，随机迭代方法的收敛速率主要针对强凸目标的最后迭代发展，即当 \(h_i\) 和 \(g\) 都是强凸且 Lipschitz 光滑时。对于非强凸目标，收敛分析通常针对平均迭代进行。然而，这设置了一个不幸的困境，因为对于反问题，迭代方法最信息丰富的收敛模式是非强凸函数 \(h_i\) 和凸的、适当的、下半连续 \(g\) 的最后（或运行中）迭代的收敛。不幸的是，并非所有方法都建立了这种收敛速率。

随机梯度方法收敛的一个重要考虑是梯度方差对步长制度的影响。具体而言，为了控制梯度噪声，SGD 需要递减的步长，这阻碍了其收敛速度。方差缩减方法解决了这个问题，它们能够足够快地减少梯度噪声，同时保持常数步长以提供更快的收敛速率。在 表 3 中，我们总结了关于 Lipschitz 光滑 \(h_i\) 和适当的、凸的、下半连续 \(g\) 的收敛速率和收敛模式的最新技术，关于平均迭代和最后迭代以及凸性性质。

3.2.4其他基于随机梯度的方法

SGD 和方差缩减算法近年来针对感兴趣的特定场景进行了许多修改和适应。例如，Loopless-SVRG 算法 [KHR20] 是 SVRG 的修改，只使用一个循环，而不是在预定义间隔处更新参考迭代点和全梯度，在每次迭代中它从 \([0,1]\) 取随机样本。如果抽样值低于预定义的概率，则更新参考点和全梯度，否则执行标准更新。SARAH [NLST17] 是 SVRG 的有偏变体，它不仅在外循环中更新全梯度，还在内循环中更新。类似的有偏方差缩减梯度方案包括 SPIDER [FLLZ18, WJZ+19]。另见 [DLS22] 中名为 SARGE 的梯度估计量，它结合了 SAGA 和 SVRG 的思想。

第 2.2.2 节 中的三项和算法也承认随机变体。例如，PD3O 或 Condat-Vũ 算法的随机变体可以与 SGD 或 SVRG 以混合搭配的方式结合，通过用相应的随机估计量替换梯度 \(\nabla h(x)\)，见 [SCMR22]。类似地，对于其他包含梯度分量的第一阶方法，可以用随机估计量替换梯度，如原始-对偶不动点算法 [ZZ21]、随机 Davis-Yin 三算子分裂方法 [YVC16]。注意，这些随机方法在原始变量上使用随机性，例如当应用于 \(\nabla h(x)\) 时。对于涉及对偶变量的方法，我们将讨论推迟到 第 3.3 节，因为它们导致不同类型的随机优化方法。

用于 SGD、SAGA 和 SVRG 的原理也已应用于最大似然期望最大化（MLEM）和 MAP，见 [CM09, CZTZ18, KWML19]，并应用于层析图像重建 [KTA+21]。

3.2.5加速 SVRG 方法

由于方差缩减梯度方法可以使用常数步长并具有与其确定性对应物相同的收敛速率，自然会问 NAG 和 FISTA 等加速方法是否可以在此设定中应用。这是一个相当具有挑战性的任务，因为 NAG 和 FISTA 对误差非常敏感 [AD15]。沿此方向的首次尝试将 FISTA 与 SVRG 梯度估计量结合；然而，当问题不是强凸时，不提供理论保证，见 [Nit14]。后来在 [LMH15] 中，基于加速 PPA 形式，提出了一种称为 Catalyst 的方法来加速求解 (19)。Catalyst 包含两个循环：外循环是加速 PPA，而内循环用随机梯度方案求解到一定精度。在 [AZ17] 中，提出了一种名为 Katyusha 的加速方案来加速 SVRG。通过结合"负动量"，展示了非强凸问题的 \(O(1/k^2)\) 收敛速率。在 [DES22] 中，提出了一种更优雅的加速方案，见算法 10，作者表明不需要负动量就能达到 \(O(1/k^2)\) 收敛速率。

近年来，向量外推方案也被采用到随机优化中，例如非线性加速 [SBd17] 方案（如 Anderson 加速及其变体）被应用于加速随机梯度方法。这类方案与拟牛顿方案（如 L-BFGS）密切相关，如 [Sci19] 所示。

3.2.6增量梯度方法

在 2010 年代机器学习随机优化算法浪潮之前，一个相关的优化算法家族已经受到越来越多的关注：增量梯度方法。与 SGD 一样，增量梯度方法也设计用于求解具有有限和结构 (19) 的目标函数。核心思想是在迭代的每一步中，算法通过预定规则（循环或固定排列）从有限和中选择一个函数，并对所选函数应用 GD [Ber11]。类似于次梯度下降和 SGD，需要仔细选择的步长来确保收敛：\(\sum_k \tau_k = \infty\)，\(\sum_k \tau_k^2 < \infty\)。这一想法用于临床 PET 成像中使用的改进块序列正则化期望最大化（BSREM）和松弛 OS 可分离抛物面代理 [PY01, AF03]。

为了规避递减步长的需要，开发了 [BHG07] 中的增量平均梯度（IAG）方案，它可以承受常数步长。该方法的基本思想是平均历史梯度，这也是随机平均梯度方法 SAG 和 SAGA（见算法 8）背后的关键原理 [SLRB17, DBLJ14]。

这种方法可以扩展到非光滑优化（例如用于 PGD）。例如，当非光滑函数可以分解为有限和时，可以选择和中的一个项并计算其近端算子。我们参考专门的综述 [Ber11] 以获取更多细节。

最后，与它们的随机对应物相比，增量梯度方法具有较弱的理论保证，如更具限制性的步长标准 [Ber11, BHG07, SLRB17]。在 PET 成像背景下的 [TAK+23] 数值研究表明，BSREM 在每个考虑的指标上都劣于随机方差缩减算法。见 第 5 节 以获取增量梯度方法与随机方法的比较。我们补充说，以循环方式访问梯度索引的增量梯度算法对于某些病态用例可能表现出指数级的坏行为 [HM93, RR12, Def15]。

3.3随机（对偶）坐标下降

在反问题的背景下，随机梯度方法通常降低测量空间中的复杂度，特别是当光滑项 \(h\) 编码数据保真度时，如 \(\frac{1}{2}\|Kx - v\|_2^2\)。另一种选择是降低优化变量空间的复杂度。

回忆 第 2.1.1 节 中的 CD，它执行迭代 \[ x^{(k+1)} = x^{(k)} - \tau_k \nabla_i h\bigl(x^{(k)}\bigr), \] 其中 \(\nabla_i h(x) = \frac{\partial h(x)}{\partial x_i} e_i\) 且 \(e_i\) 是第 \(i\) 个标准基向量。与确定性情况一样，我们可以使用（随机选择的）坐标子集，而不是使用单个坐标。

注意，方向 \(\nabla_i h(x)\) 是 \(\nabla h(x)\) 的无偏估计。此外，它也是方差缩减的。如果 \(x^{(k)} \to x^*\)，则由于 \(x^*\) 的最优性条件，我们得到 \(\nabla h(x^{(k)}) \to 0\) 和 \(\nabla_i h(x^{(k)}) \to 0\)，因此 \(\|\nabla_i h(x^{(k)}) - \nabla_i h(x^*)\|^2 \to 0\)。

为了提供计算成本的降低，CD 通常需要坐标友好的结构（如可分函数）。即，坐标更新应该显著比梯度计算便宜。因为与反问题相关的优化问题通常不具有这种结构，所以 CD 方法在此背景下的应用 largely 受限。一个例外是 [GB18]，它在 GPUs 上的大规模 CT 图像重建中使用 CD 来克服内存瓶颈。事实证明，当应用于对偶变量时，CD 对于反问题更有用，如下所述。

虽然随机性已被观察到可以改善梯度估计量的性能，与使用确定性（例如循环）梯度抽样相比，但 CD 的情况更为微妙。一方面，存在非凸问题的类别，对于这些问题循环 CD fails 收敛到平稳点 [Pow71]，但随机 CD 期望在仅几步内收敛，在其他情况下随机 CD 可以显著快于循环 CD [Wri15]。另一方面，可以设计前向算子的类别，使得循环 CD 对于相应的二次问题比随机 CD 快几倍 [GOPV17]。对于块 CD 的循环、随机和自适应选择的进一步比较，我们参考 [NLS22]。

3.3.1随机原始-对偶分裂

为简单起见，我们将限制在 \(m = 1\) 和 \(n = 0\) 的情况，即我们的模板 (17) 取形式

具有关联的鞍点问题

当将 PDHG（算法 3）应用于这个问题时，每次迭代需要为所有 \(A_i\) 求值矩阵-向量乘积以更新原始变量，为 \(A_i^*\) 求值以更新对偶变量，这在计算上可能是昂贵的。因此，在 [CERS18] 中提出改为执行对偶坐标上升。所得算法称为 SPDHG，在算法 11 中提供。

SPDHG 在 Bregman 距离意义下收敛，条件为 \(\tau\sigma < 1 / (\ell \max_i \|A_i\|^2)\)。此收敛结果可以扩展，算法可以推广到其他抽样策略和步长制度，见 [CERS18, AFC22, GDE21] 以获取更多细节。

注意，由于 \(A\) 的线性性，我们可以将原始更新重写为

注意，这与 SAGA 估计器关于 (23) 中双线性项 \(x\) 的梯度非常相似。然而，由于 \(1+\ell\) 而不是 \(\ell\) 的因子，这是一个有偏梯度估计量。其在各种应用中的大规模问题的有效性已得到证明，如具有正弦图数据的 PET [EMS19]、具有列表模式数据的 PET [SH22]、光声断层成像 [RPDC20]、CT [CDE+24, RDH21, PAD+21]、多光谱 CT [PAD+21]、并行 MRI [GDE21, GDE22] 和磁粒子成像 [ZB21]。

与任何原始-对偶算法一样，调整原始和对偶步长的比例是困难的。[ZB21] 中对此进行了实证处理，[CDE+24] 提供了理论基础。注意，问题也可以通过其他原始-对偶算法类似地处理，例如 [SSZ13, FR15, FB19]。

3.4延伸阅读

随机二阶方法。标准二阶方法在大规模问题中应用有限，由于计算目标 Hessian 的代价和不可行性。如 第 2.3 节 所讨论的，确定性拟牛顿方法使用关于曲率的信息而不计算 Hessian。相反，使用近似技术来近似 Hessian 或相应的矩阵-向量乘积。将拟牛顿方法直接扩展到随机设定可能具有挑战性，因为曲率的噪声估计（与梯度的随机估计量结合）可能 destabilize 算法。已经提出了许多方法，通过利用拟牛顿和随机梯度方法的优势来解决这个问题。

该领域的开创性工作是随机拟牛顿方法 [SYG07]，它引入了 BFGS 和 L-BFGS 的随机修改，展示了次线性收敛速率。这一限制通过 [MNJ16] 得到解决，通过随机 L-BFGS 变体 [BHNS16] 与方差缩减的结合，在温和假设下达到了线性收敛速率。

子采样拟牛顿方法，见例如 [RKDA14, EM15, XYR+16, RKM19]，随机子采样 Hessian 以进一步降低成本，并使用矩阵集中不等式来确保曲率信息得到良好保持。尽管取得了这些进展，随机拟牛顿方法仍未在反问题成像社区中广泛采用，有显著的例外 [MXC+19, HPAT22, PA21]，它们继续是一个活跃的研究领域。

前向算子 sketching。sketching 方法的基本思想是用一个更小、因此更易于管理的问题来替代原反问题。这通常针对最小二乘问题 \(\min_{x \in \mathbb{R}^d} \|Kx - v\|_2^2\) 进行研究。思想是找到一个 sketching 矩阵 \(S \in \mathbb{R}^{m \times s}\)，其中 \(m \ll s\)，使得 \[ \min_{x \in \mathbb{R}^d} \|S(Kx - v)\|_2^2 \] 的解在某种意义上与原问题的解相似，但计算成本大幅降低。常见的 sketching 技术包括随机投影、子采样和杠杆得分抽样 [Woo14, DMM06, Mah11]。然而，这种 sketching 是有代价的，因为有不可逆的信息损失 [PW16]。

迭代 Hessian sketching（IHS）。通过 sketching 求解反问题的另一种方法是 sketching Hessian。类似于牛顿法，IHS 使用 Hessian 构建原问题的二次近似，然后对其进行 sketch。更详细地，考虑最小化光滑凸目标 \(h(x)\)。在每次迭代中，sketching 矩阵 \(S \in \mathbb{R}^{s \times d}\) 用于构建 sketched Hessian 矩阵 \(H_S = (H^{1/2})^* S^* S H^{1/2}\)，其中 \(H\) 是当前迭代 \(x^{(k)}\) 处的 Hessian 矩阵。IHS 然后遵循迭代 \[ x^{(k+1)} = x^{(k)} - \tau_k H_S^{-1} \nabla h\bigl(x^{(k)}\bigr), \] 其中 \(\tau_k > 0\) 是步长。对于某些凸优化问题，IHS 比基线方法更快 [PW16, PW17]。它也可以与（加速）PGD 结合 [TGD17]，并已用于加速多线圈 MRI 重建 [OOI+24]。

虽然它们有一些明显的优势，但基于 sketching 的方法很少应用于反问题成像。这主要归因于 sketching 可能破坏前向算子中存在的任何有益结构的事实。例如，在 CT 中，Radon 变换是稀疏矩阵，我们只能使用子采样作为 sketching，导致次优性能。MRI 算子包含 FFT，这再次阻止了 sketching 的朴素应用。最近在 [OOI+24] 中，开发了一种特殊的 sketching 算子，可以保留 FFT 的快速计算结构，并应用于多线圈 MRI 重建任务。

随机镜像/Bregman 下降。近年来，SGD 已扩展到不依赖欧几里得结构的问题。随机镜像下降（SMD）是一种强大的这种推广，在每一步中可以以近端点形式写成

其中 \(\varphi\) 是镜像映射，\(\tilde{\gamma}_k\) 是 \(x^{(k)}\) 处目标梯度的估计量。SMD 通常通过非光滑 Lipschitz 连续凸函数的随机方案的视角进行研究 [CBCG04]，尽管其收敛性可以对非光滑凸目标 [GL12]、光滑目标 [Bub15] 和其他目标条件进行证明，并已应用于不适定反问题 [JLZ23]。近年来也研究了自适应步长方案 [DLLM23]。

SGD 也在 \(\mathcal{V}\) 和 \(\mathcal{U}\) 都是 Banach 空间的问题上进行了研究，通过应用类似近端点的算法 [JK23]，在 CT 中有应用。在这种情况下，选择 \(\mathcal{V}\) 和 \(\mathcal{U}\) 为某些 \(L^p\) 或 Sobolev 空间可以作为解的隐式正则化器，并允许适应不同的噪声假设。然而，选择 \(\mathcal{V}\) 和 \(\mathcal{U}\) 以优化重建质量是一个具有挑战性的问题。这可以在某种程度上通过使用模（或变指数）Banach 空间 [LKEC23] 来缓解，它可以适应空间可变噪声类型的情况。

训练神经网络的自适应步长方法。SGD 在（深度）神经网络训练背景下的一个重要进展是自适应步长随机梯度方案，如 AdaGrad [DHS11]、AdaDelta [Zei12]、RMSProp [TH12]、ADAM [KB15]，仅举数例。我们参考 [Rud16] 以获取这些方法的简要概述。

这类方法的主要动机是观察到某些变量或解的组件在训练期间的更新明显少于其他变量。更精确地，令 \(d^{(k)} = \widetilde{\nabla} h(x^{(k)})\) 表示第 \(k\) 次迭代中计算的随机梯度。然后假设对于迭代 \(k\)，存在 \(1 \leq i \leq d\)，使得 \(|d_i^{(k)}|\) 与 \(j \neq i\) 的 \(|d_j^{(k)}|\) 相比较小。在这种情况下，解的相应 entry \(x_i\) 将需要显著更长的训练时间。自适应步长方法通过使用为 \(x\) 的每个 entry 定制的步长来克服这个问题。对于对角矩阵 \(H_k = \diag(1/(|d_i^{(k)}| + \varepsilon))\)，自适应步长 SGD 可以遵循迭代

因此，如果对于某个 \(i\)，\(|d_i^{(k)}|\) 非常小，则 \(H_k\) 的相应 entry 将很大，因此 entry \(x_i^{(k+1)}\) 将使用更大的有效步长并接受更大的更新。

与标准 SGD 方案相比，一般而言，自适应步长方案有两个优势。首先，步长调优变为自动的，使得每个 entries 以相对相同的速度收敛。其次，步长选择是免调参的，例如在 (24) 中我们只需要选择 \(\tau > 0\)。在过去十年中，自适应步长 SGD 方案在训练神经网络中变得更加流行，而标准 SGD 获得了更深入的理论理解。

随机化近端算子。第 3 节中的算法通过使用低成本的梯度估计量来降低迭代复杂度。类似地，第 3.3 节中的算法在对偶空间中使用低成本的更新。这加速了反问题中的计算，因为它降低了应用算子 \(K\) 或 \(A\) 固有的计算成本。然而，计算成本的第二个来源在于近端算子的求值，特别是对于在每 次迭代中计算 \(\prox_{\tau g}\) 或 \(\prox_{\tau f}\) 的算法。由于随机算法需要更多迭代来执行相同数量的数据遍历（并利用所有测量数据），如果 \(\prox_{\tau g}\) 的求值代价高昂，这可能产生进一步的计算成本。

为了解决这个问题，近年来开发了新的方法，将近端域的随机方法。在 [CR23] 中，RandProx 算法被引入，它在原始-对偶 Davis-Yin 算法中使用随机估计量来进行近端步骤。随机估计量在其最简单的形式中由一个概率参数 \(p\) 定义，即取决于抛硬币的结果，近端算子要么应用要么不应用。ProxSkip [MMSR22] 将此方法论直接应用于近端梯度风格算法。然而，这些方法主要应用于联邦或去中心化学习，它们在反问题成像中的应用尚未得到充分研究。

随机迭代正则化。随机优化器的正则化性质是机器学习社区中已得到充分观察和研究的一个性质。对于线性反问题和最小二乘型目标，已证明 SGD [JL18, JZZ21] 和 SVRG [JZZ22] 对于先验选择的适当停止准则具有正则化性质。类似的结果已针对 SGD 的镜像和 Banach 空间变体 [JLZ23, JK23]、随机化 Kaczmarz 方法 [NT14, TL23] 和其他随机优化方案得到证明。

非线性反问题。随机梯度方法已成功应用于一些结构适当的非线性反问题。例如，随机化 Kaczmarz 方法及其受 SVRG 启发的方差缩减变体已成功应用于相位恢复问题 [XLTW22]、EIT [vdDA12, ZQTL24]、光学层析成像 [MAP20] 和其他问题。然而，一般而言，将随机梯度方法应用于非线性问题提出了额外的挑战，由于前向模型及其性质引入的复杂性。此外，非线性反问题通常需要迭代依赖的步长，通常通过回溯找到。随机算法的回溯是一个开放的研究问题，没有通用的解决方案；另见 第 4 节。这些问题通常需要针对问题的解决方案，该领域需要更多研究。

4反问题中随机优化的挑战

在前两节中，我们提供了确定性和随机优化算法的概述。在本节中，我们对随机优化方法的各种理论和实践方面进行深入讨论。

4.1随机梯度方法何时及为何快速？

如 第 3 节 开头所述，随机优化算法的目标是处理大规模优化问题。人们自然会问随机优化方法是否保证比其确定性对应物更快。以下是一些要点。

随机加速因子。反问题是否适合应用随机梯度方法，很大程度上取决于优化问题的结构。

在 [TED19, TEGD20] 中，提出了一个度量指标来估计反问题是否适合应用随机梯度方法。这个度量指标，称为随机加速因子（stochastic acceleration factor），定义为全光滑项梯度的 Lipschitz 常数与子集光滑项的 Lipschitz 常数之比，\(\Upsilon := L / L_{\max}\)。对于随机方法的理想情况，这个比值接近 \(n\)。许多具有常见数据集的机器学习问题落入这一类别，以及一些重要的反问题成像问题，如 CT 重建。然而，也有许多负面例子，其中 \(L \approx L_{\max}\)，使得 \(\Upsilon \approx 1\)，在这种情况下不推荐随机梯度方法。这些包括压缩感知反问题、图像去噪或图像修复。对于图像去模糊，\(\Upsilon\) 强烈依赖于核的大小，另见 第 5 节 中的数值示例。

在线性反问题的背景下，数据拟合通常由最小二乘给出，专门的分析揭示了随机加速的成功依赖于前向算子的结构。[TEGD20] 提供的理论分析表明，当且仅当反问题的 Hessian 具有快速衰减的谱时，反问题适合随机梯度方法。在优化社区中，通常被忽视的是随机梯度方法利用数据矩阵的固有低维结构。如果在某些问题中这种结构较弱，使用随机算法可能不会带来任何好处。

方差缩减方法何时有益？方差缩减方法的重要性高度依赖于测量中的噪声量。已证明在 [VML+19] 中，SGD 在所谓的插值制度中以线性收敛，在反问题的背景下，这本质上等同于无噪声测量制度。在这种情况下，SGD 甚至可以比方差缩减方法更快，因为它可以使用 \(1/(nL_{\max})\) [GLQ+19] 的常数步长。相比之下，方差缩减方法通常需要更小的常数步长。例如，已知最大的可证明收敛步长是 SAGA 的 \(1/(3nL_{\max})\) [DBLJ14] 和 SVRG 的 \(1/(4nL_{\max})\) [XZ14]。当测量噪声不可忽略时，朴素 SGD 需要使用递减步长以确保收敛，对于高精度解变得缓慢。在这种情况下，推荐方差缩减随机梯度方法，因为它们在这种对反问题普遍的制度中具有快速收敛。另见 第 5 节 中的数值实验。

高效部分计算。如上所述，测试随机加速因子是 SGD 对特定反问题的理论效用的必要条件。在这里，我们描述另一个与可并行性相关的测试。实际上，反问题的所有随机方法都将前向算子 \(K\) 划分为更小的块 \(K_i\)，使得 \((Kx)_i = K_i x\)。一个根本重要的问题是 \(K_i x\)（\(i = 1, \ldots, s\)）是否能与 \(Kx\) 一样快或更快地计算。随机方法本质上将问题序列化，这增加了计算开销。随机方法的实际可行性然后取决于随机加速因子和前向算子的可并行性。例如，大多数 X 射线变换的实现可以在其角度上并行化，并行 MRI 中的前向算子可以在线圈上并行化。

子集大小与计算开销。子集数量允许在确定性和高方差随机算法之间进行插值，两者都提供不同的好处和缺点。为了说明这一点，我们比较 PGD 与 SGD。在 PGD 的每次迭代中，最昂贵的操作是梯度和近端算子的求值。在成像反问题中，通常使用计算成本可忽略不计的近端算子，例如 TV 的近端算子，可以通过迭代 [BT09, RC20] 或应用深度神经网络 [HLP21, TMTS24b] 来近似。因此，SGD 的单个数据遍历具有相同的梯度计算成本（因为它仅在子集上求值），但它需要 \(n\) 次求值近端算子，这通常导致相当大的开销。因此，对于大多数应用，子集数量 \(n\) 需要保守选择以避免这种计算开销。

4.2超参数调优

参数调优一直是随机优化中一个重要且长期的挑战。由于它们的随机性质，大多数算法参数通常比其确定性对应物更难先验选择或即时调整。通常需要调优的典型算法参数包括以下内容。

步长选择。设置适当的步长是随机优化中的一个关键挑战，并一直是开发高级随机算法（如方差缩减随机算法）的驱动力。一方面，随机算法继承了确定性算法的步长问题，例如需要知道 Lipschitz 常数，这可能不可用、计算精确代价高昂或高度空间变化。另一方面，有理论保证的步长通常过于保守，在实践中，大于理论极限的步长通常工作良好并导致更快的性能。例如，在 第 5 节 的第一个实验中，SAGA 使用步长 \(1/(1.75nL_{\max})\)，在该步长下算法收敛。类似地，原始-对偶算法需要适当平衡原始和对偶步长以获得良好的性能。总体而言，随机优化算法的步长选择需要更多专门研究。

回溯。虽然线搜索方案通常可以为确定性梯度方法（如 PGD 或 FISTA [BT09]）实现，但将这些方案应用于随机梯度方法通常很困难，无论是理论上还是实践中。传统线搜索方法需要求值完整目标，这是不可取的，因为这需要使用所有数据。使用相应的随机变体在实践中可以表现良好，但缺乏理论保证 [GSBR20]。最近，在 [DM23] 中为机器学习中的随机梯度方法的步长即时调优取得了许多突破。然而，这尚未完全适应近端设定，这是反问题中普遍存在的。

子集大小。子集大小的选择通常是实践中的一项棘手任务。虽然在理论上，当选择微小的子集大小时，随机方法的计算复杂度是最优的，但在实践中还必须考虑大多数现代计算设备高度并行的事实，这意味着硬件会强烈偏向较大的子集，导致一个权衡。目前，关于提供子集大小选择实用指南的研究有限。

子集模式。选择子集模式，即前向算子的划分，对于最优应用随机梯度方法至关重要。如 [TEGD20] 所讨论的，子集模式决定了子集上的梯度 Lipschitz 常数，因此影响可容许的步长和收敛速度。子集选择的影响也在 [GDE22] 中针对并行 MRI 进行了讨论。合适子集模式的选择对于某些反问题（如全角度 2D CT）是 reasonably 良好理解的。然而，对于许多其他问题，如有限角度 CT、3D PET 等，这个问题尚未解决。

重启。加速随机梯度方法不会自动适应（受限的）强凸性，这是使用稀疏诱导正则化项时普遍存在的。没有重启，这些方法仍然显示快速收敛，但只是初始的。当接近高精度解时，算法开始在解附近振荡，减慢收敛。为了解决这个问题，已开发了自适应重启方案来加速基于动量的方法 [TGBD18, FQ20]。

热启动。在实践中，方差缩减方法如 SVRG 和 SAGA 在开始时比朴素 SGD 慢。通过用 SGD 的初始迭代替换方差缩减方法的迭代，在迭代开始时观察到更快的收敛 [TAK+23, KR17]。然而，先验地不清楚这个阶段应该持续多长时间。

4.3内存考量

这里介绍的随机算法具有非常不同的内存占用，这可能在它们在应用中的使用方式中发挥重要作用。在下文中，我们将讨论内存需求的影响，以及 SGD、SAG(A)、SVRG 和 SPDHG 的一些缓解策略，尽管类似的考虑可以应用于其他算法。

朴素 SGD不需要存储任何额外的信息，因此具有与 GD 相同的 \(O(\dim \mathcal{X})\) 内存占用，仅用于存储当前迭代。

SAG 和 SAGA在其标准形式中，需要存储计算梯度的表。表中元素的和然后需要在每次迭代中存储或重新计算。因此，标准实现的 SAG 和 SAGA 的内存开销是 \(O(n \cdot \dim \mathcal{X})\)，当 \(n\) 和 \(\dim \mathcal{X}\) 很大时，这可能大得令人望而却步。特别是在大规模 CT 中，维度 \(\dim \mathcal{X}\) 可能非常大。然而，在某些缓解场景中，内存需求可以显著降低。考虑线性模型 [DBLJ14, GSBR20]，即当 \(h_i = h \circ A_i\) 对于 \(i = 1, \ldots, n\) 和线性算子 \(A_i\) 时。这例如用于最小二乘数据保真度或 Kullback-Leibler 散度（见 第 1 节）。\(h_i\) 的梯度满足 \(\nabla h_i(x) = A_i^* \nabla h(A_i x)\)。换句话说，我们可以存储 \(\nabla h_i(A_i x)\) 而不是 \(\nabla h_i(x)\)。如原始 SAGA 中一样，每次迭代需要对 \(A_i\) 和 \(A_i^*\) 各应用一次，因此不引入额外的计算成本。修改后的 SAGA 算法在算法 12 中写出。

修改后的 SAGA 的存储需求是 \(O(\dim \mathcal{X} + \dim \mathcal{Y})\)。注意，存储的梯度和在 \(\mathcal{X}\) 中，而对偶变量 \(y_i^{(k)}\) 在 \(\mathcal{Y}_i\) 中。因此，存储所有 \(y_i^{(k)}\)（\(i = 1, \ldots, n\)）可以用 \(\dim \mathcal{Y}\) 的内存完成。此外，我们可以利用模型的线性，使得如原始 SAGA 中一样，每次迭代只需要一个算子 \(A_i\) 的伴随求值。

SAGA 的首选变体取决于成像模态和选择的子集数量。例如，在 CT 中，图像和数据大小通常相似，使得修改后的 SAGA 算法即使对于小的 \(n\) 也更可取。对于 PET，图像体积通常比数据小一个数量级，因此对于 第 1 节 中的设定，当子集数量超过 15 时，修改后的算法开始变得可取。

SVRG的标准实现（见算法 9）需要存储参考迭代点 \(\tilde{x}^{(j)}\) 和全梯度 \(\gamma^{(j)}\)，产生 \(O(\dim \mathcal{X})\) 的内存，与 \(n\) 无关，与 GD 相当。然而，这以每次迭代计算两个随机梯度为代价，一个在先前迭代点，另一个在参考迭代点。对于 SAG(A)，我们可以推导出具有相同 \(O(\dim \mathcal{X} + \dim \mathcal{Y})\) 内存占用和与修改后 SAGA 相同计算成本的 SVRG 修改版本。

SPDHG具有本质上与 PDHG 相同的内存占用：\(O(\dim \mathcal{X} + \dim \mathcal{Y})\)。取决于对偶空间，原始-对偶方法的内存占用与 GD 相比可能显著更大。

5数值实验

在本节中，我们展示一个简短的数值研究，展示本综述中讨论的代表性随机优化算法，并将其与其确定性对应物进行比较。注意，此数值研究并不旨在提供算法性能的全面和彻底比较。相反，我们的目标是展示随机算法确实对某些问题有益，并说明它们的行为。本文中所有结果的复制代码已公开可用，见文末的数据可用性声明。

5.1稀疏脉冲去模糊

我们考虑一个玩具去模糊问题，其中真实值是稀疏脉冲，前向算子使用均匀模糊核。此示例的目标是比较 PGD 和 SAGA，并研究这些算法的性能如何受模糊核大小和步长制度的影响。数值结果用 MATLAB 复制。

令 \(u^{\dagger} \in \mathbb{R}^d\)，\(d = 1000\) 为具有 20 个均匀分布非零 entries 的稀疏向量，代表真实值，\(K \in \mathbb{R}^{d \times d}\) 为大小 \(\kappa\) 的均匀核模糊的矩阵表示，具有周期性边界条件，归一化为 \(\|K\| = 1\)。我们考虑从含噪数据 \(v = Ku^{\dagger} + \eta\) 中恢复 \(u^{\dagger}\)，其中 \(\eta \in \mathbb{R}^d\) 是随机高斯噪声。我们考虑大小为 \(\kappa_j = 2j - 1\)（\(j = 1, \ldots, 100\)）的核，和 10 个不同噪声水平的随机高斯噪声 \(\eta \sim \mathcal{N}(0, \sigma^2 I)\)，在 \([10^{-4}, 10^{-1}]\) 中等距取点。

为了从观测 \(v\) 中恢复 \(u^{\dagger}\)，我们最小化 \(\Phi(x) = g(x) + h(x)\)，其中 \(x = u\)，\(g(x) = \mu\|x\|_1\) 且 \(h(x) = \frac{1}{2}\|Kx - v\|_2^2\)。对于 SAGA，目标 \(\Phi\) 被重写为有限和形式 \(\Phi(x) = g(x) + \sum_{i=1}^{n} h_i(x)\)，其中 \(n = d\)，\(h_i(x) := \frac{1}{2}(K_i x - v_i)^2\) 且 \(K_i \in \mathbb{R}^{1 \times n}\)（\(i = 1, \ldots, n\)）是 \(K\) 的第 \(i\) 行。注意 \(L = 1\)，\(\|K_i\| = 1/\sqrt{\kappa}\)（\(i = 1, \ldots, n\)），因此对于给定核大小 \(\kappa\)，\(L_{\max} = 1/\kappa\)。正则化参数选择为 \(\mu = \frac{1}{2}\|K^* v\|_{\infty}\)。

对于两种方法，我们使用两种步长选择。对于 PGD，我们使用 \(\tau = 1/L\)（记为 PGD1）和 \(\tau = 1.9/L\)（记为 PGD2）。对于 SAGA，我们使用 \(\tau = 1/(3nL_{\max})\)（记为 SAGA1）和 \(\tau = 1/(1.75nL_{\max})\)（记为 SAGA2）。

对于每个核大小和噪声水平，我们计算一个高精度解 \(x^*\)，然后运行 PGD 和 SAGA，使用停止准则 \(\|x^{(k)} - x^*\| / \|x^*\| \leq 10^{-4}\)。我们记录达到该精度所需的数据遍历次数，记为 \(N_{\text{PGD}}\) 和 \(N_{\text{SAGA}}\)，并计算比值 \(N_{\text{SAGA}} / N_{\text{PGD}}\)，对 10 个噪声水平取平均。图 3 展示了结果。

5.2模拟平行束 CT

我们在此复现了即将发表的关于 Core Imaging Library（CIL）中大型成像问题随机优化方法软件框架的数值结果。这正在 CCPi 和 SyneRBI 项目周围作为社区努力开发。复现的结果涉及大小为 \(128 \times 128\) 像素的 Shepp-Logan 体模的 CT 重建。测量遵循前向问题 \(v = Ku^{\dagger} + \eta\)，其中 \(\eta \sim \mathcal{N}(0, \sigma^2 I)\)，前向算子 \(K\) 由范围 \([0, \pi)\) 中 240 个角度的平行束 CT 定义。

我们比较确定性算法 PGD（算法 1）与 SGD（算法 7）和 SAGA（算法 8）的性能。如 第 3 节 所讨论的，求解反问题时，通常不使用完全随机的方法，在这个意义上，不是使用随机构建的批次，而是将测量数据预划分为批次。然后在每次迭代中随机选择一个划分元素。

为了恢复体模 \(u^{\dagger}\)，我们最小化 \(\Phi(x) = g(x) + h(x)\)，其中 \(x = u\)，\(g(x) = \lambda \text{TV}(x)\)，且 \(h(x) = \frac{1}{2}\|Kx - v\|_2^2\) 用于 PGD。SGD 和 SAGA 将目标重写为 \(\Phi(x) = g(x) + \sum_{i=1}^{n} h_i(x)\)，其中 \(n \in \mathbb{N}\)，\(h_i(x) = \frac{1}{2}\|K_i x - v_i\|_2^2\) 且 \(\sum_{i=1}^{n} h_i(x) = h(x)\)。\(K_i\) 算子（\(1 \leq i \leq n\)）定义为使用范围 \([(i-1)/n, \pi)\) 中 \(240/n\) 个角度的 CT 算子，从 \((i-1)/n\) 开始，每第 \(n\) 个角度取增量 \(\pi n/240\)。这有时被称为等距或交错划分。

初始图像为 \(x^{(0)} = 0\)。PGD 的步长为 \(\tau = 1/\|K\|_2^2\)。对于 SGD，我们使用递减步长计划 \(\tau_k = 1/(2nL_{\max}(1 + 0.01k/n))\)，其中 \(L_{\max} = \max_i \|K_i\|_2^2\)，遵循 [TAK+23]，满足条件 (21)。对于 SAGA，我们使用步长 \(\tau = 1/(3nL_{\max})\)。

5.3核桃的扇形束 CT 数据

在这组实验中，我们取核桃锥形束 CT 扫描的中心切片。与 第 5.2 节 一样，我们复现了 CIL 软件库中随机框架即将发表的数值结果。完整正弦图数据在范围 \([0, 2\pi)\) 中有 721 个投影。为降低计算成本，我们只考虑范围 \([0, 2\pi)\) 中的 360 个投影，角度分离为 \(1^\circ\)。重建图像大小为 \(280 \times 280\) 像素。

含噪测量是通过遵循 Beer-Lambert 定律破坏正弦图（测量数据）而生成的。我们首先通过 \(N_I = \text{Pois}(I_0 \exp(-v))\) 计算期望计数，其中 \(v\) 是加载的（无噪声）正弦图，\(I_0\) 是初始光束的强度。然后，通过首先将 \(N_I\) 中的所有 0 entries 校正为 1，然后应用线性化（对数后测量）作为 \(v_\delta = -\ln(N_I/I_0)\) 来计算含噪测量。我们考虑三个强度水平：低强度（\(I_0 = 50\)）、中等强度（\(I_0 = 250\)）和高强度（\(I_0 = 5000\)），模拟从高到低噪声。

目标函数定义为

其中 \(\alpha_I \in \{0.124, 0.045, 0.01\}\) 分别是高、中、低噪声水平的正则化参数。参考解示于图 7。我们研究基于梯度的方法（PGD、SGD、SAGA、SVRG）和原始-对偶方法（PDHG 和 SPDHG）。

5.4模拟 PET 数据的图像重建

迭代方法的收敛行为对 Lipschitz 常数的敏感性可能在 Lipschitz 常数不可用、难以估计或高度变化的重建任务中产生严重后果。一个典型情况是 PET，其中数据保真项是泊松对数似然，或等价的 Kullback-Leibler（KL）散度。这个数据保真项通常不是全局 Lipschitz 光滑的。然而，在 PET 的情况下，对于任何具有非零散射和背景事件的真实 PET 扫描，它是 Lipschitz 光滑的。然而，相应的 Lipschitz 常数是悲观的，通常在实践中不使用。尽管如此，随机（和非随机）梯度方法仍用于 PET，尽管它们的使用在原则上更具启发性，但随机方法可以展示快速和可靠的性能 [TAK+23, KTA+21, CERS18]。

在优化模板 (7) 的术语中，我们将目标写成如下，其中 \(r\) 表示背景和散射事件。

基于梯度的方法：\(g(x) = 0.1\|\nabla x\|_1 + \iota_{[0, \infty)^d}(x)\)，\(f \equiv 0\)，\(h(x) = \KL(v \mid Kx + r)\) 用于 FISTA，\(h_i(x) = \KL(v_i \mid K_i x + r_i)\) 用于 SGD、SAGA 和 SVRG。

SPDHG：\(g(x) = 0.1\|\nabla x\|_1 + \iota_{[0, \infty)^d}(x)\)，\(A_i = K_i\) 且 \(f_i(y) = \KL(v_i \mid y + r_i)\)，\(h \equiv 0\)。

6展望

优化算法是求解反问题的支柱，而随机方法在大规模应用中正变得越来越重要。我们讨论了具有工具箱特性的前沿算法：几乎所有相关算法都是由一小部分来自确定性优化工具箱和随机优化工具箱的工具组装而成。特别地，随机工具箱是为机器学习而开发的，但事实证明它对反问题同样具有极好的实用性。

在未来十年，我们预计这一趋势将继续，越来越多应用将意识到随机优化提供的潜力。此外，我们预计专门为反问题开发的随机优化工具，这些工具可能与机器学习无关。潜在的例子包括图像域中的 sketching 和非线性前向算子带来的非凸性。我们还预计可以专门为反问题发展的理论。例如，已证明 SGD 在所谓的插值制度中可以快得多 [VML+19]。在反问题的背景下，这一制度意味着没有噪声且没有正则化，这当然是不现实的。相反，可能存在专门为反问题定制的理论，解释在特定噪声制度中哪些算法最快，从而可能解释我们在 第 5 节 中进行的观察。

最后，正如我们在手稿各处的进一步阅读部分中所指出的，一条令人振奋的前进路径是反问题、优化与机器学习之间的相互交叉与融合。这可以在建模的背景下进行，例如学习正则化模型或学习前向算子或其采样，也可以在学习的背景下进行，学习求解反问题的最优算法。我们相信，当前基于优化的范式与学习显式正则化将成为反问题中的主导计算范式，而为这些应用开发高效的优化算法将继续对从业者很重要。据我们所知，这尚未在随机算法的背景下进行尝试。在我们看来，这一发展远未结束，需要更多研究来充分理解和利用反问题背景下的联系。

附录

数学预备知识与符号

贯穿本文，\(\mathcal{X}\) 和 \(\mathcal{Y}\) 表示配备内积 \(\langle \cdot, \cdot \rangle\)和关联范数 \(\|x\| = \sqrt{\langle x, x \rangle}\) 的一般实 Hilbert 空间。从上下文可以清楚使用的是哪个内积。

线性有界算子 \(A: \mathcal{X} \to \mathcal{Y}\) 的（Hermitian）伴随记为 \(A^*\)，\(I\) 是恒等算子。我们使用 \(\|A\| := \sup_{\|x\|=1} \|Ax\|\) 表示诱导算子范数。

由 \((a, b]\) 我们表示区间 \(a < x \leq b\)。如果 \(\mathcal{X} = \mathbb{R}^d\)，我们写 \(x_i\) 表示向量的第 \(i\) 个分量。对于 \(1 \leq p < \infty\)，我们表示向量 \(x \in \mathbb{R}^d\) 的 \(p\)-范数为 \(\|x\|_p = (\sum_{i=1}^{d} |x_i|^p)^{1/p}\)。\(I\) 表示恒等算子。

凸性

函数 \(f: \mathcal{X} \to \mathbb{R}_{\infty}\) 被称为凸的，如果

且如果上述不等式是严格的，则它是严格凸的。如果

则它是 \(\mu\)-强凸的。如果对于某个 \(L > 0\) 成立

则函数 \(f\) 被称为 \(L\)-Lipschitz 连续的。如果 \(f\) 是可微的且具有 \(L\)-Lipschitz 连续梯度，即

则函数 \(f\) 被称为 \(L\)-光滑的。如果 \(f\) 二阶可微且 \(\mu\)-强凸，则其 Hessian 的所有特征值都大于 \(\mu\)（即 \(\nabla^2 f(x) \succeq \mu I\)），且如果它们都被 \(L\) 所界（即 \(\nabla^2 f(x) \preceq LI\)），则它是 \(L\)-光滑的。

如果其有效域 \(\{x \in \mathcal{X} : f(x) < \infty\}\) 非空，则函数 \(f\) 被称为适当的。如果其次水平集 \(\{x \in \mathcal{X} : f(x) \leq \gamma\}\) 对每个 \(\gamma \in \mathbb{R}\) 都是闭的，则它也是闭的。

非光滑分析

对于适当的、闭的凸函数 \(f\)，其次微分定义为

次微分中的元素被称为次梯度。函数 \(f\) 的 Fenchel（或凸）共轭定义为

注意 \(f^*\) 是闭凸的，无论 \(f\) 的凸性如何。\(f\) 的双共轭定义为

集合 \(C \subset \mathcal{X}\) 的指示函数定义为

适当闭凸函数 \(f: \mathcal{X} \to \mathbb{R}_{\infty}\) 的具有参数 \(\tau > 0\) 的近端算子是映射 \(\prox_{\tau f}: \mathcal{X} \to \mathcal{X}\)，定义为

由于 \(f\) 的凸性（因此目标的强凸性），最小化点总是存在且唯一。如果 \(\mathcal{X} = \mathbb{R}^d\) 具有欧几里得范数，则许多凸和非凸函数（例如范数和指示函数）的近端算子具有闭式表达 [CP11b, Bec17]，且对于许多其他函数也是已知的 [CCCP]。这类函数被称为近端友好的。例如，非负象限 \([0, \infty)^d\) 的指示函数的近端算子具有特别简单的表达 \(\prox_{\tau \iota_{[0,\infty)^d}}(x) = \max(x, 0)\)。

A.1 当 \(AA^* = \alpha I\)（\(\alpha > 0\)）且 \(f\) 的近端算子可计算时，复合项 \(f \circ A\) 的近端算子可以显式计算。在这种情况下，\(f \circ A\) 的近端算子计算为 \[ \prox_{\tau f \circ A}(x) = x + \frac{1}{\alpha} A^*\bigl(\prox_{\tau \alpha f}(Ax) - Ax\bigr). \]

参考文献