一篇论文复现的整体思路和复现记录(二，算法篇)

竹外桃花依旧

1504人浏览 · 2024-05-08 00:50:49

竹外桃花依旧 · 2024-05-08 00:50:49 发布

典型问题形式

考虑以下凸问题：
$\begin{align} \min\limits_{x_1,x_2} \quad &f_1(x_1)+f_2(x_2) \\ s.t.\quad &A_1x_1+A_2x_2 = b \end{align}\tag{1}$
要求 $f_1$ , $f_2$ 是适当的闭凸函数，但不要求是光滑的。 $x_1\in\mathbb{R}^n$ , $x_2\in\mathbb{R}^m$ , $A_1\in\mathbb{R}^{p\times n}$ , $A_2\in\mathbb{R}^{p\times m}$ , $b\in\mathbb{R}^{p}$
问题的特点在于目标函数可以分为彼此分离的两块。但是变量都被线性约束结合在一起。常见的一些无约束和带约束的问题可以表现为这一形式。

1

可以分为两块的无约束优化问题：
$\min\limits_x\quad f_1(x)+f_2(x)$
引入新的变量 $z$ ，令 $x = z$ , 将问题转换为：
$\begin{aligned} \min\limits_{x, z} \quad& f_1(x)+f_2(z)\\ s.t.\quad&x-z = 0 \end{aligned}$

2

带线性变换的无约束优化问题：
$\min\limits_x\quad f_1(x)+f_2(Ax)$
引入新的变量 $z$ ，令 $z = A x$ , 将问题转换为：
$\begin{aligned} \min\limits_{x, z} \quad& f_1(x)+f_2(z)\\ s.t.\quad&Ax-z = 0 \end{aligned}$

3

凸集 $C\subset\mathbb{R}^n$ 上的约束优化问题
$\begin{aligned} \min\limits_{x} \quad& f(x)\\ s.t.\quad&Ax\in C \end{aligned}$
$I_C(z)$ 是集合 $C$ 的示性函数，引入约束 $z = A x$ ，则问题转化为：
$\begin{aligned} \min\limits_{x,z} \quad& f(x)+I_C(z)\\ s.t.\quad&Ax-z=0 \end{aligned}$
（Here is a problem: 示性函数是什么东西）

4

全局一致性问题
$\min\limits_x\quad\sum_{i=1}^N\phi_i(x)$
令 $x = z$ ，并将x复制 $N$ 份，分别为 $x_i$ , 那么问题转换为：
$\begin{aligned} \min\limits_{x,z} \quad &\sum_{i=1}^N\phi_i(x_i)\\ s.t.\quad&x_i-z=0 \end{aligned}$

这个例4和文章中所显示的形式上比较相似。

增广拉格朗日函数法和交替方向乘子法（ADMM）

增广拉格朗日函数法

写出前述问题（1）的增广拉格朗日函数
$\begin{align} L_{\rho}(x_1,x_2,y) = f_1(x_1)+f_2(x_2)+y^T(A_1x_1+A_2x_2-b)+\frac{\rho}{2}{\Vert{A_1x_1+A_2x_2-b}\Vert}^2_2 \end{align}\tag{2}$

其中， $\rho>0$ 是二次罚项的系数。

常见求解带约束问题的增广拉格朗日函数法为如下更新：
$\begin{align} (x_1^{k+1},x_2^{k+1}) &= \mathop{argmin}\limits_{x_1,x_2}L_{\rho}(x_1,x_2,y^k),\tag{3}\\ y^{k+1} &= y^k + \tau\rho(A_1x_1^{k+1}+A_2x_2^{k+1}-b)\tag{4} \end{align}$

交替方向乘子法： Alternating direction method of multipliers, ADMM

交替方向乘子法的基本思路：

第一步迭代（3）的同时对 $x_1$ 和 $x_2$ 进行优化有时候比较困难，而固定一个变量求解关于另一个变量的极小化问题有时候可能比较简单，因此可以考虑：对 $x_1$ 和 $x_2$ 交替求极小值：

迭代格式可以总结如下：

$\begin{align} x_1^{k+1} &= \mathop{argmin}\limits_{x_1}L_{\rho}(x_1,x_2^k, y^k)\tag{5}\\ x_2^{k+1} &= \mathop{argmin}\limits_{x_1}L_{\rho}(x_1^{k+1},x_2, y^k)\tag{6}\\ y^{k+1} &= y^k + \tau\rho(A_1x_1^{k+1}+A_2x_2^{k+1}-b)\tag{7} \end{align}$

其中， $\tau$ 为步长，其取值范围通常为 $(0,\frac{1+\sqrt{5}}{2}]$

原问题最优性条件

因为 $f_1,f_2$ 均为闭凸函数，约束为线性函数，所以当Slater条件成立的时候，可以使用凸优化问题的KKT条件作为交替方向乘子法的收敛准则。问题（1）的拉格朗日函数为
$\begin{aligned} L(x_1,x_2,y) = f_1(x_1)+f_2(x_2)+y^T(A_1x_1+A_2x_2-b) \end{aligned}$
根据最优性条件定理，若 $x_1^*,x_2^*$ 为问题（1）的最优解， $y^*$ 为对应的拉格朗日乘子，则满足以下条件：
$\begin{align} 0\in\partial_{x_1}L(x_1^*, x_2^*,y^*) &= \partial f_1(x^*_1)+A_1^T y^*,\tag{8a}\\ 0\in\partial_{x_2}L(x_1^*, x_2^*,y^*) &= \partial f_2(x^*_2)+A_2^T y^*,\tag{8b}\\ A_1x^*_1+A_2x^*_2 &= b\tag{8c} \end{align}$
其中，前两者成为原始可行性条件，后两者称为对偶可行性条件。

ADMM单步迭代最优性条件

$x_2$ 的更新步骤：
$\begin{align*} x_2^k = \mathop{argmin}\limits_x{f_2(x)+\frac{\rho}{2}{\Vert{A_1x_1^k+A_2x-b+\frac{y^{k-1}}{\rho}}\Vert}^2} \end{align*}$
根据最优性条件不难推出，
$\begin{align*} 0\in\partial f_2(x_2^k)+A^T_2[y^{k-1}+\rho(A_1x_1^k+A_2x^k_2-b)] \end{align*}\tag{9}$
当 $\tau=1$ 时候，根据（7）可知上式方括号中的表达式就是 $y^k$ ，最终有：
$\begin{align*} 0\in\partial f_2(x_2^k)+A^T_2y^k \end{align*}$
由 $x_1$ 的更新公式
$\begin{align*} x_1^k=\mathop{argmin}\limits_x\{f_1(x)+\frac{\rho}{2}{\Vert{A_1x+A_2x_2^{k-1}-b+\frac{y^{k-1}}{\rho}}\Vert}^2\} \end{align*}$
假设子问题能够精确求解，根据最优性条件，
$\begin{align*} 0\in\partial f_1(x_1^k)+A_1^T[\rho(A_1x_1^k+A_2x_2^{k-1}-b)+y^{k-1}] \end{align*}$
根据ADMM的第三式（7）取 $\tau=1$ 有
$\begin{align*} 0\in\partial f_1(x_1^k)+A_1^T(y^k+\rho A_2(x_2^{k-1}-x_2^k)） \end{align*}\tag{10}$
对比（8a)可知多出来的项为 $A_1^TA_2(x_2^{k-1}-x_2^k)$ 。因此要检测对偶可行性只需要检测残差
$s^k = A^T_1A_2(x_2^{k-1}-x_2^k)$
综上，当 $x_2$ 更新取到精确解且 $\tau=1$ 时候，判断ADMM收敛只需要检测前述两个残差 $r^k, s^k$ 是否充分小：
$\begin{align*} 0\approx||r^k||=\Vert A_1x_1^k+A_2x_2^k-b\Vert\quad 原始可行性\\ 0\approx||s^k||=\Vert A^T_1A_2(x_2^{k-1}-x_2^k)\Vert\quad 对偶可行性\\ \end{align*}\tag{11}$

常见问题和技巧

线性化

线性化技巧使用近似点项对子问题目标函数进行二次近似。
不失一般性，我们考虑第一个子问题，即：
$\begin{align*} \min_{x_1}f_1(x_1)+\frac{\rho}{2}{\Vert{A_1x_1-v^k}\Vert}^2 \end{align*}\tag{12}$
其中： $v^k = b-A_2x_2^k-\frac{1}{\rho}y^k$
当子问题目标函数可微的时候，线性化将问题（12）变为：
$\begin{align*} x^{k+1}_1 = \mathop{argmin}\limits_{x_1}{(\nabla f_1(x_1^k)+\rho A_1^T(A_1x_1^k-v^k))^Tx_1+\frac{1}{2\eta_k}{\Vert x_1-x_2^2\Vert}^2_2} \end{align*}$
其中， $\eta_k$ 是步长参数，这等价于做一步梯度下降。
当目标函数不可微的时候，可以考虑只将二次项线性化，即：
$\begin{align*} x^{k+1}_1 = \mathop{argmin}\limits_{x_1}{(f_1(x_1^k)+\rho A_1^T(A_1x_1^k-v^k))^Tx_1+\frac{1}{2\eta_k}{\Vert x_1-x_2^2\Vert}^2_2} \end{align*}$
这等价于做一步近似点梯度步。

缓存分解

如果目标函数中含有二次函数，例如 $f_1(x_1) = \frac{1}{2}{\Vert Cx_1-d\Vert}^2_2$ ，那么针对 $x_1$ 的更新（5）等价于求解线性方程组
$(C^TC+\rho A^T_1A_1)x_1=C^Td+\rho A_1^Tv^k$
虽然子问题有显式解，但是每步求解的复杂度仍然比较高，这时候可以考虑用缓存分解的方法，首先对 $C^TC+\rho A^T_1A_1$ 进行Cholesky分解，并缓存分解的结果，在每步迭代中，只需要求解简单的三角形方程组
当 $\rho$ 发生更新的时候，就要重新进行分解，特别地，当 $C^TC+\rho A^T_1A_1$ 一部分容易求逆，另一部分是低秩的情形时，可以使用矩阵辅助求逆公式（SMW）进行求逆。

优化转移

有时候为了方便求解子问题，可以用一个性质好的矩阵 $D$ 近似二次项 $A_1^TA_1$ ，此时子问题（12）替换为
$x_1^{k+1} = \mathop{argmin}\limits_{x_1}{f_1(x_1)+\frac{\rho}{2}{\Vert{A_1x_1-v^k}\Vert}^2 +\frac{\rho}{2}(x_1-x^k)^T(D-A_1^TA_1)(x_1-x^k)}$
这种方法也被称为优化转移。
通过选取合适的 $D$ ，当计算 $\mathop{argmin}\limits_{x_1}\{f_1(x_1)+\frac{\rho}{2}x_1^TDx_1\}$ 明显比计算 $\mathop{argmin}\limits_{x_1}\{f_1(x_1)+\frac{\rho}{2}x_1^TA_1TA_1x_1\}$ 要容易的时候，优化转移可以极大地简化子问题的计算。特别地，当 $D=\frac{\eta_k}{\rho}I$ 时，优化转移等价于做单步的近似点梯度步。

二次罚项系数的动态调节

原始可行性和对偶可行性分别用 $r^k\Vert$ 和 $s^k\Vert$ 度量。

求解过程中二次罚项系数 $\rho$ 太大会导致原始可行性 $r^k\Vert$ 下降很快，但是对偶可行性 $s^k\Vert$ 下降很慢；二次罚项系数太小，则会有相反的效果，这样都会导致收敛比较慢或者得到的解的可行性很差。
一个自然的想法是在每次迭代的时候动态调节惩罚系数 $\rho$ 的大小，从而使得原始可行性和对偶可行行能够以比较一致的速度下降到0.一个简单有效的方式是令
$\rho^{k+1} = \begin{cases} \gamma_p\rho^k, & \Vert r^k\Vert>\mu \Vert s^k\Vert \\ \frac{\rho^k}{\gamma_p}, & \Vert s^k\Vert>\mu \Vert r^k\Vert \\ \rho^k, &其他 \end{cases}$
其中， $\mu>1$ , $\gamma_p>1$ , $\gamma_d>1$ 是参数，常见的选择
为 $\mu=10, \gamma_p=\gamma_d=2$ .在迭代过长中将原始可行性 $r^k||$ 和对偶可行性 $s^k||$ 保持在彼此的 $\mu$ 倍内，如果发现 $r^k||$ 或 $s^k||$ 下降过慢就应该相应增大或者减小二次罚项系数 $\rho^k$ .

多块问题的ADMM

考虑有多块变量的情形
$\begin{align} \min_{x_1,x_2,\dots,x_N}\quad &f_1(x_1)+f_2(x_2)+\dots+f_N(x_N),\\ s.t. \quad &A_1x_1+A_2x_2+\dots+A_Nx_N =b. \end{align}\tag{13}$
这里， $f_i(x_i)$ 是闭函数， $x_i\in \mathbb{R}^{n_i}$ , $A_i\in \mathbb{R}^{m\times n_i}$
同样写出增广拉格朗日函数 $L_{\rho}(x_1, x_2, \dots, x_N, y)$ ，相应的多块ADMM迭代格式为
$\begin{align*} x_1^{k+1} &= \mathop{argmin}\limits_{x_1}L_{\rho}(x, x^{k}_2, \dots, x^{k}_N, y^k)\\ x_2^{k+1} &= \mathop{argmin}\limits_{x_1}L_{\rho}(x^{k+1}_1, x, \dots, x^k_N, y^k)\\ \dots \\ x_N^{k+1} &= \mathop{argmin}\limits_{x_1}L_{\rho}(x^{k+1}_1, x^{k+1}_2, \dots, x_N, y^k)\\ y^{k+1} &= y^k + \tau \rho (A_1x_1^{k+1}+A_2x_2^{k+1}+\dots+A_Nx_N^{k+1}-b) \end{align*}$
其中， $\tau\in(0, \frac{\sqrt{5}+1}{2})$ 为步长参数。

应用举例

限于篇幅和内容需要，部分应用举例已折叠。

Lasso 问题的Primal形式

Lasso 问题的对偶形式

广义Lasso问题

半定规划问题

ADMM求解半定规划问题

稀疏逆协方差矩阵估计

矩阵分离问题

图像去噪模型

全变差去模糊方法

ADMM 求解图像去模糊问题及实例

全局一致性优化问题

原问题为：
$\begin{align*} \min_{x_i,z}\quad &\sum_{i=1}^N\phi_i(x_i),\\ s.t. \quad &x_i-z=0, i = 1,2,\dots,N. \end{align*}$
增广拉格朗日函数：
$L(x_1,\dots, x_N, z, y_1,\dots, y_N) = \sum_{i=1}^N\phi_i(x_i) +\sum_{i=1}^N y_i^T(x_i-z)+\frac{\rho}{2}\sum_{i=1}^N{\Vert x_i-z\Vert}^2$
固定 $z^k$ , $y_i^k$ ，更新 $x_i$ 的公式为

$\begin{align} x_i^{k+1} = \mathop{argmin}\limits_{x}\{{\phi_i(x)+\frac{\rho}{2}{\Vert{ x_i-z+\frac{y_i^k}{\rho}}\Vert}^2}\} \end{align}\tag{21}$

注意：虽然表面上看，增广拉格朗日函数有 $(N + 1)$ 个变量块，但本质上还是两个变量块，这是因为在更新某个 $x_i$ 时候并未利用其他 $x_i$ ，所有 $x_i$ 可以看成一个整体。相应地，所有的乘子 $y_i$ 也可以看成一个整体。
迭代式（21）的具体计算依赖于 $\psi_i$ 的形式，在一般情况下，更新 $x_i$ 的表达式为
$x_i^{k+1}=prox_{\phi_i/\rho}(z^k-y_i^k/\rho)$
固定 $x_{i}^{k+1}$ , $y_i^k$ ，问题关于 $z$ 时二次函数，因此可以直接写出显式解：
$z^{k+1} = \frac{1}{N}\sum_{i=1}^N(x_i^{k+1}+y_i^k/\rho).$
综上，该问题的交替方向乘子法迭代格式为
$\begin{align*} x_i^{k+1}&=prox_{\phi_i/\rho}(z^k-y_i^k/\rho), , i = 1,2,\dots,N.\\ z^{k+1} &= \frac{1}{N}\sum_{i=1}^N(x_i^{k+1}+y_i^k/\rho).\\ y_i^{k+1} &= y_i^k+\tau\rho(x_i^{k+1}-z^{k+1}), i = 1,2,\dots,N. \end{align*}$

与上一篇的关系

因此，在上一篇文章中所提到的
$\begin{align*} \text{Minimized} \quad & \{-\sum_{j=0}^n L_j(x_j)\}\\ \text{subject to} \quad & y-z_j= 0\quad j = 0,\dots,n \end{align*}\tag{23}$

应用ADMM解决这一问题：
$\begin{align*} y^{k+1}&:=\mathop{argmin}\limits_y\{\langle\sum_{j=0}^np_j^k, y\rangle+\frac{\rho}{2}\sum_{j=0}^n{\Vert y-z_j^k\Vert}^2 \}\\ z^{k+1}_j&:=\mathop{argmin}\limits_{z_j}\{-L_j(z_j)-\langle p_j^k, z_j\rangle+\frac{\rho}{2}\sum_{j=0}^n{\Vert y^{k+1}-z_j\Vert}^2 \}\\ p^{k+1}_j&:=p^{k}_j+\rho(y^{k+1}-z^{k+1}_j) \end{align*}$

迭代过程为：
$y^{k+1} = \frac{1}{n+1}\sum_{n+1}^nz_j^k+-\frac{1}{(n+1)\rho}\sum_{j=0}^np_j^k$
$\begin{align*} \{x_i^*,i\in\mathcal{A}_j \} :=\mathop{argmin}\limits_{x_i \in \chi_i, i\in\mathcal{A}_j}\left\{\sum_{i\in\mathcal{A}_j}g_i(x_i)+\frac{\rho}{2}{\left\Vert\left[ \begin{matrix} max\{0, \lambda_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in A_j \psi_i(x_i)+\lambda_{p_j^k}})\}\\ \nu_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in A_j}A_i(x_i)+\lambda_{p_j^k}) \end{matrix} \right]\right\Vert}^2\right\} \end{align*}$
$z_j = \left[ \begin{matrix} \lambda_{zj}\\ \nu_{zj} \end{matrix} \right] = \left[ \begin{matrix} max\{0, \lambda_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in \mathcal{A}_j \psi_i(x_i)+\lambda_{p_j^k}})\}\\ \nu_{y^{k+1}}+\frac{1}{\rho}(\sum_{i\in \mathcal{A}_j}\mathcal{A}_i(x_i)+\lambda_{p_j^k}) \end{matrix} \right]$
$\begin{aligned} p^{k+1}_j:=p^{k}_j+\rho(y^{k+1}-z^{k+1}_j) \end{aligned}$