带约束优化问题与增广拉格朗日方法

引言

我们的目标还是为了解决 病态的逆问题。为了解决此类问题，我们需要引入一些约束和正则化手段，以获得有意义的解。这里的约束条件包括：

1. 正值性约束(Non-negativity constraint)：要求求得的解 x 的所有元素均为非负数，即 x_i \geq 0。
2. 支持性约束(Support constraint)：要求信号 x 在特定的范围（时间域或空间域）内为非零，在其他区域则应为零或接近零。

为实现这些约束，我们通常采用以下方法：

**增广拉格朗日方法 (Augmented Lagrangian) ** + 分布式交替方向乘子法 (ADMM)

在有约束的优化问题中，为了在不破坏问题可解性和结构性的前提下处理约束条件，常用增广拉格朗日方法 (Augmented Lagrangian)。其思想是：

• 原始优化问题通常可表示为最小化目标函数 J(x)，同时满足约束条件 g(x)=0 或 h(x) \leq 0。

• 增广拉格朗日方法将约束条件通过拉格朗日乘子(Lagrange multipliers)嵌入目标函数中，同时添加二次惩罚项，使得当约束不满足时会导致目标函数有额外惩罚，鼓励算法在迭代过程中逐步满足约束。

• 最终，我们得到一个新的增广拉格朗日函数 L(x, \lambda)，其中 \lambda 是拉格朗日乘子。通过交替更新 x 和 \lambda，可以逼近原问题的解。

ADMM (交替方向乘子法) 是一种适合大规模并行或分布式计算的优化求解方法。它的原理是在增广拉格朗日框架下，将原问题分解为多个子问题，每个子问题仅涉及一部分变量，从而更易求解。具体步骤是：

• 将原问题中的变量分裂为两个（或多个）组，比如 x 和 z。

• 在每次迭代中，先固定一个组的变量，求解关于另一个组变量的子问题；然后固定另一个组变量，求解关于前一个组的子问题。

• 在每个子问题求解完成后，更新拉格朗日乘子以确保最终解满足约束。

ADMM 结合 拉格朗日乘子法 和 增广拉格朗日法，能够分解大规模优化问题为易处理的子问题，通过交替更新不同变量，逐步逼近最优解。

梯度方法和分裂技术

梯度方法依赖目标函数的梯度信息来决定搜索方向，从而逐步逼近最优解。对于二次可微的目标函数，梯度方法（如最速下降法、共轭梯度法）能够高效求解。对于更复杂、更难处理的函数，梯度方法也可以与其它技巧相结合。

分裂技术 (Splitting) 将复杂的优化问题分解为多个子问题。通过将约束或正则化项分别分配给不同的子问题，我们在每一步解决较为简化的子问题，然后合并这些子问题的解来更新整体解。它与 ADMM 本身的思想相同，因此经常配合使用。

卷积/去卷积问题的模型描述

在处理卷积/去卷积问题时，模型描述为：

其中 y 是观测值，H 是卷积算子，x 是原始信号，e 是噪声，\ast 表示卷积运算。目标是通过去卷积和降噪恢复信号 x，即 \hat{x}(y)。

因为 H 通常是已知的（例如已知系统的冲激响应 h），理论上我们希望从 y 中“反演”出 x，这称为反卷积或逆问题。然而，由于卷积本身容易导致频域成分的混合且可能存在较小奇异值（特别是当 H 不可逆或近似不可逆时），使得问题病态。小的噪声 e 都会在解 x 中被放大。

为了解决这一病态问题，我们引入正则化项限制解的形状和特性，从而得到更稳定的解。

带正则化的优化问题

在给定噪声数据 y 的情形下，我们的目标是找到一个合适的 x 来最小化以下代价函数：

• 参数 \mu 控制数据拟合项和正则化项之间的权衡
• D 为梯度算子，用于捕捉信号变化的梯度信息。

通过加上一项 \mu \| D x \|^2 的惩罚，可以鼓励 x 的平滑性。这在反卷积中至关重要，因为病态的反问题往往会将小噪声放大成剧烈的高频振荡。

信号恢复的解

在信号恢复的优化问题中，我们定义了代价函数：

目标是求解：

通过最优化求导的方式，可以得到问题的正常方程（normal equation）。具体过程如下：

将 J_{\text{PLS}}(x) 展开并对 x 求梯度（设为零以求极值点）：

对 x 求梯度：

整理得：

将 (H^T H + \mu D^T D) 拼起来，那么：

在这个方程中，只要 (H^T H + \mu D^T D) 是可逆的，就可以直接求得闭式解：

此时，我们已经给出了一个明确的最小二乘结合正则化的求解公式，但对高维数据而言，直接求逆 (H^T H + \mu D^T D) 算不了。

基于快速傅里叶变换 (FFT) 的高效求解

在许多实际应用中，卷积算子 H 以及梯度算子 D 都可以用循环结构（circulant structure）近似。循环矩阵有一个重要性质：它们可以通过快速傅里叶变换（FFT）对角化。将卷积问题转到频域后，我们可以利用 FFT 将卷积变为点积，从而将矩阵求逆问题转化为对角元素的求解问题。

对 H 和 D 做循环近似，将其表示为对应的循环矩阵。在频域中 H 与 D 分别对应的对角矩阵为 \Lambda_h 与 \Lambda_d，它们的元素为：

• 这里 \overset{\circ}{h}_n 和 \overset{\circ}{d}_n 分别是 H 和 D 在频域对应频率点 n 的值。

在频域下，对应的代价函数的解可以表示为：

因为 (\Lambda_h^\dagger \Lambda_h + \mu \Lambda_d^\dagger \Lambda_d) 是对角矩阵，其逆矩阵非常容易求得。对角矩阵的求逆只需要对每个对角元素求逆，并不需要复杂的矩阵分解。于是：

得到频域解 \overset{\circ}{\hat{x}} 后，通过逆FFT 即可回到时域（或空间域）得到 \hat{x}_{\text{PLS}}。

我们在反问题实际运用中已经详细解释过此过程，其具体数学表达解释和应用参考 Winer - Hunt 方法理论及实现

此方法的核心是将原本难以直接求逆的高维矩阵，通过频域的对角化处理来简化矩阵求解过程，这个方法是最小二乘和正则化相结合的典型，同时利用 FFT 技术优化了计算性能。

不同求解方法的比较

我们上面用循环近似矩阵和傅立叶基来逃避了 (\Lambda_h^\dagger \Lambda_h + \mu \Lambda_d^\dagger \Lambda_d)^{-1} 这个矩阵求逆的问题，当然还有一些其他方法。

1. 莽夫流派，直接计算

   明知山有虎偏向虎山行，低维问题没问题，高维问题不可能解的了

2. 线性系统求解算法

   • 高斯消元法（Gauss, Gauss-Jordan）：逐步消去未知数，适合小规模问题。

   • 替代法（Substitution）：利用已知变量逐步消去未知数。

   • 三角分解法（Triangularisation）：将矩阵分解为上下三角矩阵，以加速求解步骤。适合较大但仍需计算资源有限的问题。

3. 数值优化

   • 梯度下降及其变种（如共轭梯度法）：不必显式求逆，可以利用梯度信息逐步逼近最优解，适用于目标函数可导的优化问题。
   • 逐像素更新（Pixel-wise, Pixel by Pixel）：对图像类问题中特定像素逐步优化，适合高维度且稀疏问题的场景。

4. 对角化方法

   如本文所述，将矩阵转化为循环矩阵形式简化求解过程，再通过快速傅里叶变换（FFT）实现对角化，进而简化求逆。

5. 特殊算法（尤其针对一维情况）

   • 递归最小二乘法（Recursive Least Squares）：动态调整参数，适合实时信号处理问题。
   • 卡尔曼平滑器或滤波器（Kalman Smoother or Filter）：提供动态系统状态的最优估计，可随时间递归更新。

约束条件的具体表达

• 正值性约束：

  C_p : \forall p \in \mathcal{M}, \ x_p \geq 0

  对应到矩阵形式时，如果要将其纳入增广拉格朗日或ADMM框架，可以表示为 x \geq 0，这可以看做是一个非负约束的凸集。

• 支持性约束（掩模）：

  C_s : \forall p \in \bar{\mathcal{S}}, \ x_p = 0

  这里 \mathcal{S} 是信号的支持区域（例如已知信号只存在于图像的某个子块中），\bar{\mathcal{S}} 则是其补集。在掩模外强制 x_p=0 可以视为在某些维度上严格约束解。

优化问题的约束最小化解

约束最小化问题的总体框架

优化目标函数为：

线性约束条件为：

• 支持区域外像素值强制为零：

  x_p = 0, \quad \text{for } p \in \bar{\mathcal{S}}

• 像素值必须为非负：

  x_p \geq 0, \quad \text{for } p \in \mathcal{M}

凸性质：

目标函数的凸性：

这是一个二次型目标函数，且通常矩阵 H^T H 与 D^T D 都是半正定的。当 \mu > 0，J_{\text{PLS}}(x) 严格凸，从而保证唯一解的存在。

约束集合的凸性：

• 支持性约束 x_p = 0 \, (\forall p \in \bar{\mathcal{S}}) 对应线性等式约束，是一个凸集合。
• 非负性约束 x_p \geq 0 对应正正交(positive orthant)，也是一个凸集合。

由于目标函数和约束集合都是凸的，并且目标函数是严格凸的，所以优化问题存在唯一解。

我们综合上述三条内容，可得:

正值性约束的单变量与多变量分析

正值性约束的单变量分析

先考虑最简单的情况，即单变量约束问题：

• 无约束解：直接取 t = \bar{t}，这是目标函数的最小值点。

• 有约束解：

  • 若 \bar{t} \geq 0，则约束不活跃，解仍为 \hat{t} = \bar{t}。
  • 若 \bar{t} < 0，则为满足 t \geq 0，解被修正为 \hat{t} = 0。
  • 最终形式为：
  \hat{t} = \max[0, \bar{t}]

正值性约束的双变量分析

考虑有两个变量的情形：

无约束时，解可以通过求导得到。但是当约束加入后，如果无约束解出现 t_1<0 或 t_2<0 的情况，就需要将其投影回正半空间 t_1 \geq 0,t_2 \geq 0。与单变量情形类似，这是一种多维投影（项目）操作。实际求解时可使用坐标下降、投影梯度等方法逐步修正解，使之满足正值性约束。

带等式约束优化问题的求解

我们之前写过优化问题定义为：

即在我们考虑的优化问题中，有一部分像素的取值是被严格约束的：这些像素必须为零（例如在某些区域外不允许有信号）。为了处理这个约束，我们想办法把问题简化成一个不包含这些受约束像素的子问题，然后再回到原始问题中构建最终解，这就是选择矩阵 T的作用。称之为截断。

截断操作

引入一个选择矩阵 T，用于选取未被约束的像素：

• T 是一个 M \times N 的矩阵，M < N，这意味着在所有 N 个像素中，我们只选择了 M 个像素作为自由变量（其他固定为零的像素忽略）。
• 当我们用 T 作用在 x 上时，T x 会返回一个只有 M 个分量的向量，表示这些「未被约束」像素的值。
• 同理，如果我们有一个长度为 M 的向量 \bar{x}，通过 x = T^t \bar{x} 就能把它「扩展回」原本长度为 N 的空间，在没有被选择的位置填上0（这正好满足了那些等式约束：要求某些位置的像素为零）。

通过选择矩阵 T，原本的优化问题就变成了一个较小的「无约束」问题：我们只在较小维度的空间里找那 M 个未被约束像素的最优解。等求出 \bar{x} 后，再用 x = T^t \bar{x} 把解映射回原始空间，就能得到满足约束的全局解。

闭式解的推导

1. 原始问题：

   \hat{x} = \arg \min_x \| y - Hx \|^2 + \mu \| D x \|^2

   受限于：某些坐标必须为0。

2. 用 x = T^t \bar{x} 替换后，我们的优化目标就只针对 \bar{x} 来写：

   \bar{J}_{\text{PLS}}(\bar{x}) = \| y - H T^t \bar{x} \|^2 + \mu \| D T^t \bar{x} \|^2.

   这就等于在一个更小的维度（M 维而非 N 维）上做优化。

3. 闭式解的形式：

   \hat{\bar{x}} = \arg \min_{\bar{x} \in \mathbb{R}^M} \bar{J}_{\text{PLS}}(\bar{x})

   对这个较小问题求导、求极值后，最终我们可以得到 \bar{x} 的闭式解，即：

   \hat{\bar{x}} = \left[ T H^t H T^t + \mu T D^t D T^t \right]^{-1} T H^t y

4. 回填原始变量：

   \hat{x} = T^t \hat{\bar{x}}

   拿到 \hat{\bar{x}} 后，再通过 x = T^t \hat{\bar{x}} 把它扩展回去就得到最终的 \hat{x}。得到最终解：

   \hat{x} = T^t \left[ T \left( H^t H + \mu D^t D \right) T^t \right]^{-1} T H^t y

选择矩阵 T的本质是选出未被约束的维度来处理，可以将约束问题简化为更小的未约束问题，而且在很多实际应用中，受严格等式约束的那些像素（或数据点）往往只占据整个空间中的很小一部分，把它们从求解问题中剔除出去，可以大大降低计算量。

举个通俗易懂的例子，就像在一个杂乱的「箱子」里选出要处理的零件（定义选择矩阵 T），在干净整齐的「桌子」上（低维子空间）修理好零件（求解无约束优化问题），再把修好的零件放回原来的「箱子」适当位置（用 T^t 回填），这可比在箱子里直接修零件简单多了。

利用拉格朗日方法求解带等式约束的优化问题

上述的理论没有任何问题，问题在于实践中求不了这一大串的逆矩阵，在 Wiener Hunt 的正则化问题中，我们采用了循环近似矩阵和对角矩阵的傅立叶基分解，但是不适用于这里的带约束问题。所以我们就需要一种新方案来逼近最优解——也就是用拉格朗日方法（包括增广拉格朗日和迭代算法）来处理。

问题描述

原始（无约束）目标函数：

首先，我们的目标函数仍然是：

等式约束：

除此之外，我们还有一类约束：某些像素 x_p 必须为零。

用矩阵形式表示为：

拉格朗日形式：

为将约束纳入目标函数，我们引入一个拉格朗日乘子 \ell，把约束条件 嵌入 到目标里：

当优化这个拉格朗日函数时，\ell 会自动调整，确保最终解满足 \tilde{T}x=0

闭式解

理论上，通过对 \mathcal{L}(x,\ell) 同时对 x 和 \ell 求偏导并设为零，我们可以得到一个闭式解的表达式。

定义:

解为：

这个理论非常漂亮，但是实际求解不好用，因为要算 Q^{-1} 本身就够麻烦，后面还嵌套了 (\tilde{T} Q^{-1} \tilde{T}^t)^{-1}，更是没法算

迭代算法求解：从闭式转向迭代

由于闭式解在高维场景下难以直接实现，我们转而求助迭代方法，这也是常用方法。通过迭代，我们一步步逼近最优解，不要求一次到位。

1. 初始化拉格朗日乘子 \ell^{[0]} 为某个初值。

2. 在每次迭代中：

   • 更新 x^{[k+1]}：
   x^{[k+1]} = \arg \min_x \mathcal{L}(x, \ell^{[k]}) = (H^tH + \mu D^tD)^{-1}(H^t y - \tfrac{1}{2}\tilde{T}^t \ell^{[k]})

   这一步本质上是求解无约束问题的优化，多了一个项 \tilde{T}^t \ell^{[k]} 来让 x 满足约束。

   \ell^{[k+1]} = \ell^{[k]} + 2 \rho \tilde{T} x^{[k+1]}
   • 更新拉格朗日乘子 \ell^{[k+1]}：
     \ell^{[k+1]} = \ell^{[k]} + 2 \rho \tilde{T} x^{[k+1]}

   在这一步中，如果 x^{[k+1]} 没满足 \tilde{T}x=0，\ell^{[k+1]} 就会朝着强制 x 满足约束的方向改变，下次迭代会进一步推进 x 向可行解靠拢。

我们从 闭式解 转到了 迭代解，仅仅是因为闭式解没法理论计算而已

迭代算法扩展：增广拉格朗日方法

上述方法的理论很好，但是在实际实现中，收敛速度可能不理想或者过程不够稳健。为此，我们可以再给他优化一下，在拉格朗日框架中加入一个额外的二次罚项，引入松弛变量，这就是增广拉格朗日(Augmented Lagrangian)方法。

增广拉格朗日形式

原始目标函数与约束保持不变：

增广拉格朗日函数：

这里的 \rho \|\tilde{T}x\|^2 项在 x 偏离可行域时会产生额外的惩罚，让迭代向满足约束的方向更快地收敛。

迭代更新

• 更新 x^{[k+1]}：

  x^{[k+1]} = (H^tH + \mu D^tD + \rho \tilde{T}^t \tilde{T})^{-1} (H^t y - \tfrac{1}{2}\tilde{T}^t \ell^{[k]})

  相比之前多了 \rho \tilde{T}^t \tilde{T} 这项，进一步强调约束。

• 更新 \ell^{[k+1]}：

  \ell^{[k+1]} = \ell^{[k]} + 2 \rho \tilde{T} x^{[k+1]}

  同理，如果 x^{[k+1]} 未严格满足约束，\ell 会更剧烈地纠正下一次迭代的解。

引入松弛变量

我们引入松弛变量 s 来进一步简化约束处理方式。通过让 x 和 s 分离，我们的约束变成了 x=s，然后再强制 s=0 在特定位置。增广拉格朗日方法中加入 s，让问题的分解更为直接。

用 s 来表示约束位置的值（最终要求 s_p=0 在受约束的像素上）：

迭代更新算法：

• 更新 x^{[k+1]}：

• 更新 s_p^{[k+1]}：

  s_p^{[k+1]} = \begin{cases} x_p^{[k+1]} + \ell_p^{[k]}/(2\rho) & \text{for } p \in \mathcal{S} \\ 0 & \text{for } p \in \bar{\mathcal{S}} \end{cases}

  对受约束像素设为0，对非受约束像素更新为 x_p^{[k+1]} + \ell_p^{[k]}/(2\rho)

• 更新 \ell^{[k+1]}：

  \ell^{[k+1]} = \ell^{[k]} + 2 \rho \left( x^{[k+1]} - s^{[k+1]} \right)

通过这种方式，我们把 约束等于零 的问题转化为 两个变量必须相等 的问题，通过增广拉格朗日方法使得 x 和 s 相同，同时还能方便地直接在 s 上施加简单的条件（如 s=0）

迭代算法：利用 ADMM 迭代求解

ADMM (Alternating Direction Method of Multipliers) 方法通过分块优化，将问题拆分为针对 x 和 s 的两个子问题交替求解，再更新拉格朗日乘子 \ell。主要迭代步骤如下：

更新 x： 固定 s 和 \ell，求解最小化 \mathcal{L}_\rho(x, s^{[k]}, \ell^{[k]}) 对 x 的优化问题。

得到更新公式为：

当 H 和 D 具有卷积结构时，可以借助快速傅里叶变换(FFT)对角化的属性来快速求解这个方程。

更新 s： 有了新的 x^{[k+1]} 后，将其代入优化问题最小化 \mathcal{L}_\rho(x^{[k+1]}, s, \ell^{[k]}) 对 s。

对于 s 而言，我们需要满足：

• 当 p \in \bar{\mathcal{S}} 时，s_p = 0（等式约束）。
• 当 p \in \mathcal{S} 时，s_p \geq 0（不等式约束），并且希望 s_p 尽量接近 x_p^{[k+1]} + \ell_p^{[k]}/(2\rho)。

综合考虑后：

这里的 \max\{0,\cdot\} 操作保证了非负性约束；对于被约束为零的像素，直接设为0即可。

更新拉格朗日乘子 \ell：

最后，更新乘子以反映 x 和 s 一致性程度：

若 x_p^{[k+1]} 和 s_p^{[k+1]} 不一致，\ell_p^{[k+1]} 会在下一次迭代中推动 x 或 s 更好地满足约束条件。

方法特点

• 多重约束：该方法能够同时处理等式和不等式约束，通过对 s 的简单操作（如截断到非负值或直接设为0）实现复杂约束条件。
• 可扩展性：ADMM 适合大规模问题，因为每一步的子问题求解相对简单，且可以利用FFT加速 x 更新步骤的运算。

总结

问题背景

在许多图像与信号处理任务中，我们面临所谓的病态逆问题。以去卷积为例：我们观测到的是卷积后且带有噪声的信号 y = Hx + e，其中 H 为卷积算子、x 为原始信号。要从 y 中恢复 x，这属于逆问题。然而，在实际中，H 通常不可逆或近似不可逆，并且观测信号中包含噪声，导致问题严重病态：微小的噪声会在反演过程中被无限放大。

为避免这种不稳定性，我们引入正则化，例如加入 \mu \| Dx \|^2 项（D 为梯度算子）来约束解的平滑性，使问题有所谓的“稳健”解。然而，仅有正则化还不够，如果想让解 x 更具物理意义或先验知识（比如信号非负、仅在特定区域有值等），还必须加入约束条件。其中主要包括：

• 正值性约束：x_p \geq 0，确保解不出现物理上不合理的负值。
• 支持性约束：x_p=0 对于特定区域外的像素，确保信号只存在于已知的支持范围内。

综合这些，问题变为带约束的二次优化问题。

基本优化模型的建立

在加入正则化后，我们的优化目标函数为：

不考虑约束条件时，J_{\text{PLS}} 的无约束最优解有闭式形式：

若能求出 (H^T H + \mu D^T D)^{-1}，就能直接得到 \hat{x}_{\text{PLS}}​。在简单、低维情形下这可行。为提高效率，对于卷积算子 H 和梯度算子 D，我们可利用 FFT 将其对角化，将求逆复杂度降至对角元素求逆，从而快速获得无约束情况下的解。这为去卷积问题在无约束、仅有简单正则化的条件下提供了高效解法。

但是加入严格的等式约束（支持性）和不等式约束（非负性）时，仅靠 FFT 对角化无法直接求得被约束优化问题的闭式解。这是一个严格凸的二次优化问题，具有唯一解。我们需要在解中同时满足两类约束，从而将问题定义为：

该问题同时具备凸目标和凸约束，理论上有唯一最优解。但是，直接求解依旧困难。

增广拉格朗日与 ADMM

面对不可直接求逆的复杂约束问题，我们求助于增广拉格朗日方法和ADMM（交替方向乘子法）。

在有约束的优化问题中，这两个工具让我们不用一次性求出难以处理的逆矩阵，而是通过迭代求解，逐步接近最优解。每一轮迭代中，我们只需求解相对简单的子问题，或者做简单的投影与截断（如确保 x_p\ge0 简单地截断负值为0）

截断操作与选择矩阵：减少问题维度

对于等式约束 x_p=0，我们引入选择矩阵 T 将已知为零的维度剔除，从而在缩减维度的子空间中求解无约束子问题。这样大大减少了计算量和复杂度。在求解完成后，再通过 T^t 将解映射回原始空间，自动满足这些等式约束。

这一过程相当于剔除了不必要的元素，只在较小维度上计算，然后将结果填回原问题维度中。

从闭式解转向迭代求解

虽然理论上可以写出增广拉格朗日形式的闭式解，但在高维情况下不现实。于是，我们采用迭代算法：

• 每一步解决形如 (H^T H + \mu D^T D + \rho \tilde{T}^t\tilde{T})^{-1} 这类线性系统，可以用FFT加速子问题求解。
• 更新拉格朗日乘子，推动解满足约束。
• 如果有不等式约束（如非负性），在更新松弛变量（或辅助变量 s）时简单投影到可行域（如 \max\{0, \cdot\}）即可。

增广拉格朗日 + 松弛变量 + ADMM 实现全面解决方案

引入松弛变量 s 将约束（比如 x_p=0 或 x_p\ge0）转化为对 s 的简单限制（例如 s_p=0 或 s_p \ge0），同时让 x 与 s 保持一致。这样分而治之的策略，配合增广拉格朗日的二次惩罚项，提升了算法的收敛速度和稳定性。

最后以 ADMM 迭代为主干架构，我们获得了一个可扩展、可并行、适用于大规模问题的求解框架：

1. 更新 x：利用FFT等技术快速求解子问题。
2. 更新 s：对每个像素执行简单的非负投影或赋零操作即可轻松满足约束。
3. 更新 \ell（拉格朗日乘子）：不断约束 x 与 s 的一致性，从而逐步逼近满足所有约束条件的全局最优解。

这一整套方法不仅限于简单的正值性和支持性约束，还可扩展至更多类型的约束与更复杂的正则化项。利用 ADMM 的灵活分块思想，我们能处理更广泛的信号和图像反演问题，如三维重建、多光谱图像处理、时间序列分析等。