来源论文: https://arxiv.org/abs/2606.27086v1 生成时间: Jun 26, 2026 16:13

GPU加速优越化算法（Superiorization）在约束物理问题中的应用与 SupPy 框架深度解析

0. 执行摘要 (Executive Summary)

在计算物理、医学物理以及量子化学等诸多科学与工程计算领域，高维约束优化问题（High-Dimensional Constrained Optimization Problems）的求解构成了底层计算的核心。传统的求解策略通常将问题建模为严格的约束优化问题，即在可行域 $C$ 内最小化某个目标函数 $f(x)$：

$$\min_{x \in C} f(x)$$

然而，在实际的物理应用中，这种传统方法往往面临三大严峻挑战：

  1. 约束集的不相容性（Inconsistency）：由于物理测量噪声、建模误差或过于严苛的物理指标，多个约束条件交集可能为空（即 $C = \emptyset$），导致标准优化算法（如内点法、SQP法）不收敛或直接失效。
  2. 计算复杂度的阶跃：在超高维空间中，严谨地逼近全局最优解需要消耗极其庞大的计算资源。而在许多工程场景下，物理上合理的“次优解”或“优越解”已完全满足实际需求。
  3. 软偏好的硬性建模冲突：诸如去噪、正则化等辅助目标往往属于“软偏好”（Soft Preferences），将其作为硬性约束或与主目标融合进行多目标优化，会导致调参过程极其繁琐。

优越化方法（Superiorization Method, SM） 作为一种介于可行性寻优（Feasibility-seeking）与严格约束优化（Constrained Optimization）之间的全新计算范式，巧妙地解决了上述痛点。优越化方法的核心思想是：利用基本可行性寻优算法（Basic Algorithm）所具备的“有界扰动韧性”（Bounded Perturbation Resilience, BPR），在寻优迭代过程中交织注入目标函数的负梯度（或子梯度）扰动。 最终，算法不仅能够收敛到可行域（或其最近邻域），而且其最终收敛点的目标函数值 $f(x)$ 显著低于未经扰动的可行性寻优算法所得到的解。也就是说，它不以寻找绝对的数学全局最优点为目标，而是低成本地寻找一个“既可行又在物理上更优”的解。

为了将这一先进方法推向主流科学计算，Tobias Becher 等人开发了 SupPy——一个开源、模块化的 Python 原生工具箱。SupPy 实现了 CPU（基于 NumPy 和 SciPy）与 GPU（基于 CuPy）的双端无缝加速，支持多种投影算法（如 DROP、Kaczmarz/ART、EMR、EL、CG）与扰动策略的自由组合。本文将对 SupPy 的理论基础、算法细节、关键 Benchmark 表现进行深度技术剖析，并探讨其在量子化学电子结构计算与几何优化中的潜在应用价值。

1. 核心科学问题，理论基础，技术难点与方法细节

1.1 核心科学问题：可行性与最优化的权衡

在物理反问题（Inverse Problems）中，我们往往拥有一组观测数据 $b$ 和系统矩阵 $A$，需要重建物理实体 $x$。为了使重建结果符合物理先验（例如非负性、能量守恒、边界范围等），我们定义了一系列闭凸约束集 $\{C_i\}_{i=1}^m$。可行性寻优的目标是寻找一个点 $x^* \in C := \bigcap_{i=1}^m C_i$。

若在此基础上引入一个刻画物理品质（如图像总变分 TV、能量泛函）的目标函数 $f(x)$，传统的做法是求解约束最优化问题。然而，当 $m$ 极大、物理噪声高时，直接求解不仅极其缓慢，还可能因为 $C = \emptyset$ 导致无法得出任何物理上可用的结果。优越化方法提出的科学问题是：能否在不破坏可行性寻优算法收敛性的前提下，动态地引导迭代轨迹，使其向目标函数 $f(x)$ 更小的区域靠拢？

1.2 理论基础：有界扰动韧性（Bounded Perturbation Resilience）

优越化方法的数学合理性严格建立在**有界扰动韧性（BPR）**之上。设可行性寻优算法由算子 $T: \mathbb{R}^n \to \mathbb{R}^n$ 定义，其产生的无扰动迭代序列为：

$$x^{k+1} = T(x^k)$$

如果对于任意满足如下条件的有界扰动序列 $\{\beta_k v^k\}_{k=0}^\infty$：

  1. $\sum_{k=0}^\infty \beta_k < \infty$ （步长级数收敛）
  2. 序列 $\{v^k\}_{k=0}^\infty$ 有界

扰动迭代序列：

$$y^{k+1} = T(y^k + \beta_k v^k)$$

仍然能够像无扰动序列一样收敛到可行集 $C$（或在不相容情况下收敛到距离各约束集最近的邻域），则称该算法 $T$ 具有有界扰动韧性。这一性质为我们在寻找可行解的过程中“顺便”优化 $f(x)$ 提供了完美的理论保证。

1.3 数学方法细节

1.3.1 投影算子与松弛技术

对于线性不等式约束 $\langle a^i, x \rangle \le b_i$，定义约束集 $C_i = \{x \in \mathbb{R}^n \mid \langle a^i, x \rangle \le b_i\}$。点 $x$ 在 $C_i$ 上的正交投影算子 $P_{C_i}(x)$ 定义为：

$$P_{C_i}(x) = \begin{cases} x - \frac{\langle a^i, x \rangle - b_i}{\|a^i\|^2} a^i, & \text{if } \langle a^i, x \rangle > b_i \\ x, & \text{if } \langle a^i, x \rangle \le b_i \end{cases}$$

为了加速收敛或平滑迭代轨迹，引入松弛因子 $\lambda \in [0, 2]$ 与步长函数 $\sigma(x)$，得到推广的投影算子：

$$P_{C_i, \lambda, \sigma}(x) = x + \lambda \sigma(x) (P_{C_i}(x) - x)$$

当 $\sigma(x) > 1$ 时，此算子被称为外推投影（Extrapolation Projection）。

1.3.2 邻近度度量（Proximity Measure）

在约束集不相容（即 $C = \emptyset$）时，算法无法找到绝对可行解。此时通过半加权距离平方和来度量当前点与约束集族的邻近程度（Proximity）：

$$\mathcal{P}(x) = \frac{1}{2} \sum_{i=1}^m w_i \|P_{C_i}(x) - x\|^2$$

其中 $w_i > 0$ 且 $\sum_{i=1}^m w_i = 1$。SupPy 将其推广为更加通用的 $p$-范数形式的幂和度量：

$$\mathcal{P}_p(x) = \sum_{i=1}^m w_i (\mathcal{P}_i(x))^p$$

1.3.3 分裂可行性问题（Split Feasibility Problem, SFP）

在诸多复杂物理建模中（如放射治疗计划设计），约束条件往往存在于不同的物理空间。输入空间中的向量 $x \in C \subseteq \mathbb{R}^n$，经过线性物理算子 $A \in \mathbb{R}^{m \times n}$ 转换到输出空间 $y = Ax \in Q \subseteq \mathbb{R}^m$。分裂可行性问题的目标是寻找到一个 $x^*$ 满足：

$$x^* \in C = \bigcap_{i=1}^s C_i \quad \text{且} \quad A x^* \in Q = \bigcap_{j=1}^t Q_j$$

SupPy 采用了 Byrne 提出的经典 CQ 算法来求解此问题：

$$x^{k+1} = P_C (x^k + \gamma A^T (P_Q(A x^k) - A x^k))$$

其中 $\gamma \in (0, 2/L)$，$L$ 为矩阵 $A^T A$ 的最大特征值。为了避免昂贵的特征值计算，实际中常用 Frobenius 范数上限 $L^* \ge L$ 代替：

$$L^* := \|A\|_F^2 = \sum_{i=1}^m \sum_{j=1}^n |a_{ij}|^2$$

1.4 优越化算法伪代码剖析

如下所示为 SupPy 框架所采用的优越化通用控制流。其最精妙之处在于双层循环结构：内层循环通过 $n_{red}$ 步的“函数减小步骤”生成扰动点，仅在目标函数值切实降低（$f(z) \le f(x^{k,n})$）时才接纳扰动；外层循环则执行标准的物理可行性投影。

算法 1: 优越化控制流 (Superiorization Method)
输入: 初始点 x0, 可行性寻优算子 T, 目标函数 f, 内层扰动步数 n_red
1: k <- 0
2: x^k <- x0
3: while 未达到终止条件 do
4:     n <- 0  // 进入扰动阶段 (Perturbation Phase)
5:     x^{k,n} <- x^k
6:     while n < n_red do
7:         loop <- True
8:         j <- 0
9:         while loop do
10:            计算方向向量 v^k (通常为正规化子梯度: g_f(x^{k,n}) / ||g_f(x^{k,n})||)
11:            计算候选点: z = x^{k,n} + beta_{k,j} v^k
12:            j <- j + 1
13:            if f(z) <= f(x^{k,n}) then
14:                n <- n + 1
15:                x^{k,n} = z
16:                loop <- False
17:            end if
18:        end while
19:    end while
20:    x^{k+1} = T(x^{k,n})  // 进入可行性寻优步骤 (Feasibility-Seeking Phase)
21:    k <- k + 1
22: end while

1.5 技术难点：步长衰减、重启动（Restart）与自适应机制

要想使上述扰动算法成功收敛，扰动的步长 $\beta_k$ 必须严格控制。SupPy 实现了两种核心步长控制策略：

  1. 幂律衰减与重启动（Power Law with Restarts）：

    $$\beta_k = \gamma \alpha^\ell$$

    其中 $\gamma$ 为基础步长核，$\alpha \in (0, 1)$ 为衰减因子，$\ell$ 为随迭代次数递增的整数。若 $\beta_k$ 衰减过快，扰动在后期将完全失效。因此，SupPy 引入了重启动技术（Restart）：每隔一定步数（如50次迭代），重置 $\ell$，但为了保持 BPR 的数学收敛约束，每次重启动后 $\ell$ 的起点逐步增加（例如首次重启动后起点设为1，第二次设为2，以此类推），确保整体级数 $\sum \beta_k$ 收敛。

  2. 自适应子梯度投影（Adaptive Subgradient Projection）： Nikazad 等人提出的自适应选择策略，无需手动设定繁琐的衰减参数：

    $$\beta_k := \begin{cases} 0, & \text{if } f(x^k) \le f_{ref} \\ \frac{f(x^k) - f_{ref}}{\|g_f(x^k)\|}, & \text{if } f(x^k) > f_{ref} \end{cases}$$

    其中 $f_{ref}$ 为动态更新的目标函数参考值，$g_f(x^k)$ 为 $x^k$ 处的子梯度。该方法既能保证有界扰动韧性，又能在目标函数显著高于参考值时提供更强烈的扰动纠偏。

2. 关键 Benchmark 体系、计算所得数据与性能展示

SupPy 在三个截然不同的物理反问题场景上进行了严格的基准测试，充分展示了其异构加速与物理重构的优越性能。

2.1 Benchmark 1: 二维地震波层析成像 (Seismic Tomography)

2.1.1 物理建模背景

重建穿越地质介质的地震波传播速度分布。该问题被建模为大型稀疏线性方程组：

$$A x \approx b$$

其中 $A \in \mathbb{R}^{8192 \times 65536}$（由 AIR Tools II 中的 seismicwavetomo 生成，模拟 64 个震源与 128 个接收器，分辨率为 $256 \times 256$ 像素），$b$ 为实测的地震波传播时间。测试集分为两组：无噪 clean 数据，以及带有 5% 高斯噪声的 noisy 数据。重构算法采用 DROP（对角松弛正交投影），扰动策略包含 L1 正则、L2 正则以及全变分（TV）正则。

2.1.2 关键性能数据与物理重建结果

  • 计算速度提升：在相同测试机上，使用传统的 MATLAB CPU 版 AIR Tools II 运行 1000 次 DROP 迭代需耗时约 100 秒；而使用 SupPy 的 GPU（CUDA）加速版本，计算时间骤降至 10 秒以内，实现了 10倍以上 的绝对加速。
  • 重构误差分析（以相对误差 $\epsilon = \|x^k - x_T\|/\|x_T\|$ 评估，其中 $x_T$ 为地质真实模型）：
    • 无噪数据（Clean Data）：如图 1（左）所示，随着迭代进行，无扰动可行性寻优、L1-优越化以及 TV-优越化的相对误差均平滑下降。经过 1000 次迭代后，无扰动（DROP）的误差为 0.1510，TV-优越化降至 0.1504。L2 优越化在 clean 数据中表现不佳，1000次迭代后误差滞留在 0.2001 左右。
    • 多噪数据（Noisy Data）：如图 1（右）所示，噪声引起了明显的过拟合现象。未经优越化的 DROP 算法在第 300 次迭代附近达到误差极小值，随后相对误差反弹（由于过拟合噪声）。而 TV-优越化 和 L1-优越化 则展现出极强的去噪能力，在限制物理噪声的同时，相对误差持续保持在更低水平（TV-优越化极小误差为 0.1943，而标准 DROP 极小误差为 0.1972）。L2-优越化受噪声干扰严重，误差保持在 0.2119 处。这直观证明了非平滑正则化（如 TV 和 L1）在优越化架构中对反问题求解的巨大价值。

2.2 Benchmark 2: 低剂量 CT 图像重建 (Low-Dose CT Reconstruction)

2.2.1 物理建模背景

采用标准的 LoDoPaB-CT 数据集，图像分辨率为 $362 \times 362$。系统矩阵 $A$ 维度高达 $(513000) \times (131044)$。此基准测试对比了四种不同的可行性寻优算法：

  1. ART (Algebraic Reconstruction Technique, 顺序 Kaczmarz 算法)
  2. EL (Extrapolated Landweber, 外推 Landweber 算法)
  3. EMR (Error Minimizing Relaxation, 最小化误差松弛算法)
  4. CG (Conjugate Gradient, 共轭梯度法)

评估了这四种算法在“仅进行可行性寻优”与“实施 TV 优越化”两种模式下的表现，同时考察了 GPU 对大规模稀疏矩阵乘法的加速效果。

2.2.2 性能与数据定量表征

下表汇总了 CT 图像重建的关键数据（提取自论文 Table 4）：

算法类型 (Alg.)硬件端300次迭代 寻优耗时 $t_{\text{feas}}$ (s)300次迭代 寻优极小误差 $\epsilon_{\min}$300次迭代 优越化耗时 $t_{\text{sup}}$ (s)300次迭代 优越化极小误差 $\epsilon_{\min}$自适应方差终止 寻优迭代步数 $n_{\text{feas}}$自适应方差终止 优越化迭代步数 $n_{\text{sup}}$自适应方差终止 优越化最终误差 $\epsilon_{\text{stop}}$
ARTCPU14471.015781.014714411.392
EMRGPU CPU10.0 2350.09317.3 3480.0661091090.075
ELGPU CPU7.2 1820.32115.1 2770.316——0.661
CGGPU CPU9.4 2290.09015.5 0.0680.0681151140.069

2.2.3 数据深度解读

  1. GPU 硬件级加速表现：对于 EMR 算法，单纯可行性寻优在 CPU 上需要 235秒，而 GPU 上仅需 10.0秒，加速比达 23.5倍；对于优越化模式，CPU 耗时 348秒，GPU 仅需 17.3秒，加速比达 20.1倍。这充分体现了 CuPy 库对海量投影运算的强大硬件压榨能力。
  2. 优越化带来的图像质量飞跃：以 EMR 为例，300次迭代下，不加优越化的重建误差为 0.093，加入 TV-优越化后相对误差大幅降至 0.066。如图 5、图 6 所示，TV 优越化不仅大幅降低了图像中的条纹和散斑噪声，还极其完美地保留了生物软组织的解剖学边缘，极大地提高了低剂量条件下的重构品质。
  3. 算法健壮性与不相容处理：EL 和顺序 ART 在不相容系统上表现极差，ART 在 300 步时的误差高达 1.0。相比之下，EMR 和 CG 算法表现出极高稳定性，并且配合自适应方差终止标准（Variance stopping criterion），能在不明显牺牲重构质量的前提下提早终止迭代（如 EMR 在第 109 步即成功终止，耗时缩短至 GPU 端的 6.7秒，误差仅为 0.075）。

2.3 Benchmark 3: 放射治疗计划设计 (Radiotherapy Treatment Planning)

2.3.1 物理建模背景

在强约束物理设计中，我们需要决定 9 个方向的入射光子束流强度 $w \in \mathbb{R}^n_+$，使得在复杂人体组织内的剂量分布满足医生设定的临床硬性指标（如肿瘤靶区 Target 的处方剂量，以及对核心危及器官 Core、Surrounding Body 的剂量上限约束）。该问题极其复杂，由于物理剂量学的天然不相容性，这组约束集合通常为空集（$C = \emptyset$）。传统的非线性规划求解器（如基于内点法的 IPOPT）在这种不相容约束下往往无法收敛到可行解甚至彻底罢工。

测试靶区分别采用了马蹄形物理模型（Horseshoe-shaped phantom）以及实际的头颈部患者临床数据（Head-and-neck patient），基于 CORT 公开数据集。目标函数为最小化正常机体组织（Body）的平均剂量 $f_{Mean}$。约束指标如表 2 所示：

  • 靶区约束：$D_{95\%} \ge 2\text{ Gy}$，$D_{5\%} \le 2.17\text{ Gy}$，且满足剂量区间限制 $1.93\text{ Gy} \le d \le 2.27\text{ Gy}$；
  • 核心保护区约束：$D_{5\%} \le 0.33\text{ Gy}$。

2.3.2 性能与物理剂量分布评估

  • 计算鲁棒性与收敛速度：
    • 对于马蹄形模型，IPOPT 求解器由于硬性约束冲突，在迭代中挣扎并宣告失败，无法给出临床可接受计划。而 SupPy 优越化方法在 GPU 上迭代 10000 步（仅耗时 300秒，而纯可行性寻优耗时 105秒）即给出了物理上完美契合的临床剂量分布。
    • 如图 8 所示，通过 DVH（剂量-体积直方图）对比，优越化计划（实线）在完美包覆肿瘤靶区（Target）的同时，使周围健康机体（BODY）的受照剂量曲线显著向左下方平移（表示更低的辐射损伤）。具体而言，相比于单纯的可行性寻优，优越化算法成功在不妥协任何硬性肿瘤剂量的指标下，将机体（BODY）的健康受照剂量降低了近 15% 到 20%。这直接印证了优越化算法在处理不相容物理设计问题时的绝对统治力。

3. 代码实现细节、复现指南与开源生态

3.1 SupPy 架构设计细节

SupPy 采用了极富远见的解耦设计（Decoupling Design）。它将可行性寻优模块（projections）与扰动算法模块（perturbations）完全独立。这种高内聚、低耦合的设计使得物理学家能够自由组装、测试全新的算法。其类层次继承关系如下：

  • BaseProjection：作为所有投影算子的虚基类，向高层框架统一提供 project(x)、proximity(x) 等关键接口。
  • CustomProjection：允许用户通过简单的 lambda 函数定制任意复杂的非线性投影操作。
  • Superiorization：主控制类，接收任意 BaseProjection 子类实例以及 BasePerturbation 策略实例，驱动核心迭代控制流。

3.2 典型代码实例深度剖析

3.2.1 简单双球体可行性寻优 (Listing 1 剖析)

import numpy as np
from suppy.projections import BallProjection, SequentialProjection

# 1. 定义两个相互交叠的物理球体约束集，分别设定球心坐标与半径
center_1, center_2 = np.array([1.0, 1.0]), np.array([-0.5, -1.0])
radius_1, radius_2 = 2.0, 1.0

# 2. 实例化球体投影算子
ball_1 = BallProjection(center_1, radius_1)
ball_2 = BallProjection(center_2, radius_2)

# 3. 组合为顺序投影算法 (Sequential Projection)
joined_projection = SequentialProjection([ball_1, ball_2])

3.2.2 高阶物理优越化重构实现 (Listing 2 剖析)

以下代码展示了如何快速构建一个优越化求解器。其目标是：在寻找上述两球体交集的可行空间内，使系统向能量极低值点（即原点，目标函数 $f(x) = x \cdot x$）靠拢。

import numpy as np
from suppy.projections import BallProjection, SequentialProjection
from suppy.perturbations import PowerSeriesGradientPerturbation
from suppy.superiorization import Superiorization

# 1. 定义物理能级目标函数 (f(x) = x^T x) 及其解析梯度 (g_f(x) = 2x)
func_1 = lambda x: x @ x
grad_1 = lambda x: 2.0 * x

# 2. 几何约束模型构建
center_1, center_2 = np.array([1.2, 0.0]), np.array([0.0, 1.4])
radius = 1.0
ball_1 = BallProjection(center_1, radius)
ball_2 = BallProjection(center_2, radius)
proj = SequentialProjection([ball_1, ball_2])

# 3. 实例化带幂律衰减的有界梯度扰动器
pert = PowerSeriesGradientPerturbation(func_1, grad_1)

# 4. 组装优越化核心引擎
sup = Superiorization(proj, pert)

# 5. 执行对比测试
x0 = np.array([2.5, 1.5])
# 运行纯物理投影重构
x_proj = proj.solve(x0, storage=True)
# 运行优越化重构 (在10步迭代内注入目标函数扰动)
x_sup = sup.solve(x0, 10, storage=True)

3.3 无缝异构硬件切换原理

SupPy 的一大设计亮点是零感知硬件切换。其底层并不显式暴露复杂的 CUDA 内存拷贝代码，而是通过对输入数据类型的动态分发（Dynamic Dispatching）实现。例如：

  • 如果传入 x0 的类型为 numpy.ndarray，框架内部会自动调用 NumPy 及 SciPy 稀疏矩阵加速运算，所有任务运行在 CPU 核心上。
  • 如果传入 x0 的类型为 cupy.ndarray，框架检测到 GPU 统一内存指针后，会将所有中间状态矩阵运算无缝切换到 CuPy，使数十万个线程在 GPU Stream 上异步执行。

3.4 快速复现指南与开源链接

  • 官方源码库：项目托管于 GitHub，遵循开源 BSD 3-Clause 协议。点击访问 SupPy 官方 GitHub 仓库。
  • 版本说明：论文中的基准数据均基于 v0.4.0 专用分支运行。其 Zenodo 归档 DOI 为 10.5281/zenodo.20085222。
  • 安装步骤：
    # 1. 创建干净的物理科学计算虚拟环境
conda create -n suppy_env python=3.10 -y
conda activate suppy_env

# 2. 安装 NumPy, SciPy (针对CPU端)
pip install numpy scipy matplotlib

# 3. 依据您的 CUDA 版本安装 CuPy (以 CUDA 12.x 为例)
pip install cupy-cuda12x

# 4. 克隆并安装 SupPy 库
git clone https://github.com/DKFZ-OpenMedPhys/SupPy.git
cd SupPy
pip install -e .

4. 关键引用文献与学术局限性批判

4.1 关键引用文献

在优越化算法和计算物理领域，本研究立足于以下里程碑式著作：

  1. [1] Censor, Y., Davidi, R., Herman, G. T. Perturbation resilience and superiorization of iterative algorithms. Inverse Problems, 2010. （奠定了优越化算法数学完备性与有界扰动韧性的理论根基）
  2. [17] Censor, Y., Elfving, T., Herman, G. T., Nikazad, T. On Diagonally Relaxed Orthogonal Projection Methods. SIAM Journal on Scientific Computing, 2008. （提出了著名的 DROP 算法，是地震重构中可行性寻优的核心）
  3. [18] Nikazad, T., Abbasi, M., Elfving, T. Error minimizing relaxation strategies in Landweber and Kaczmarz type iterations. Journal of Inverse and Ill-posed Problems, 2017. （提出了高鲁棒性的 EMR 算法）
  4. [30] Byrne, C. L. Iterative oblique projection onto convex sets and the split feasibility problem. Inverse Problems, 2002. （奠定了多空间分裂可行性求解的数学框架）
  5. [29] Aragón-Artacho, F. J., Censor, Y., Gibali, A. The superiorization method with restarted perturbations for split minimization problems. Applied Mathematics and Computation, 2023. （引入了重启动扰动机制，成功解决了长周期迭代中步长过度衰减的问题）

4.2 SupPy 框架与优越化方法的局限性批判

尽管 SupPy 框架及优越化算法在处理物理反问题时展现出惊人的效果，但站在严谨的学术与量子化学研究视角，该方法依然存在若干急需攻克的底层技术及理论缺陷：

4.2.1 物理约束的高阶慢收敛陷阱 (DVC Bottleneck)

在放射治疗计划（IMRT）中，剂量-体积约束（DVC） 具有高度的非凸与非线性特征。SupPy 目前在处理 SFP 问题时，其对于 DVC 约束的投影收敛极其缓慢。正如论文所述，使用 DVC 约束优越化在 GPU 上生成一个马蹄形计划需长达 4.5分钟。当我们将 DVC 约束替换为简单的线性最大剂量约束（Linear max-dose constraints）时，在几乎完全相同的重构剂量质量下，收敛时间暴降至 50秒。这说明当前的投影框架在面对高阶复杂的非平滑、非凸约束集合时，其投影流的速度受到了极大的限制。

4.2.2 超参数寻优的“炼丹化”窘境

优越化算法的性能高度依赖于超参数的设定，包括步长衰减核 $\gamma$、衰减因子 $\alpha$、单步内层扰动迭代数 $n_{red}$ 以及重启动间隔 $n_{restart}$。虽然论文中设计了自适应策略，但不同物理场景的数学泛函形态迥异，自适应参考值 $f_{ref}$ 的衰减速度无法被严格解析证明。这导致在探索新物理问题时，研究人员必须耗费大量精力进行超参数网格搜索。这对于亟需高确定性、高精度的科学计算而言，是一项不可忽视的实用壁垒。

4.2.3 严格数学收敛阶（Convergence Rate）的缺失

需要澄清的是：优越化方法的理论仅证明了**“收敛到可行集”以及“最终解的目标函数值不差于纯可行性解”**。至今，学术界仍未给出一个严谨的、能够证明优越化解与绝对全局最优点（Global Optimum）之间数学代数距离的紧致边界（Tight Bound）。这种“保障可行，顺便优化”的哲学，虽然在医学物理和地震去噪中绰绰有余，但在需要绝对最低能量极值点、极度敏感的量子化学（如过渡态能量计算）中，很难让人放心地彻底取代严格的约束优化求解器。

4.2.4 自定义算子的 CUDA 级性能损失

虽然 SupPy 依赖 CuPy 实现了底层数组级的 GPU 计算，但是当用户派生 CustomProjection 来定义极其复杂的非线性物理投影时，CuPy 无法将其自动融合（Kernel Fusion）为一个单一的 CUDA Kernel。大量的中间数组创建以及内存搬运会导致频繁的 GPU 寄存器与显存交互，这极大地拖累了运算速度。目前该工具箱亟需引入即时编译技术（如 Numba CUDA 或 CuPy RawKernel）以真正释放极致的高性能计算能力。

5. 其他必要补充：优越化方法在量子化学与电子结构计算中的潜在应用与理论拓展

尽管本文的研究背景聚焦于医学物理与地球物理成像，但其底层数学内核——“高维空間约束寻优 + 能量正则泛函削减”——与量子化学（Quantum Chemistry, QC）及电子结构理论（Electronic Structure Theory） 存在天然且深邃的契合性。本节专门为量子化学科研工作者，探讨优越化方法（SM）在现代计算化学领域的四大潜在应用方向，并推演其理论公式。

5.1 在约束密度泛函理论（Constrained DFT, CDFT）中的应用

5.1.1 物理问题背景

约束密度泛函理论（CDFT）在研究分子内电荷转移、非绝热耦合、极化子激发态等领域应用广泛。其数学目标是：在严格限定某些特定空间区域（例如给体 Donor 区域 $D$ 与受体 Acceptor 区域 $A$）电荷差（或自旋多重度）的约束下，寻找使电子结构系统自由能最低的电子密度 $\rho(r)$。传统的 CDFT 求解使用拉格朗日乘子法（Lagrange Multipliers）：

$$\mathcal{W}[\rho, V_c] = E_{DFT}[\rho] + V_c \left( \int_D \rho(r) dr - \int_A \rho(r) dr - N_c \right)$$

通过对势能 $V_c$ 进行不断的牛顿法迭代来调整势垒。当分子体系庞大且电子能级极其接近时，该多重外循环经常陷入难以收敛的境地（SCF 振荡）。

5.1.2 优越化求解方案设计

我们可以将电子结构体系的可行性集合与能量优化完全拆解：

  1. 约束集合定义（Feasibility-seeking）： 定义凸约束集 $C$，其要求三维空间电荷差严格等于处方值 $N_c$：

    $$C = \left\{ \rho(r) \in \mathbb{L}^1(\mathbb{R}^3) \ \middle| \ \int_D \rho(r) dr - \int_A \rho(r) dr = N_c, \ \int \rho(r) dr = N_{total}, \ \rho(r) \ge 0 \right\}$$

    对电荷约束集的投影 $P_C$ 可以极快地通过在实空间网格（Real-space Grids）上对电荷进行空间线性平移与正规化（Scaling）来实现，计算成本极低，且永远不会产生发散问题。

  2. 能量优越化（Superiorization）： 将目标函数 $f(\rho)$ 设为原生的密度泛函能量 $E_{DFT}[\rho]$。其负子梯度正好对应于系统的一阶单粒子哈密顿量有效势：

    $$-g_f(\rho) = - \frac{\delta E_{DFT}}{\delta \rho(r)} = - v_{eff}(r)$$

    在迭代过程中，我们无需去精确求解繁琐的 $V_c$ 拉格朗日势。我们只需要先对当前电荷密度进行实空间势场的步长扰动：

    $$\rho_{perturbed} = \rho^k - \beta_k \frac{v_{eff}(r)}{\|v_{eff}\|}$$

    然后将这组扰动后的电荷密度投影回满足严格电荷差约束的可行集 $C$ 中。这不仅在数学上保证了电荷分布永远合法（严格满足电荷差），而且能以极低的计算代价将系统平稳引导至最低能量状态。这彻底规避了拉格朗日势振荡发散的世纪难题。

5.2 轨道正交化约束与自洽场（SCF）快速收敛

5.2.1 物理问题背景

在 Hartree-Fock (HF) 或 DFT 迭代求解 Kohn-Sham 方程中，单粒子轨道系数矩阵 $C$ 必须在每一步都严格满足空间自正交条件：

$$C^T S C = I$$

其中 $S$ 为重叠积分矩阵。传统的自洽场（SCF）方法使用对角化哈密顿量矩阵（F diagonalization）来隐式满足此条件。但在过渡金属复合物、大体系激发态及低能隙系统中，SCF 的对角化极易出现能级交叉、反复振荡甚至不收敛现象。

5.2.2 优越化求解方案设计

利用优越化方法，我们可以直接抛弃全对角化（Diagonalization-free）：

  • 可行性寻优步骤：定义可行集 $C_{ortho} = \{C \mid C^T S C = I\}$。点在 $C_{ortho}$ 上的投影就是经典的 Gram-Schmidt 正交化、对称正交化（Löwdin Orthogonalization）或奇异值分解（SVD）投影，这些几何投影在数值上极度稳定且计算代价极低。
  • 优越化扰动步骤：设定目标函数为 Hartree-Fock 能量泛函 $E_{HF}(C)$。其负梯度方向对应于 Fock 矩阵施加在当前系数上的梯度： $$-g_f(C) = - [F(C), C]$$ 在优越化架构下，我们将当前系数矩阵沿 Fock 梯度方向平移一小步（减小电子能级），随后直接通过 Gram-Schmidt 投影拉回正交可行域： $$C^{k+1} = P_{ortho} \left( C^k - \beta_k [F(C^k), C^k] \right)$$ 这种方法不涉及昂贵的 Fock 矩阵完整对角化，通过局部梯度与几何正交化投影的交织，能够极其平滑地将电子轨道推向基态，特别适合超大规模体系的线性标度（$O(N)$）求解。

5.3 固体能带结构设计中的禁带宽度最大化控制

在计算材料化学中，设计具有特定带隙（Band Gap）的晶体结构是一项重要任务。设晶胞中的势场排布为 $V(r)$，约束集规定了势场在物理空间上的平滑度、总积分上限和对称性：

$$C_{V} = \left\{ V(r) \ \middle| \ \int |\nabla V|^2 dr \le μ, \ \int V dr = V_0 \right\}$$

我们的目标是最大化带隙 $E_g(V) = E_{LUMO}(V) - E_{HOMO}(V)$（等价于最小化 $-E_g(V)$）。利用 SupPy 架构，研究人员可以把带隙函数设为优越化目标 $f(V) = -E_g(V)$，晶胞势场的空间物理约束作为可行集。在迭代重构晶格势场时，一方面确保势场的物理合法性（通过投影满足 $C_V$），另一方面通过带隙梯度动态微调势场：

$$V^{k+1} = P_{C_V} \left( V^k + \beta_k \frac{\nabla_V E_g(V^k)}{\|\nabla_V E_g(V^k)\|} \right)$$

这种新颖的晶体设计范式，能够极具效率地从无到有“雕刻”出符合目标物理性质的新型光电半导体材料。

5.4 几何优化（Geometry Optimization）中的硬性键长/键角约束

在研究多肽大分子的构象寻找、分子对接（Molecular Docking）或者特定配位复合物的几何优化时，我们往往需要保持某些关键化学键（如配体环、刚性芳香环）的键长、键角或二面角严格不变（例如冻结某些内坐标），而在剩余自由度上最小化分子势能。这种带约束的分子动力学/几何优化，在笛卡尔坐标系下通常需要借助极其复杂的拉格朗日乘子算法（如 SHAKE、RATTLE 算法）。

利用 SupPy 工具箱的数学思想：

  • 将键长/键角刚性不变的几何限制视为凸（或局部凸）投影集。例如，限制原子 $i, j$ 间距为 $d_0$，其正交投影仅仅是沿着它们的连线方向拉伸或缩短，使其距离等于 $d_0$，投影计算成本为零；
  • 分子力场下的总势能（Potential Energy）设为目标函数 $f(X)$，其负梯度 $-g_f(X)$ 对应于经典的原子受力（Atomic Forces）。

在每一步几何优化中，原子首先沿着牛顿力场给出的受力方向移动一小步（能量衰减），随后直接通过快速空间几何投影将各原子拉回到满足刚性化学键的限制位置上。这种方法的计算复杂度仅与原子数呈线性关系，且永远不会发生 SHAKE 算法在强力场扰动下因不收敛而导致的系统崩溃（Explosion）。这为大规模高分子材料与生物物理大分子的经典/量子力学混杂模拟（QM/MM）提供了全新的高效几何优化思路。

6. 总结

Tobias Becher 等人开发的 SupPy 框架，不仅为医学图像重建和地球物理反问题提供了一款高效率、高拓展性的 GPU 异构加速利器，更向计算物理与化学界展示了优越化算法（Superiorization）在应对高度复杂、不相容约束及高维寻优时的独特魅力。虽然其在非凸、高度非平滑约束下的数学理论和收敛性能仍有待深入挖掘，但对于那些不执着于绝对数学全局最优点、追求在物理可行性和高效计算之间达到极佳平衡的科学研究而言，优越化方法无疑开辟了一条全新的康庄大道。