最优化学习笔记： 交替方向乘子法（2）

8.6 交替方向乘子法（续1）

8.6.4 应用举例

本小节给出一些交替方向乘子法的应用实例. 在实际中，大多数问题并不直接具有问题 (8.6.1) 的形式. 我们需要通过一系列拆分技巧将问题化成 ADMM 的标准形式，同时要求每一个子问题尽量容易求解. 注意，对同一个问题可能有多种拆分方式，不同方式导出的最终算法可能差异巨大，应当选择最容易求解的拆分方式.

1. LASSO 问题
LASSO 问题为
$$\min \mu \Vert x{\Vert }_1+\frac{1}{2}\Vert Ax-b{\Vert }^2$$

这是典型的无约束复合优化问题，将其写成 ADMM 标准问题形式：
$$\begin{array}{rl}\underset{x,z}{\min }& f(x)+h(z)\\ s.t.& x=z\end{array}$$

其中 $f(x)=\frac{1}{2}\Vert Ax-b{\Vert }^2$, $h(z)=\mu \Vert z{\Vert }_1$.
对于此问题，交替方向乘子法迭代格式为
$$\begin{array}{rl}& x^{k+1}=\underset{x}{\mathrm{argmin}}\left\{\frac{1}{2}\Vert Ax-b{\Vert }^2+\frac{\rho }{2}\Vert x-z^k+\frac{1}{\rho }{y}^k{\Vert }_2^2\right\}=(A^{\mathrm{T}}A+\rho I{)}^{-1}(A^{\mathrm{T}}b+\rho z^k-y^k)\\ & z^{k+1}=\underset{z}{\mathrm{argmin}}\left\{\mu \Vert z{\Vert }_1+\frac{\rho }{2}\Vert x^{k+1}-z+\frac{1}{\rho }{y}^k{\Vert }^2\right\}={\mathrm{prox}}_{(\mu /\rho )\Vert \cdot {\Vert }_1}\left(x^{k+1}+\frac{1}{\rho }{y}^k\right)\\ & y^{k+1}=y^k+\tau \rho (x^{k+1}-z^{k+1})\end{array}$$

注意，因为 $\rho >0$, 所以 $A^{\mathrm{T}}A+\rho I$ 总是可逆的. $x$ 迭代本质上是计算一个岭回归问题 (${\ell }_2$ 范数平方正则化的最小二乘问题); 而对 $z$ 的更新为 ${\ell }_1$ 范数的邻近算子，同样有显式解. 在求解 $x$ 迭代时，若使用固定的罚因子 $\rho$, 我们可以缓存矩阵 $A^{\mathrm{T}}A+\rho I$ 的初始分解，从而减小后续迭代中的计算量. 需要注意的是，在 LASSO 问题中，矩阵 $A\in {\mathbb{R}}^{m\times n}$ 通常有较多的列 (即 $m\ll n$), 因此 $A^{\mathrm{T}}A\in {\mathbb{R}}^{n\times n}$ 是一个低秩矩阵，二次罚项的作用就是将 $A^{\mathrm{T}}A$ 增加了一个正定项. 该 ADMM 主要运算量来自更新 $x$ 变量时求解线性方程组，复杂度为 $\mathcal{O}(n^3)$ (若使用缓存分解技术或 SMW 公式则可进一步降低每次迭代的运算量).

接下来考虑 LASSO 问题的对偶问题
$$\begin{array}{rl}\min & b^{\mathrm{T}}y+\frac{1}{2}\Vert y{\Vert }^2\\ s.t.& \Vert A^{\mathrm{T}}y{\Vert }_{\infty }⩽\mu \end{array}(\mathrm{8.6.26})$$

将 $\Vert A^{\mathrm{T}}y{\Vert }_{\infty }⩽\mu$ 变成示性函数放在目标函数中，并引入约束 $A^{\mathrm{T}}y+z=0$, 可以得到如下等价问题：
$$\begin{array}{rl}\min & \underset{f(y)}{\underbrace{b^{\mathrm{T}}y+\frac{1}{2}\Vert y{\Vert }^2}}+\underset{h(z)}{\underbrace{I_{\Vert z{\Vert }_{\infty }⩽\mu }(z)}}\\ s.t.& A^{\mathrm{T}}y+z=0\end{array}(\mathrm{8.6.27})$$

对约束 $A^{\mathrm{T}}y+z=0$ 引入乘子 $x$, 对偶问题的增广拉格朗日函数为
$$L_{\rho }(y,z,x)=b^{\mathrm{T}}y+\frac{1}{2}\Vert y{\Vert }^2+I_{\Vert z{\Vert }_{\infty }⩽\mu }(z)-x^{\mathrm{T}}(A^{\mathrm{T}}y+z)+\frac{\rho }{2}\Vert A^{\mathrm{T}}y+z{\Vert }^2$$

在这里故意引入符号 $x$ 作为拉格朗日乘子，$x$ 恰好对应的是原始问题的自变量.
以下说明如何求解每个子问题. 当固定 $y,x$ 时，对 $z$ 的更新即向无穷范数球 { z ∣ ∥ z ∥ ∞ ⩽ μ } \{z|\|z\|_\infty\leqslant\mu\} {z∣∥z∥∞​⩽μ} 做欧几里得投影，即将每个分量截断在区间 $[-\mu ,\mu ]$ 中；当固定 $z,x$ 时，对 $y$ 的更新即求解线性方程组
$$(I+\rho A{A}^{\mathrm{T}})y=A(x^k-\rho z^{k+1})-b$$

因此得到 ADMM 迭代格式为
$$\begin{array}{rl}& z^{k+1}={\mathcal{P}}_{\Vert z{\Vert }_{\infty }⩽\mu }\left(\frac{x^k}{\rho }-A^{\mathrm{T}}{y}^k\right)\\ & y^{k+1}=(I+\rho A{A}^{\mathrm{T}}{)}^{-1}(A(x^k-\rho z^{k+1})-b)\\ & x^{k+1}=x^k-\tau \rho (A^{\mathrm{T}}{y}^{k+1}+z^{k+1})\end{array}$$

虽然 ADMM 应用于对偶问题也需要求解一个线性方程组，但由于 LASSO 问题的特殊性 $(m\ll n)$, 求解 $y$ 更新的线性方程组需要的计算量是 $\mathcal{O}(m^3)$, 使用缓存分解技巧后可进一步降低至 $\mathcal{O}(m^2)$, 这大大小于针对原始问题的 ADMM.
对原始问题，另一种可能的拆分方法是
$$\begin{array}{rl}\min & \underset{f(x)}{\underbrace{\mu \Vert x{\Vert }_1}}+\underset{h(z)}{\underbrace{\frac{1}{2}\Vert z{\Vert }^2}}\\ s.t.& z=Ax-b\end{array}$$

可以写出其增广拉格朗日函数为
$$L_{\rho }(x,z,y)=\mu \Vert x{\Vert }_1+\frac{1}{2}\Vert z{\Vert }^2+y^{\mathrm{T}}(Ax-b-z)+\frac{\rho }{2}\Vert Ax-b-z{\Vert }^2$$

但是如果对这种拆分方式使用 ADMM, 求解 $x$ 的更新本质上还是在求解一个 LASSO 问题！这种变形方式将问题绕回了起点，因此并不是一个实用的方法.

用同第 6.2 节中一样的 $A$ 和 $b$, 并取 $\mu =10^{-}3$, 分别使用 ADMM 求解原始问题和对偶问题，这里取 $\tau =1.618$, 原始问题和对偶问题的参数 $\rho$ 分别为 0.01 和 100, 终止条件设为 $|f(x^k)-f(x^{k-1})|<10^{-8}$ 和最大迭代步数 2000. 此外，对于原始问题的 ADMM, 我们还添加终止条件 $\Vert x^k-z^k\Vert <10^{-10}$; 相应地，对于对偶问题，额外的终止条件取为 $\Vert A^{\mathrm{T}}{y}^k+z^k\Vert <10^{-10}$. 算法结果见图8.9. 这里的 ADMM 没有使用连续化策略来调整 $\mu$, 因此可以看出 ADMM 相对其他算法的强大之处．此外，对于这个例子，求解原问题需要的迭代步数较少，但求解对偶问题每一次迭代所需要的时间更短，综合来看 ADMM 求解对偶问题时更快．

import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['Microsoft YaHei']

def f_prox(x, threshold):
    return np.sign(x) * np.maximum(np.abs(x) - threshold, 0)
def solve_primal(A, b, mu, rho, tau, max_iter=2000, tol=1e-8):
    m, n = A.shape
    x = np.zeros((n, 1))
    z = np.zeros((n, 1))
    y = np.zeros((n, 1))
    k = 0
    history = []
    while k < max_iter:
        x_new = np.linalg.solve(A.T @ A + rho * np.eye(n), A.T @ b + rho * z - y)
        z_new = f_prox(x_new + y / rho, mu / rho)
        y_new = y + tau * rho * (x_new - z_new)
        obj_val = mu * np.linalg.norm(x, ord=1) + 0.5 * np.linalg.norm(A @ x - b)**2
        history.append(obj_val)
        if k > 0 and (np.abs(history[-1] - history[-2]) < tol or np.linalg.norm(x_new - z_new) < tol):
            break
        x = x_new
        z = z_new
        y = y_new
        k += 1
    return history, k
def solve_dual(A, b, mu, rho, tau, max_iter=2000, tol=1e-8):
    m, n = A.shape
    y = np.zeros((m, 1))
    z = np.zeros((n, 1))
    x = np.zeros((n, 1))
    k = 0
    history = []
    while k < max_iter:
        z_new = np.maximum(np.minimum(x / rho - A.T @ y, mu), -mu)
        y_new = np.linalg.solve(np.eye(m) + rho * A @ A.T, A @ (x - rho * z_new) - b)
        x_new = x - tau * rho * (A.T @ y_new + z_new)
        obj_val = mu * np.linalg.norm(x, ord=1) + 0.5 * np.linalg.norm(A @ x - b)**2
        history.append(obj_val)
        if k > 0 and (np.abs(history[-1] - history[-2]) < tol or np.linalg.norm(A.T @ y + z) < tol):
            break
        x = x_new
        z = z_new
        y = y_new
        k += 1
    return history, k

np.random.seed(42)
m = 512
n = 1024
A = np.random.randn(m, n)
u = sp.random(n, 1, density=0.1).toarray()
b = np.dot(A, u)
mu = 1e-3
tau = 1.618
rho_primal = 0.01
rho_dual = 100
primal_history, k_primal = solve_primal(A, b, mu, rho_primal, tau)
dual_history, k_dual = solve_dual(A, b, mu, rho_dual, tau)
plt.figure(figsize=(10, 6))
# 创建对数坐标轴的折线图
plt.semilogy((np.array(dual_history) - dual_history[-1]) / (dual_history[-1]), label='ADMM解对偶问题')
plt.semilogy((np.array(primal_history) - primal_history[-1]) / (primal_history[-1]), linestyle='--', label='ADMM解原问题')
plt.xlabel('迭代次数')
plt.ylabel('$(f(x^k) - f^*) / f^*$')
plt.legend()
plt.xlim(0)
plt.grid(False)
plt.show()

2. 广义 LASSO 问题
广义 LASSO 问题的定义为
$$\underset{x}{\min }\mu \Vert Fx{\Vert }_1+\frac{1}{2}\Vert Ax-b{\Vert }^2(\mathrm{8.6.28})$$

在 LASSO 问题中，增加 $\Vert x{\Vert }_1$ 项是要保证 $x$ 的稀疏性，而对许多问题，$x$ 本身不稀疏，但在某种变换下是稀疏的. 一个重要的例子是当 $F\in {\mathbb{R}}^{(n-1)\times n}$ 是一阶差分矩阵
$$F_{ij}=\{\begin{array}{ll}1,& j=i+1\\ -1,& j=i\\ 0,& \text{其他}\end{array}$$

且 $A=I$ 时，广义 LASSO 问题为
$$\underset{x}{\min }\frac{1}{2}\Vert x-b{\Vert }^2+\mu \sum _{i=1}^{n-1}|x_{i+1}-x_i|$$

这个问题就是图像去噪问题的 TV 模型; 当 $A=I$ 且 $F$ 是二阶差分矩阵时，问题 (8.6.28) 被称为一范数趋势滤波.
下面介绍如何使用 ADMM 求解问题 (8.6.28)．通过引入约束 $Fx=z$，将问题 (8.6.28) 写为交替方向乘子法所对应问题的形式：
$$\begin{array}{rl}\underset{x,z}{\min }& \frac{1}{2}\Vert Ax-b{\Vert }^2+\mu \Vert z{\Vert }_1\\ s.t.& Fx-z=0\end{array}(\mathrm{8.6.29})$$

引入乘子 $y$, 其增广拉格朗日函数为
$$L_{\rho }(x,z,y)=\frac{1}{2}\Vert Ax-b{\Vert }^2+\mu \Vert z{\Vert }_1+y^{\mathrm{T}}(Fx-z)+\frac{\rho }{2}\Vert Fx-z{\Vert }^2$$

此问题的 $x$ 迭代是求解方程组
$$(A^{\mathrm{T}}A+\rho F^{\mathrm{T}}F)x=A^{\mathrm{T}}b+\rho F^{\mathrm{T}}\left(z^k-\frac{y^k}{\rho }\right)$$

而 $z$ 迭代依然通过 ${\ell }_1$ 范数的邻近算子.
因此交替方向乘子法所产生的迭代为
$$\begin{array}{rl}& x^{k+1}=(A^{\mathrm{T}}A+\rho F^{\mathrm{T}}F{)}^{-1}\left(A^{\mathrm{T}}b+\rho F^{\mathrm{T}}\left(z^k-\frac{y^k}{\rho }\right)\right)\\ & z^{k+1}={\mathrm{prox}}_{(\mu /\rho )\Vert \cdot {\Vert }_1}\left(F{x}^{k+1}+\frac{y^k}{\rho }\right)\\ & y^{k+1}=y^k+\tau \rho (F{x}^{k+1}-z^{k+1})\end{array}$$

在这个问题中，ADMM 迭代的主要计算量仍在 $x$ 更新上. 由于图像去噪问题涉及的变量数量可能有上百万，这一步迭代就要涉及求解一个百万量级的线性方程组. 在这种数据规模下不适合再使用矩阵分解的方式求解方程组，而应当注意系数矩阵的特殊结构. 对于全变差去噪问题，$A^{\mathrm{T}}A+\rho F^{\mathrm{T}}F$ 是三对角矩阵，所以此时 $x$ 迭代可以在 $\mathcal{O}(n)$ 的时间复杂度内解决; 对于图像去模糊问题，$A$ 是卷积算子，则利用傅里叶变换可将求解方程组的复杂度降低至 $\mathcal{O}(n\log n)$; 对于一范数趋势滤波问题，$A^{\mathrm{T}}A+\rho F^{\mathrm{T}}F$ 是五对角矩阵，所以 $x$ 迭代仍可以在 $\mathcal{O}(n)$ 的时间复杂度内解决.

五对角矩阵是一种特殊的方阵，其除了主对角线外，还有上方和下方各有两条对角线。换句话说，五对角矩阵中除了主对角线外，只有相邻的第一条对角线和第二条对角线上有非零元素。
$$\mathbf{A}=\left(\begin{array}{cccccc}{a}_1& b_1& c_1& & & \\ d_2& a_2& b_2& c_2& & \\ e_3& d_3& a_3& \ddots & \ddots & \\ & e_4& \ddots & \ddots & & c_{n-2}\\ & & \ddots & & a_{n-1}& b_{n-1}\\ & & & e_n& d_n& a_n\end{array}\right)$$

3. 逆协方差矩阵估计
第三章介绍了概率图模型和逆协方差矩阵估计问题．该问题的基本形式是
$$\underset{X}{\min }⟨S,X⟩-\ln \det X+\mu \Vert X{\Vert }_1(\mathrm{8.6.30})$$

其中 $S$ 是已知的对称矩阵，通常由样本协方差矩阵得到. 变量 $X\in {\mathcal{S}}_{++}^n,\Vert \cdot {\Vert }_1$ 定义为矩阵所有元素绝对值的和. 接下来说明如何应用交替方向乘子法.
目标函数由光滑项和非光滑项组成，因此引入约束 $X=Z$ 将问题的两部分分离：
$$\begin{array}{rl}\min & \underset{f(X)}{\underbrace{⟨S,X⟩-\ln \det X}}+\underset{h(Z)}{\underbrace{\mu \Vert Z{\Vert }_1}}\\ s.t.& X=Z\end{array}$$

引入乘子 $U$ 作用在约束 $X-Z=0$ 上，可得增广拉格朗日函数为
$$L_{\rho }(X,Z,U)=⟨S,X⟩-\ln \det X+\mu \Vert Z{\Vert }_1+⟨U,X-Z⟩+\frac{\rho }{2}\Vert X-Z{\Vert }_F^2$$

注意，针对矩阵情形应当使用 $F$ 范数替换 ${\ell }_2$ 范数作为增广拉格朗日函数法的罚项. 接下来分别写出 ADMM 子问题的显式解. 首先，固定 $Z^k,U^k$, 则 $X$ 子问题是凸光滑问题，对 $X$ 求矩阵导数并令其为零，
$$S-X^{-1}+U^k+\rho (X-Z^k)=0$$

这是一个关于 $X$ 的矩阵方程，可以求出满足上述矩阵方程的唯一正定的$X$ 为
$$X^{k+1}=Q\mathrm{Diag}(x_1,x_2,\cdots ,x_n)Q^{\mathrm{T}}$$

其中 $Q$ 包含矩阵 $S-\rho Z^k+U^k$ 的所有特征向量，$x_i$ 的表达式为
$$x_i=\frac{-d_i+\sqrt{d_i^2+4\rho }}{2\rho }$$

$d_i$ 为矩阵 $S-\rho Z^k+U^k$ 的第 $i$ 个特征值. 其次，固定 $X^{k+1},U^k$, 则 $Z$ 的更新为矩阵形 ${\ell }_1$ 范数的邻近算子. 最后是常规的乘子更新.

4. 矩阵分离问题
考虑矩阵分离问题 (3.8.2)：
$$\begin{array}{rl}\underset{X,S}{\min }& \Vert X{\Vert }_{\ast }+\mu \Vert S{\Vert }_1\\ s.t.& X+S=M\end{array}(\mathrm{8.6.31})$$

其中 $\Vert \cdot {\Vert }_{\ast }$ 与 $\Vert \cdot {\Vert }_1$ 分别表示矩阵核范数与 ${\ell }_1$ 范数. 引入乘子 $Y$ 作用在约束 $X+S=M$ 上，可以得到此问题的增广拉格朗日函数
$$L_{\rho }(X,S,Y)=\Vert X{\Vert }_{\ast }+\mu \Vert S{\Vert }_1+⟨Y,X+S-M⟩+\frac{\rho }{2}\Vert X+S-M{\Vert }_F^2(\mathrm{8.6.32})$$

在第 $(k+1)$ 步，交替方向乘子法分别求解关于 $X$ 和 $S$ 的子问题来更新得到 $X^{k+1}$ 和 $S^{k+1}$. 对于 $X$ 子问题，
$$\begin{array}{rl}{X}^{k+1}& =\underset{X}{\mathrm{argmin}}{L}_{\rho }(X,S^k,Y^k)\\ & =\underset{X}{\mathrm{argmin}}\left\{\Vert X{\Vert }_{\ast }+\frac{\rho }{2}{\Vert X+S^k-M+\frac{Y^k}{\rho }\Vert }_F^2\right\}\\ & =\underset{X}{\mathrm{argmin}}\left\{\frac{1}{\rho }\Vert X{\Vert }_{\ast }+\frac{1}{2}{\Vert X+S^k-M+\frac{Y^k}{\rho }\Vert }_F^2\right\}\\ & =U\mathrm{Diag}({\mathrm{prox}}_{(1/\rho )\Vert \cdot {\Vert }_1}(\sigma (A)))V^{\mathrm{T}}\end{array}$$

其中 $A=M-S^k-\frac{Y^k}{\rho },\sigma (A)$ 为 $A$ 的所有非零奇异值构成的向量并且 $U$Diag$(\sigma (A))V^{\mathrm{T}}$ 为 $A$ 的约化奇异值分解. 对于 $S$ 子问题，
$$\begin{array}{rl}{S}^{k+1}& =\underset{S}{\mathrm{argmin}}{L}_{\rho }(X^{k+1},S,Y^k)\\ & =\underset{S}{\mathrm{argmin}}\left\{\mu \Vert S{\Vert }_1+\frac{\rho }{2}{\Vert X^{k+1}+S-M+\frac{Y^k}{\rho }\Vert }_F^2\right\}\\ & ={\mathrm{prox}}_{(\mu /\rho )\Vert \cdot {\Vert }_1}\left(M-X^{k+1}-\frac{Y^k}{\rho }\right)\end{array}$$

对于乘子 $Y$, 依然使用常规更新，即
$$Y^{k+1}=Y^k+\tau \rho (X^{k+1}+S^{k+1}-M)$$

那么, 交替方向乘子法的迭代格式为
$$\begin{array}{rl}{X}^{k+1}& =U\mathrm{Diag}({\mathrm{prox}}_{(1/\rho )\Vert \cdot {\Vert }_1}(\sigma (A)))V^{\mathrm{T}}\\ S^{k+1}& ={\mathrm{prox}}_{(\mu /\rho )\Vert \cdot {\Vert }_1}\left(M-L^{k+1}-\frac{Y^k}{\rho }\right)\\ Y^{k+1}& =Y^k+\tau \rho (X^{k+1}+S^{k+1}-M)\end{array}$$

5. 全局一致性优化问题
第 8.6.1 小节介绍了全局一致性优化问题
$$\underset{x}{\min }\sum _{i=1}^N{\varphi }_i(x)$$

并给出了一个拆分方式
$$\begin{array}{rl}\underset{x_i,z}{\min }& \sum _{i=1}^N{\varphi }_i(x_i)\\ s.t.& x_i-z=0,i=1,2,\cdots ,N\end{array}$$

其增广拉格朗日函数为
$$L_{\rho }(x_1,x_2,\cdots ,x_N,z,y_1,y_2,\cdots ,y_N)=\sum _{i=1}^N{\varphi }_i(x_i)+\sum _{i=1}^N{y}_i^{\mathrm{T}}(x_i-z)+\frac{\rho }{2}\sum _{i=1}^N\Vert x_i-z{\Vert }^2$$

固定 $z^k,y_i^k$, 更新 $x_i$ 的公式为
$$x_i^{k+1}=\underset{x}{\mathrm{argmin}}\left\{{\varphi }_i(x)+\frac{\rho }{2}{\Vert x-z^k+\frac{y_i^k}{\rho }\Vert }^2\right\}(\mathrm{8.6.33})$$

在这里注意，虽然表面上看增广拉格朗日函数有 $(N+1)$ 个变量块，但本质上还是两个变量块. 这是因为在更新某 $x_i$ 时并没有利用其他 $x_i$ 的信息，所有 $x_i$ 可以看成一个整体. 相应地，所有乘子 $y_i$ 也可以看成一个整体. 迭代式 (8.6.33) 的具体计算依赖于 ${\varphi }_i$ 的形式，在一般情况下更新 $x_i$ 的表达式为
$$x_i^{k+1}={\mathrm{prox}}_{{\varphi }_i/\rho }\left(z^k-\frac{y_i^k}{\rho }\right)$$

固定 $x_i^{k+1},y_i^k$, 问题关于 $z$ 是二次函数，因此可以直接写出显式解：
$$z^{k+1}=\frac{1}{N}\sum _{i=1}^N\left(x_i^{k+1}+\frac{y_i^k}{\rho }\right)$$

综上，该问题的交替方向乘子法迭代格式为
$$\begin{array}{rl}& x_i^{k+1}={\mathrm{prox}}_{{\varphi }_i/\rho }\left(z^k-\frac{y_i^k}{\rho }\right),i=1,2,\cdots ,N\\ & z^{k+1}=\frac{1}{N}\sum _{i=1}^N\left(x_i^{k+1}+\frac{y_i^k}{\rho }\right)\\ & y_i^{k+1}={y_i}^k+\tau \rho ({x_i}^{k+1}-z^{k+1}),i=1,2,\cdots ,N\end{array}$$

6. 非凸集合上的优化问题
非凸集合上的优化问题可以表示为
$$\begin{array}{rl}\min & f(x)\\ s.t.& x\in S\end{array}(\mathrm{8.6.34})$$

其中 $f$ 是闭凸函数，但 $S$ 是非凸集合. 利用例 8.21 的技巧，对集合 $S$ 引入示性函数 $I_S(z)$ 并做拆分，可得到问题 (8.6.34) 的等价优化问题：
$$\begin{array}{rl}\min & f(x)+I_S(z)\\ s.t.& x-z=0\end{array}(\mathrm{8.6.35})$$

其增广拉格朗日函数为
$$L_{\rho }(x,z,y)=f(x)+I_S(z)+y^{\mathrm{T}}(x-z)+\frac{\rho }{2}\Vert x-z{\Vert }^2$$

由于 $f$ 是闭凸函数，固定 $z,y$ 后对 $x$ 求极小就是计算邻近算子：
$$x^{k+1}={\mathrm{prox}}_{f/\rho }\left(z^k-\frac{y^k}{\rho }\right)$$

固定 $x,y$, 对 $z$ 求极小实际上是到非凸集合上的投影问题：
$$z^{k+1}=\underset{z\in S}{\mathrm{argmin}}\frac{1}{2}{\Vert z-\left(x^{k+1}+\frac{y^k}{\rho }\right)\Vert }^2={\mathcal{P}}_S\left(x^{k+1}+\frac{y^k}{\rho }\right)$$

一般来说，由于 $S$ 是非凸集合，计算 $P_{\mathrm{s}}$ 是比较困难的 (例如不能保证存在性和唯一性), 但是当 $S$ 有特定结构时，到 $S$ 上的投影可以精确求解.
(1) 基数：如果 S = { x ∣ ∥ x ∥ 0 ⩽ c } S=\{x\mid\|x\|_0\leqslant c\} S={x∣∥x∥0​⩽c}, 其中 $\Vert \cdot {\Vert }_0$ 表示 ${\ell }_0$ 范数，即非零元素的数目，那么计算任意向量 $v$ 到 $S$ 中的投影就是保留 $v$ 分量中绝对值从大到小排列的前 $c$ 个，其余分量变成 0. 假设 $v$ 的各个分量满足
$$|v_{i_1}|⩾|v_{i_2}|⩾\cdots ⩾|v_{i_n}|$$

则投影算子可写成
( P S ( v ) ) i = { v i , i ∈ { i 1 , i 2 , ⋯   , i c } , 0 , 其他 . (\mathcal{P}_S(v))_i=\begin{cases}v_i,&i\in\{i_1,i_2,\cdots,i_c\},\\[2ex]0,&\text{其他}.\end{cases} (PS​(v))i​=⎩ ⎨ ⎧​vi​,0,​i∈{i1​,i2​,⋯,ic​},其他.​

(2) 低秩投影：当变量 $x\in {\mathbb{R}}^{m\times n}$ 是矩阵时，一种常见的非凸约束是在低秩矩阵空间中进行优化问题求解，即 S = { x ∣ r a n k ( x ) ⩽ r } S=\{x\mid \mathrm{rank}(x)\leqslant r\} S={x∣rank(x)⩽r}. 此时到低秩矩阵空间上的投影等价于对 $x$ 做截断奇异值分解. 设 $x$ 的奇异值分解为
x = ∑ i = 1 min ⁡ { m , n } σ i u i v i T x=\sum_{i=1}^{\min\{m,n\}}\sigma_iu_iv_i^\mathrm{T} x=i=1∑min{m,n}​σi​ui​viT​

其中 σ 1 ⩾ ⋯ ⩾ σ min ⁡ { m , n } ⩾ 0 \sigma_1\geqslant\cdots\geqslant\sigma_{\min\{m,n\}}\geqslant0 σ1​⩾⋯⩾σmin{m,n}​⩾0 为奇异值，$u_i,v_i$ 为对应的左右奇异向量，则 ${\mathcal{P}}_S$ 的表达式为
$${\mathcal{P}}_S(x)=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^{\mathrm{T}}$$

(3) 布尔 (Bool) 约束：如果限制变量 $x$ 只能在 { 0 , 1 } \{0,1\} {0,1} 中取值，即 S = { 0 , 1 } n S=\{0,1\}^n S={0,1}n, 容易验证 ${\mathcal{P}}_S(v)$ 就是简单地把 $v$ 的每个分量 $v_i$ 变为 0 和 1中离它更近的数，一种可能的实现方式是分别对它们做四舍五入操作.

7. 非负矩阵分解和补全

非负矩阵分解和补全可以看作非负矩阵分解问题和低秩矩阵补全问题的结合，即已知一个非负矩阵的部分元素，求其非负分解.
假设我们有从一个非负、秩为 $r$ 的矩阵 $M\in {\mathbb{R}}^{m\times n}$ 中采样的部分元素 M i , j , ( i , j ) ∈ Ω ⊂ { 1 , 2 , ⋯   , m } × { 1 , 2 , ⋯   , n } M_{i,j},(i,j)\in\Omega\subset\{1,2,\cdots,m\}\times\{1,2,\cdots,n\} Mi,j​,(i,j)∈Ω⊂{1,2,⋯,m}×{1,2,⋯,n}, 目标是要找到非负矩阵 $X\in {\mathbb{R}}^{m\times q},Y\in {\mathbb{R}}^{q\times n}$ 使得 $\Vert M-XY{\Vert }_F^2$ 极小. 一般地，因为数据和应用问题的不同，$q$ 可以等于、小于或者大于 $r$. 定义矩阵 $P\in {\mathbb{R}}^{m\times n}$：
$$P_{ij}=\{\begin{array}{ll}1,& (i,j)\in \Omega \\ 0,& \text{其他}\end{array}$$

则非负矩阵分解和补全问题可以写成如下形式：
$$\begin{array}{rl}\underset{X,Y}{\min }& \Vert P\odot (XY-M){\Vert }_F^2\\ s.t.& X_{ij}⩾0,Y_{ij}⩾0,\forall i,j\end{array}$$

注意，这个问题是非凸的．
为了利用交替方向乘子法的优势，我们考虑如下等价形式：
$$\begin{array}{rl}\underset{U,V,X,Y,Z}{\min }& \frac{1}{2}\Vert XY-Z{\Vert }_F^2\\ s.t.& X=U,Y=V\\ & U⩾0,V⩾0\\ & P\odot (Z-M)=0\end{array}$$

对于等式约束 $X=U$ 和 $Y=V$ 分别引入拉格朗日乘子 $\Lambda$ 和 $\Pi$, 并保留非负约束 $U⩾0,V⩾0$ 和观测值约束 $P\odot (Z-M)=0$, 则增广拉格朗日函数为
$$\begin{array}{rl}& L_{\alpha ,\beta }(X,Y,Z,U,V,\Lambda ,\Pi )\\ =& \frac{1}{2}\Vert XY-Z{\Vert }_F^2+⟨\Lambda ,X-U⟩+⟨\Pi ,Y-V⟩+\frac{\alpha }{2}\Vert X-U{\Vert }_F^2+\frac{\beta }{2}\Vert Y-V{\Vert }_F^2\end{array}(\mathrm{8.6.36})$$

应用交替方向乘子法，可得
$$\begin{array}{rl}& X^{k+1}=\underset{X}{\mathrm{argmin}}{L}_{\alpha ,\beta }(X,Y^k,Z^k,U^k,V^k,{\Lambda }^k,{\Pi }^k)\\ & Y^{k+1}=\underset{Y}{\mathrm{argmin}}{L}_{\alpha ,\beta }(X^{k+1},Y,Z^k,U^k,V^k,{\Lambda }^k,{\Pi }^k)\\ & Z^{k+1}=\underset{P\odot (Z-M)=0}{\mathrm{argmin}}{L}_{\alpha ,\beta }(X^{k+1},Y^{k+1},Z,U^k,V^k,{\Lambda }^k,{\Pi }^k)\\ & U^{k+1}=\underset{U⩾0}{\mathrm{argmin}}{L}_{\alpha ,\beta }(X^{k+1},Y^{k+1},Z^{k+1},U,V^k,{\Lambda }^k,{\Pi }^k)\\ & V^{k+1}=\underset{V⩾0}{\mathrm{argmin}}{L}_{\alpha ,\beta }(X^{k+1},Y^{k+1},Z^{k+1},U^{k+1},V,{\Lambda }^k,{\Pi }^k)\\ & {\Lambda }^{k+1}={\Lambda }^k+\tau \alpha (X^{k+1}-U^{k+1})\\ & {\Pi }^{k+1}={\Pi }^k+\tau \beta (Y^{k+1}-V^{k+1})\end{array}$$

将子问题求解，可得
$$\begin{array}{rl}& X^{k+1}=(Z^k(Y^k{)}^{\mathrm{T}}+\alpha U^k-A^k)(Y^k(Y^k{)}^{\mathrm{T}}+\alpha I{)}^{-1}\\ & Y^{k+1}=((X^{k+1}{)}^{\mathrm{T}}{X}^{k+1}+\beta I{)}^{-1}((X^{k+1}{)}^{\mathrm{T}}{Z}^k+\beta V^k-{\Pi }^k)\\ & Z^{k+1}=X^{k+1}{Y}^{k+1}+P\odot (M-X^{k+1}{Y}^{k+1})\\ & U^{k+1}={\mathcal{P}}_{+}\left(X^{k+1}+\frac{{\Lambda }^k}{\alpha }\right)\\ & V^{k+1}={\mathcal{P}}_{+}\left(Y^{k+1}+\frac{{\Pi }^k}{\beta }\right)\\ & {\Lambda }^{k+1}={\Lambda }^k+\tau \alpha (X^{k+1}-U^{k+1})\\ & P{i}^{k+1}=P{i}^k+\tau \beta (Y^{k+1}-V^{k+1})\end{array}$$

其中 ( P + ( A ) ) i j = max ⁡ { a i j , 0 } (\mathcal{P}_{+}(A))_{ij}=\max\{a_{ij},0\} (P+​(A))ij​=max{aij​,0}. 注意, 该格式为多块交替方向乘子法, 其收敛性可能要要较强假设.

8. 多块交替方向乘子法的反例
这里给出一个多块交替方向乘子法的例子，并且从数值上说明若直接采用格式 (8.6.25)，则算法未必收敛．
考虑最优化问题
$$\begin{array}{rl}\min & 0\\ s.t.& A_1{x}_1+A_2{x}_2+A_3{x}_3=0\end{array}(\mathrm{8.6.37})$$

其中 $A_i\in {\mathbb{R}}^3,i=1,2,3$ 为三维空间中的非零向量，$x_i\in \mathbb{R},i=1,2,3$ 是自变量. 问题 (8.6.37) 实际上就是求解三维空间中的线性方程组，若 $A_1,A_2,A_3$ 之间线性无关，则问题 (8.6.37) 只有零解. 此时容易计算出最优解对应的乘子为 $y=(0,0,0{)}^{\mathrm{T}}$.
现在推导多块交替方向乘子法的格式. 问题 (8.6.37) 的增广拉格朗日函数为
$$L_{\rho }(x,y)=0+y^{\mathrm{T}}(A_1{x}_1+A_2{x}_2+A_3{x}_3)+\frac{\rho }{2}\Vert A_1{x}_1+A_2{x}_2+A_3{x}_3{\Vert }^2$$

当固定 $x_2,x_3,y$ 时，对 $x_1$ 求最小可推出
$$A_1^{\mathrm{T}}y+\rho A_1^{\mathrm{T}}(A_1{x}_1+A_2{x}_2+A_3{x}_3)=0$$

整理可得
$$x_1=-\frac{1}{\Vert A_1{\Vert }^2}{A}_1^{\mathrm{T}}\left(\frac{y}{\rho }+A_2{x}_2+A_3{x}_3\right)$$

可类似地计算 $x_2,x_3$ 的表达式，因此多块交替方向乘子法的迭代格式可以写为
$$\begin{array}{rl}& x_1^{k+1}=-\frac{1}{\Vert A_1{\Vert }^2}{A}_1^{\mathrm{T}}\left(\frac{y^k}{\rho }+A_2{x}_2^k+A_3{x}_3^k\right)\\ & x_2^{k+1}=-\frac{1}{\Vert A_2{\Vert }^2}{A}_2^{\mathrm{T}}\left(\frac{y^k}{\rho }+A_1{x}_1^{k+1}+A_3{x}_3^k\right)\\ & x_3^{k+1}=-\frac{1}{\Vert A_3{\Vert }^2}{A}_3^{\mathrm{T}}\left(\frac{y^k}{\rho }+A_1{x}_1^{k+1}+A_2{x}_2^{k+1}\right)\\ & y^{k+1}=y^k+\rho (A_1{x}_1^{k+1}+A_2{x}_2^{k+1}+A_3{x}_3^{k+1})\end{array}(\mathrm{8.6.38})$$

对此问题而言，罚因子取不同值仅仅是将乘子 $y^k$ 缩放了常数倍，所以罚因子 $\rho$ 的任意取法 (包括动态调节) 都是等价的.
在数值实验中我们不妨取 $\rho =1$.
格式 (8.6.38) 的收敛性与 $A_i,i=1,2,3$ 的选取有关. 为了方便，令 $A=[A_1,A_2,A_3]$, 以及 $x=(x_1,x_2,x_3{)}^{\mathrm{T}}$, 并选取 $A$ 为
$$\tilde{A}=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 2\\ 1& 2& 2\end{array}\right]\text{或}\hat{A}=\left[\begin{array}{ccc}1& 1& 2\\ 0& 1& 1\\ 0& 0& 1\end{array}\right]$$

在迭代中，自变量初值选为(1,1,1), 乘子选为(0,0,0).

import numpy as np
import matplotlib.pyplot as plt

#A_tilde = np.array([[1, 1, 1], [1, 1, 2], [1, 2, 2]])
A1_tilde = np.array([1, 1, 1])
A2_tilde = np.array([1, 1, 2])
A3_tilde = np.array([1, 2, 2])
# A_hat = np.array([[1, 1, 2], [0, 1, 1], [0, 0, 1]])
A1_hat = np.array([1, 1, 2])
A2_hat = np.array([0, 1, 1])
A3_hat = np.array([0, 0, 1])

# Define the norm square of Ai
norm_A1_tilde_sq = np.linalg.norm(A1_tilde)**2
norm_A2_tilde_sq = np.linalg.norm(A2_tilde)**2
norm_A3_tilde_sq = np.linalg.norm(A3_tilde)**2

norm_A1_hat_sq = np.linalg.norm(A1_hat)**2
norm_A2_hat_sq = np.linalg.norm(A2_hat)**2
norm_A3_hat_sq = np.linalg.norm(A3_hat)**2

def admm(A1, A2, A3, norm_A1_sq, norm_A2_sq, norm_A3_sq, num_iterations=100):
    x1, x2, x3 = 1, 1, 1
    y = np.zeros(3)
    rho = 1
    x_norms = []
    y_norms = []
    for k in range(num_iterations):
        x1 = -1 / norm_A1_sq * A1.T.dot(y / rho + A2 * x2 + A3 * x3)
        x2 = -1 / norm_A2_sq * A2.T.dot(y / rho + A1 * x1 + A3 * x3)
        x3 = -1 / norm_A3_sq * A3.T.dot(y / rho + A1 * x1 + A2 * x2)
        y = y + rho * (A1 * x1 + A2 * x2 + A3 * x3)   
        x_norms.append(np.linalg.norm([x1, x2, x3]))
        y_norms.append(np.linalg.norm(y))
    return x_norms, y_norms

x_norms_tilde, y_norms_tilde = admm(A1_tilde, A2_tilde, A3_tilde, norm_A1_tilde_sq, norm_A2_tilde_sq, norm_A3_tilde_sq)
x_norms_hat, y_norms_hat = admm(A1_hat, A2_hat, A3_hat, norm_A1_hat_sq, norm_A2_hat_sq, norm_A3_hat_sq)

plt.figure(figsize=(6, 8))
# tilde A
plt.subplot(2, 1, 1)
plt.plot(x_norms_tilde, label='$||x||$')
plt.plot(y_norms_tilde, linestyle='--', label='$||y||$')
plt.title('ADMM Convergence for $\widetilde{A}$')
plt.xlabel('Iterations')
plt.xlim(0)
plt.legend()

# hat A
plt.subplot(2, 1, 2)
plt.plot(x_norms_hat, label='$||x||$')
plt.plot(y_norms_hat, linestyle='--', label='$||y||$')
plt.title('ADMM Convergence for $\widehat{A}$')
plt.xlabel('Iterations')
plt.xlim(0)
plt.yscale('log')
plt.legend()
plt.tight_layout()

上图记录了在两种不同 $A$ 条件下，$x$ 和 $y$ 的 ${\ell }_2$ 范数随迭代的变化过程. 可以看到当 $A=\tilde{A}$ 时数值结果表明迭代是发散的，而当 $A=\hat{A}$ 时数值结果表明迭代是收敛的. 另一个比较有趣的观察是，在下子图 中 $\Vert x\Vert$ 与 $\Vert y\Vert$ 并不是单调下降的，而是在下降的同时有规律地振荡.

参考教材《最优化：建模、算法与理论》