ADMM原理及应用

1. ADMM原理

1.1. 数学形式

我们要解的优化问题长这样：
$$\underset{x\in {\mathbf{R}}^n,z\in {\mathbf{R}}^m}{\min }f(x)+g(z)\text{subject to}Ax+Bz=c.$$
这意味着我们想同时让$f(x)$ 和 $g(z)$尽可能小，但它们又要满足一个线性约束 $Ax+Bz=c$。

• $x$ 和 $z$分别是不同的变量，
• $f$ 和$g$都是凸函数（这保证了算法更容易收敛），
• $A$ 和 $B$ 是已知的矩阵，
• $c$ 是已知的常量。

为什么要分成 $x$ 和 $z$ 两组变量？通常是因为 $f$ 和$g$ 可能各自有特殊结构，例如一个是稀疏正则项（$L_1$ 范数），另一个可能是平方和之类的简单函数。要是把它们混在一起难以统一求解，就可以“拆开”来做。

1.2. 传统“乘子法”和它的不足

在最经典的拉格朗日乘子法里，我们会先把约束放进一个“增强过”的目标函数里，称为“增广拉格朗日函数”（augmented Lagrangian），它一般长这样：
$$L_{\rho }(x,z,y)=f(x)+g(z)+y^T(Ax+Bz-c)+\frac{\rho }{2}\Vert Ax+Bz-c{\Vert }_2^2,$$
其中 $y$ 是所谓的“对偶变量”或者“拉格朗日乘子”，$\rho >0$ 是个参数（它决定惩罚力度）。

在传统乘子法里，每次迭代要同时对 $x$ 和 $z$ 都做一个“联合最小化”（joint minimization）：
$$(x^{k+1},z^{k+1})=\arg \underset{x,z}{\min }{L}_{\rho }(x,z,y^k),$$
然后再更新对偶变量
$$y^{k+1}=y^k+\rho (A{x}^{k+1}+B{z}^{k+1}-c).$$
这样做当然可以，但如果 $f$ 和 $g$的形式比较复杂，或者维度较大，那这个“联合最小化”就不好算（可能很耗时，或者甚至求不出来）。

1.3. ADMM 的核心思想：分步做

ADMM（Alternating Direction Method of Multipliers，“交替方向乘子法”）最主要的特色就是不再让 $x$ 和 $z$ 同时做大的联合求解，而是“先算 $x$，再算 $z$”的交替方式。它的三步更新如下：

1. $x$-更新：固定住旧的 $z^k$ 和 $y^k$，只对 $x$​ 做一个最优更新：
   $$x^{k+1}=\arg \underset{x}{\min }{L}_{\rho }(x,z^k,y^k).$$
   由于 $z^k、y^k$ 不变，这一步就是在一个比较“简化了”的函数里找最优 $x$。

2. $z$-更新：拿到更新后的$x^{k+1}$，再固定它和$y^k$，对 $z$​ 做最优更新：
   $$z^{k+1}=\arg \underset{z}{\min }{L}_{\rho }(x^{k+1},z,y^k).$$

3. 对偶变量 $y$ 更新：有了新的 $x^{k+1}$ 和 $z^{k+1}$，再更新对偶变量 $y$​：
   $$y^{k+1}=y^k+\rho (A{x}^{k+1}+B{z}^{k+1}-c).$$

这样一来，每一步都只是对一个变量做最小化，问题规模往往更小，如果 $f(x)$ 和 $g(z)$ 还是那种可以分开处理的“友好”函数，求解起来也更容易、更快。

1.4. Scaled Form of ADMM

在 Scaled Form 里，我们把对偶变量 换一种等价的表示。为了方便，我们先定义一个所谓的残差（residual）：
$$r=Ax+Bz-c.$$
同时，定义“缩放后的对偶变量”（scaled dual variable）$u$​：
$$u=\frac{1}{\rho }y.$$
这样，原先的项 $y^{T}r+\frac{\rho }{2}\Vert r{\Vert }^2$可以用 $u$​ 来重写成
$$\underset{\text{重新打包}}{\underbrace{(\rho /2)\Vert r+u{\Vert }_2^2}}-\underset{\text{校正}}{\underbrace{(\rho /2)\Vert u{\Vert }_2^2}}.$$
这个结论的得来类似于高中学到的凑平方项。

“用 $u=\frac{1}{\rho }y$ 替换后，线性和二次项可以组合到一个“$\frac{\rho }{2}\Vert r+u{\Vert }^2$”形式里，看起来更整洁。”

在这个新的记号下，ADMM 的迭代过程可以写成（省去常数项）：

1. $x$-更新：
   $$x^{k+1}=\arg \underset{x}{\min }(f(x)+\frac{\rho }{2}\Vert Ax+B{z}^k-c+u^k{\Vert }_2^2).$$

2. $z$-更新：
   $$z^{k+1}=\arg \underset{z}{\min }(g(z)+\frac{\rho }{2}\Vert A{x}^{k+1}+Bz-c+u^k{\Vert }_2^2).$$

3. $u$-更新：
   $$u^{k+1}=u^k+(A{x}^{k+1}+B{z}^{k+1}-c).$$

你会发现，现在对偶更新变得很简单，直接是对旧的 $u$ 加上“残差”$r=Ax+Bz-c$​。而在 unscaled form 里则是
$$y^{k+1}=y^k+\rho r^k.⟺u^{k+1}=u^k+r^k,\text{因为 }u=\frac{1}{\rho }y.$$
也就是说，两个形式做的事情完全一样，只是在对偶变量上做一个因子 $\frac{1}{\rho }$ 的缩放。

1.5. 迭代过程中主要检查的两大残差

两个残差（residual）：

1. 主（primal）残差
   $$r^{k+1}=A{x}^{k+1}+B{z}^{k+1}-c.$$
   这是用来度量“原约束  $Ax+Bz=c$”在迭代第 $k+1$ 步时的偏差。若$r^{k+1}$ 越接近 0，说明越接近可行。

2. 对偶（dual）残差
   $$s^{k+1}=\rho A^{T}B(z^{k+1}-z^k).$$
   这是用来度量对偶可行性。若$s^{k+1}$ 越接近 0，说明越接近对偶可行。

在 ADMM 中，你会看到在每次迭代完 $x^{k+1},z^{k+1}$​ 后，会“顺手”计算这两个残差，用来判断收敛程度。

1.6. 怎么设置停止准则(Stopping Criteria)？

一个常见且实用的做法就是直接对主残差和对偶残差设置阈值：
$$\Vert r^k{\Vert }_2\le {\epsilon }_{\text{pri}}\text{and}\Vert s^k{\Vert }_2\le {\epsilon }_{\text{dual}},$$

• ${\epsilon }_{\text{pri}}>0$（主可行性余量）
• ${\epsilon }_{\text{dual}}>0$（对偶可行性余量）

只要这两个残差都小于各自阈值，就认定收敛。

但是，上述设置也会存在如下问题：

• 如果问题中 $\Vert c\Vert$、$\Vert A{x}^k\Vert$、$\Vert B{z}^k\Vert$ 值很大，仅仅用一个小绝对值判断误差会显得“吹毛求疵”，或者数值上不稳定；
• 如果问题中变量本身特别小，仅用相对误差也可能不够；

所以可以采用下面一个典型设置：
ε pri = p ε abs + ε rel max ⁡ { ∥ A x k ∥ 2 , ∥ B z k ∥ 2 , ∥ c ∥ 2 } , ε dual = n ε abs + ε rel ∥ A T y k ∥ 2 , \varepsilon_{\text{pri}} =\sqrt{p}\,\varepsilon_{\text{abs}} \;+\; \varepsilon_{\text{rel}}\, \max\,\{\,\|A x^k\|_2,\;\|B z^k\|_2,\;\|c\|_2\},\\ \varepsilon_{\text{dual}} =\sqrt{n}\,\varepsilon_{\text{abs}} \;+\; \varepsilon_{\text{rel}}\, \|A^T y^k\|_2, εpri​=p ​εabs​+εrel​max{∥Axk∥2​,∥Bzk∥2​,∥c∥2​},εdual​=n ​εabs​+εrel​∥ATyk∥2​,

• ${\epsilon }_{\text{abs}}$（absolute tolerance）是绝对误差上限，
• ${\epsilon }_{\text{rel}}$（relative tolerance）是相对误差系数，
• $\sqrt{p}$ 和$\sqrt{n}$ 分别考虑了残差向量所在的维度（主问题中是 ${\mathbf{R}}^p$，对偶中是 ${\mathbf{R}}^n$）。
• 通常在应用中，会选${\epsilon }_{\text{rel}}=10^{-3}$或 $10^{-4}$​ 之类，视需求和数值规模而定。

1.7. 自适应调整罚参数 $\rho$（又称“变步长”技巧）

在标准 ADMM 中，$\rho$（也叫增广拉格朗日里的“罚参数”）一般是个固定常数。它控制着算法对原约束违背 $r^k=A{x}^k+B{z}^k-c$ 的惩罚力度：

• $\rho$ 大：更严厉地惩罚主可行性的违背，因此更容易让$r^k$ 变得较小，但对偶残差 $s^k$可能相对变大。
• $\rho$ 小：相对减少对主可行性的惩罚，往往会导致$r^k$ 可能大一点，但好处是$s^k$ 会小一些。

如果在迭代过程中主残差和对偶残差差距太大，会让 ADMM 收敛变慢。自适应调节$\rho$ 的目标就是让主、对偶残差同时保持在一个“类似量级”，这样往往能得到更均衡、更快的收敛。

常用方案是：在第 $k+1$ 步开始前，根据上一轮迭代的残差 $\Vert r^k\Vert$ 和 $\Vert s^k\Vert$ 的大小比较来决定 ${\rho }^{k+1}$​ 的取值：
$${\rho }^{k+1}:=\{\begin{array}{ll}{\tau }_{\text{incr}}{\rho }^k,& \text{如果 }\Vert r^k{\Vert }_2>\mu \Vert s^k{\Vert }_2,\\ {\rho }^k/{\tau }_{\text{decr}},& \text{如果 }\Vert s^k{\Vert }_2>\mu \Vert r^k{\Vert }_2,\\ {\rho }^k,& \text{否则.}\end{array}$$
这里：

• $\mu >1$ 通常是一个常数（比如 $\mu =10$），用来判断两种残差谁大很多；
• ${\tau }_{\text{incr}}>1$和 ${\tau }_{\text{decr}}>1$ 是用来放大或缩小 $\rho$ 的倍数（如 ${\tau }_{\text{incr}}=2,{\tau }_{\text{decr}}=2$）。

根据这个公式：

• 当“主残差大”（ $\Vert r^k\Vert >\mu \Vert s^k\Vert$) 时，就把$\rho$ 增大 ${\tau }_{\text{incr}}$r 倍，以便在后续迭代中对主可行性的违背更严厉些，让 $r^{k+1}$ 尽可能降下来；
• 当“对偶残差大” ($\Vert s^k\Vert >\mu \Vert r^k\Vert$) 时，就把 $\rho$ 减小${\tau }_{\text{decr}}$ 倍，让对偶可行性得到更多纠正，力图让 $s^{k+1}$ 得到控制；
• 否则就保持 $\rho$ 不变，说明主、对偶残差的量级相对平衡，没必要动 $\rho$。

这样做的核心思路是：“让 $\Vert r^k\Vert$ 和 $\Vert s^k\Vert$ 保持在一个大致相当的范围”。当它们差距太大时，就通过调节 $\rho$ 来“补救”。

需要注意的细节

Scaled Form 中的对偶变量要重缩放
如果你用的是 Scaled Form（即 $u^k=\frac{1}{\rho }{y}^k$ 这种），当 $\rho$ 改变时，就要相应地更新 $u^k$​：
$$u^{k+1}=\frac{1}{{\rho }^{k+1}}{y}^{k+1}.$$
在代码实现里，这往往意味着：如果你决定把 ${\rho }^k$ 改为${\rho }^{k+1}=\frac{1}{2}{\rho }^k$，就要把 $u^k$ 乘以 2（因为 $\frac{1}{\rho }$ 这一项翻倍了）。
否则会导致对偶变量与新的$\rho$​ 不一致、破坏算法正确性。

1.8. Over-relaxation (过松弛)

在 ADMM 的 $z$-update 和 $y$-update 里，常常会出现表达式 $A{x}^{k+1}$ 。Over-relaxation (过松弛) 指的是用下面这种线性组合替代$A{x}^{k+1}$：
$${\alpha }^{k}A{x}^{k+1}-(1-{\alpha }^k)(B{z}^k-c),$$
其中 ${\alpha }^k\in (0,2)$ 是一个松弛因子。当 ${\alpha }^k>1$ 时称为 over-relaxation (过松弛)；当 ${\alpha }^k<1$ 称为 under-relaxation (欠松弛)。

实验现象：

• 在许多应用里，选择 ${\alpha }^k$ 稍大于 1（比如 1.5 到 1.8 之间）可以加快收敛。它有点类似在某些迭代算法里“加一点动量”或“加速项”的感觉。
• 但过度 over-relaxation 又可能导致不稳定，所以一般不会取太大。

2. ADMM应用

2.1. ADMM求解最小绝对偏差（Least Absolute Deviations）

2.1.1. 数学形式

在回归或拟合问题中，常见的目标是最小化
$$\Vert Ax-b{\Vert }_2^2,$$
即最小二乘（Least Squares）。但有时候数据中会有较大的“离群点”（outliers），这时最小二乘可能会被几个特别极端的数据“拉偏”，并不能很好地反映整体趋势。最小绝对偏差（Least Absolute Deviations） 的想法，是用
$$\Vert Ax-b{\Vert }_1$$
来衡量误差。因为 $L_1$-范数对极端值的敏感程度不如二范数大，从而得到更鲁棒的（robust）结果。

2.1.2. 将问题转换成ADMM形式

引入变量 $z$​，令
$$z=Ax-b$$
于是原问题就可以写成
$$\begin{array}{rl}& \underset{x,z}{\min }\Vert z{\Vert }_1,\\ & \text{subject to}Ax-z=b.\end{array}$$
为了使用 ADMM，需要将目标函数写成 $f(x)+g(z)$ 加一个线性约束，这里我们定义

• $f(x)=0$，它“管”着变量 $x$ 但没有额外的目标值（因为真正的目标全在 $\Vert z{\Vert }_1$ 里）；
• $g(z)=\Vert z{\Vert }_1$，这是要最小化的部分；
• 约束 $Ax-z=b$ 则成为 ADMM 中要处理的那条等式。

2.1.3. ADMM迭代步骤

1. $x$​-update

$$\begin{array}{rl}{x}^{k+1}& =\arg \underset{x}{\min }f(x)+\frac{\rho }{2}\Vert Ax-b-z^k+u^k{\Vert }_2^2\\ & =\arg \underset{x}{\min }0+\frac{\rho }{2}\Vert Ax-b-z^k+u^k{\Vert }_2^2\\ & =\arg \underset{x}{\min }\Vert Ax-b-z^k+u^k{\Vert }_2^2.\end{array}$$

这基本上就是一个最小二乘问题：
$$\underset{x}{\min }\Vert Ax-(b+z^k-u^k){\Vert }_2^2.$$
如果 $A^{T}A$​ 可逆，那么可以直接用
$$x^{k+1}=(A^{T}A{)}^{-1}{A}^T(b+z^k-u^k).$$
这一步也可以用因式分解（比如对 $A^{T}A$ 做一次 Cholesky 分解等【见附录1】），后续迭代中重复使用，从而加速计算。做完Cholesky 分解后得到的$x$更新公式如下
$$x^{k+1}=R^{-1}(R^T{)}^{-1}{A}^T(b+z^k-u^k).$$

2. $z$​-update

$$z^{k+1}=\arg \underset{z}{\min }\Vert z{\Vert }_1+\frac{\rho }{2}\Vert A{x}^{k+1}-b-z+u^k{\Vert }_2^2.$$

这个子问题的解可以用**软阈值（soft thresholding）**算出来,即
$$z^{k+1}={\mathrm{S}}_{1/\rho }(A{x}^{k+1}-b+u^k),$$
其中 ${\mathrm{S}}_{\alpha }(\cdot )$是元素逐个的软阈值算子（对每个分量做 $z_i=\text{sign}(w_i)\max (|w_i|-\alpha ,0)$​)。关于软阈值算子的详细介绍，可以参考我这篇文章。

3. $u$​-update

$$u^{k+1}=u^k+(A{x}^{k+1}-b-z^{k+1}).$$

这一步是标准的拉格朗日乘子更新，会在迭代中不断逼近约束 $Ax-z=b$。

注意：

在ADMM介绍中，我们也提到，在 ADMM 的 $z$-update 和 $u$-update 里，常常会出现表达式 $A{x}^{k+1}$。Over-relaxation (过松弛) 指的是用下面这种线性组合替代 $A{x}^{k+1}$：
$${\alpha }^{k}A{x}^{k+1}-(1-{\alpha }^k)(B{z}^k-c),$$
其中 ${\alpha }^k\in (0,2)$ 是一个松弛因子。当 ${\alpha }^k>1$ 时称为 over-relaxation (过松弛)；当 ${\alpha }^k<1$ 称为 under-relaxation (欠松弛)。在许多应用里，选择 ${\alpha }^k$ 稍大于 1（比如 1.5 到 1.8 之间）可以加快收敛。

如果采用这一替代，那么$z$和$u$的更新就变为
$$\widehat{A{x}^{k+1}}={\alpha }^{k}A{x}^{k+1}-(1-{\alpha }^k)(B{z}^k-c),$$

$$z^{k+1}={\mathrm{S}}_{1/\rho }(\widehat{A{x}^{k+1}}-b+u^k),$$

$$u^{k+1}=u^k+(\widehat{A{x}^{k+1}}-b-z^{k+1}).$$

2.1.4. 停止准则

首先计算两个残差：

1. 主残差

$$\begin{array}{rl}{r}^{k+1}& =A{x}^{k+1}+B{z}^{k+1}-c\\ & =A{x}^{k+1}-z^{k+1}-b.\end{array}$$

2. 对偶残差
   $$\begin{array}{rl}{s}^{k+1}& =\rho A^{T}B(z^{k+1}-z^k)\\ & =-\rho A^T(z^{k+1}-z^k).\end{array}$$

然后根据停止准则的典型设置，有
ε pri = p ε abs + ε rel max ⁡ { ∥ A x k ∥ 2 , ∥ B z k ∥ 2 , ∥ c ∥ 2 } , ε dual = n ε abs + ε rel ∥ A T y k ∥ 2 , \varepsilon_{\text{pri}} =\sqrt{p}\,\varepsilon_{\text{abs}} \;+\; \varepsilon_{\text{rel}}\, \max\,\{\,\|A x^k\|_2,\;\|B z^k\|_2,\;\|c\|_2\},\\ \varepsilon_{\text{dual}} =\sqrt{n}\,\varepsilon_{\text{abs}} \;+\; \varepsilon_{\text{rel}}\, \|A^T y^k\|_2, εpri​=p ​εabs​+εrel​max{∥Axk∥2​,∥Bzk∥2​,∥c∥2​},εdual​=n ​εabs​+εrel​∥ATyk∥2​,
那么当
$$\Vert r^k{\Vert }_2\le {\epsilon }_{\text{pri}}\text{and}\Vert s^k{\Vert }_2\le {\epsilon }_{\text{dual}},$$
满足时，就可以输出结果了。

2.1.5 Matlab程序和例子

ADMM_lad函数如下

function [x, history] = ADMM_lad(A, b, rho, alpha)
% ADMM_lad   Least Absolute Deviations fitting via ADMM
% 
% [x, history] = ADMM_lad(A, b, rho, alpha)
%
% Solves the following problem via ADMM:
%   minimize ||Ax - b||_1
%
% Input:
%   A      - m x n matrix of coefficients
%   b      - m x 1 vector of observations
%   rho    - augmented Lagrangian parameter
%   alpha  - over-relaxation parameter (1.0 <= alpha <= 1.8, typical value: 1.5)
%
% Output:
%   x       - solution vector of size n x 1
%   history - structure containing iteration details:
%             - objective value at each iteration
%             - primal and dual residual norms
%             - tolerances for convergence criteria
%
% Algorithm Overview:
%   The goal is to solve the Least Absolute Deviations (LAD) problem:
%       minimize ||Ax - b||_1,
%   which is robust to outliers compared to the least squares approach. This is achieved
%   by introducing an auxiliary variable z and reformulating the problem as:
%       minimize ||z||_1,
%       subject to Ax - z = b.
%
%   Using the Alternating Direction Method of Multipliers (ADMM), the augmented Lagrangian is:
%       L(x, z, u) = ||z||_1 + (rho / 2) * ||Ax - z - b + u||_2^2,
%   where u is the scaled dual variable. The optimization proceeds in the following steps:
%     1. x-update: Solve a least squares problem for x:
%          x^{k+1} = argmin_x (1/2)||Ax - b - z^k + u^k||_2^2.
%          This step is implemented efficiently using the Cholesky factorization of A^T A.
%     2. z-update: Use the soft-thresholding operator to minimize ||z||_1:
%          z^{k+1} = S_{1/\rho}(Ax^{k+1} - b + u^k),
%          where S_{\kappa}(\cdot) is the soft-thresholding operator.
%     3. u-update: Update the dual variable to enforce the constraint Ax - z = b:
%          u^{k+1} = u^k + (Ax^{k+1} - z^{k+1} - b).
%   These steps are repeated until convergence criteria based on primal and dual residuals are met.
%
% Reference:
%   Boyd et al., "Distributed Optimization and Statistical Learning via the ADMM" (2011)%% Constants and initializationt_start = tic;               % Start timerMAX_ITER = 1000;             % Maximum number of iterationsABSTOL   = 1e-4;             % Absolute tolerance for convergenceRELTOL   = 1e-2;             % Relative tolerance for convergence[m, n] = size(A);            % Dimensions of the problem% Initialize variablesx = zeros(n, 1);             % Optimization variable xz = zeros(m, 1);             % Auxiliary variable zu = zeros(m, 1);             % Dual variable (scaled Lagrange multiplier)% Precompute Cholesky factorization of A^T A for efficient least squaresR = chol(A' * A);% Print table header for iteration logsfprintf('%-5s %-12s %-12s %-12s %-12s %-12s\n',...'Iter', 'r_norm', 'eps_pri', 's_norm', 'eps_dual', 'Objective');%% ADMM iterationsfor k = 1:MAX_ITER% Step 1: x-update (solve least squares problem)% Solve (A^T A) x = A^T (b + z - u)x = R \ (R' \ (A' * (b + z - u)));% Step 2: z-update (soft-thresholding)z_old = z;  % Save previous value of zAx_hat = alpha * A * x + (1 - alpha) * (z_old + b);  % Over-relaxation stepz = shrinkage(Ax_hat - b + u, 1 / rho);  % Soft-thresholding operator% Step 3: u-update (dual variable update)u = u + (Ax_hat - z - b);  % Update scaled Lagrange multiplier% Diagnostics and reportinghistory.objval(k)  = norm(z, 1);  % Objective function valuehistory.r_norm(k)  = norm(A * x - z - b);  % Primal residual normhistory.s_norm(k)  = norm(-rho * A' * (z - z_old));  % Dual residual norm% Convergence toleranceshistory.eps_pri(k) = sqrt(m) * ABSTOL + RELTOL * max([norm(A * x), norm(z), norm(b)]);history.eps_dual(k) = sqrt(n) * ABSTOL + RELTOL * norm(rho * A' * u);% Print iteration detailsfprintf('%-5d %-12.4f %-12.4f %-12.4f %-12.4f %-12.2f\n', k, ...history.r_norm(k), history.eps_pri(k), ...history.s_norm(k), history.eps_dual(k), history.objval(k));% Check for convergenceif (history.r_norm(k) < history.eps_pri(k) && ...history.s_norm(k) < history.eps_dual(k))break;endend% Print total computation timefprintf('Total time: %.2f seconds\n', toc(t_start));
endfunction y = shrinkage(a, kappa)
% Soft-thresholding operatory = sign(a) .* max(abs(a) - kappa, 0);
end

一个简单的调用例子

clear; close all; clc;% Example of ADMM for Least Absolute Deviations (LAD) fitting
rng(0);  % Set random seed for reproducibility% Problem dimensions
m = 500;  % Number of rows (observations)
n = 100;  % Number of columns (features)% Generate random data
A = randn(m, n);          % Coefficient matrix
x0 = 10 * randn(n, 1);    % True signal
b = A * x0;               % Observed data without noise% Add large outliers to simulate real-world noisy data
idx = randsample(m, ceil(m/50));  % Randomly pick ~2% of rows
b(idx) = b(idx) + 1e2 * randn(size(idx));  % Add large noise% Call ADMM LAD solver
rho = 1;   % ADMM parameter
alpha = 1; % Over-relaxation parameter
[x, history] = ADMM_lad(A, b, rho, alpha);% Plot results
figure('Position', [500, 200, 800, 600]);  % Set figure size% Plot original signal and reconstructed signal
subplot(2, 1, 1);
plot(1:n, x0, 'b-', 'LineWidth', 1.5); hold on;
plot(1:n, x, 'r--', 'LineWidth', 1.5);
xlabel('Index', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Value', 'FontSize', 12, 'FontWeight', 'bold');
title('Original Signal vs. Reconstructed Signal', 'FontSize', 14, 'FontWeight', 'bold');
legend({'True Signal (x_0)', 'Reconstructed Signal (x)'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);  % Adjust axis font size and line width% Plot convergence history
subplot(2, 1, 2);
semilogy(1:length(history.objval), history.objval, 'k-', 'LineWidth', 1.5); hold on;
semilogy(1:length(history.r_norm), history.r_norm, 'r-', 'LineWidth', 1.5);
semilogy(1:length(history.s_norm), history.s_norm, 'b-', 'LineWidth', 1.5);
xlabel('Iteration', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Log Scale', 'FontSize', 12, 'FontWeight', 'bold');
title('Convergence History', 'FontSize', 14, 'FontWeight', 'bold');
legend({'Objective Value', 'Primal Residual', 'Dual Residual'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);  % Adjust axis font size and line width% Display results in the console
fprintf('Original Signal (x0):\n');
disp(x0(1:10));  % Display first 10 values
fprintf('Reconstructed Signal (x):\n');
disp(x(1:10));   % Display first 10 values

运行结果为

2.2. ADMM求解基追踪（Basis Pursuit）问题

2.2.1 数学形式

Basis Pursuit 的数学形式如下：
$$\underset{x}{\min }\Vert x{\Vert }_1\text{subject to }Ax=b,$$
其中：

• $x\in {\mathbf{R}}^n$ 是优化变量；
• $A\in {\mathbf{R}}^{m\times n},b\in {\mathbf{R}}^m$；
• 通常情况下，$m<n$，即这是一个欠定问题（未知数多于方程数）。

2.2.2 将问题转换为 ADMM 形式

为了用 ADMM 求解，需要将原问题拆分成两个部分：

1. 第一部分处理 ${\ell }_1$ 范数最小化；
2. 第二部分处理等式约束 $Ax=b$。

引入一个辅助变量 $z$​，重写为：
$$\begin{array}{rl}& \underset{x}{\min }\Vert z{\Vert }_1,\\ & \text{subject to }x-z=0,Ax=b.\end{array}$$
将这个问题写成 ADMM 框架：
$$\underset{x,z}{\min }f(x)+\Vert z{\Vert }_1\text{subject to }x-z=0,$$
其中：

• $f(x)$ 是约束 $Ax=b$ 的指示函数：
  $$f(x)=\{\begin{array}{ll}0,& \text{if }Ax=b,\\ +\infty ,& \text{otherwise}.\end{array}$$

Augmented Lagrangian（增广拉格朗日函数）：
$$L(x,z,u)=\Vert z{\Vert }_1+\frac{\rho }{2}\Vert x-z+u{\Vert }_2^2,$$
其中 $u$ 是拉格朗日乘子，$\rho >0$ 是增广拉格朗日参数。

2.2.3 ADMM迭代步骤

1. $x$​-update:
   $$x^{k+1}=\arg \underset{x}{\min }f(x)+\frac{\rho }{2}\Vert x-z^k+u^k{\Vert }_2^2.$$

由于 $f(x)$是约束 $Ax=b$​ 的指示函数，优化约束为：
$$\underset{x}{\min }\frac{\rho }{2}\Vert x-z^k+u^k{\Vert }_2^2\text{subject to }Ax=b.$$
这个问题可以写成：
x k + 1 = Π { x ∣ A x = b } ( z k − u k ) , x^{k+1} = \Pi_{\{x | Ax = b\}}(z^k - u^k), xk+1=Π{x∣Ax=b}​(zk−uk),
其中 Π { x ∣ A x = b } \Pi_{\{x | Ax = b\}} Π{x∣Ax=b}​ 表示投影到约束集 { x ∣ A x = b } \{x | Ax = b\} {x∣Ax=b}。

通过计算可得【见附录2】
$$x^{k+1}=\left(I-A^T(A{A}^T{)}^{-1}A\right)(z^k-u^k)+A^T(A{A}^T{)}^{-1}b,$$

2. $z$​-update:
   $$z^{k+1}:=\arg \underset{z}{\min }\Vert z{\Vert }_1+\frac{\rho }{2}\Vert x^{k+1}-z+u^k{\Vert }_2^2.$$

这一步是一个标准的 ${\ell }_1$​-范数最小化问题，解可以通过 软阈值（soft-thresholding） 算子计算：
$$z^{k+1}=S_{1/\rho }(x^{k+1}+u^k),$$
其中 $S_{\kappa }(\cdot )$​是软阈值操作：
$$S_{\kappa }(a)=\text{sign}(a)\cdot \max (|a|-\kappa ,0).$$

3. $u$-update:

$$u^{k+1}=u^k+x^{k+1}-z^{k+1}.$$

2.2.4 Matlab程序和例子

ADMM_bp的函数

function [z, history] = ADMM_bp(A, b, rho, alpha)
% ADMM_bp   Basis Pursuit solver using ADMM
% 
% [x, history] = ADMM_bp(A, b, rho, alpha)
%
% This function solves the Basis Pursuit problem:
%     minimize ||x||_1
%     subject to Ax = b
%
% Basis Pursuit is widely used for finding sparse solutions to underdetermined
% systems of linear equations. This problem is reformulated in ADMM form as:
%     minimize f(x) + ||z||_1
%     subject to x - z = 0
% where:
%     - f(x) is the indicator function of the affine constraint Ax = b:
%          f(x) = 0 if Ax = b, +inf otherwise.
%
% The augmented Lagrangian for this problem is:
%     L(x, z, u) = ||z||_1 + (rho/2)||x - z + u||_2^2
% where u is the scaled dual variable.
%
% ADMM iteratively updates variables as follows:
%  1. x-update: Solve
%         x^{k+1} = argmin_x f(x) + (rho/2)||x - z^k + u^k||_2^2
%         Using projection onto {x | Ax = b}, the explicit solution is:
%         x^{k+1} = (I - A^T(AA^T)^-1A)(z^k - u^k) + A^T(AA^T)^-1b
%
%  2. z-update: Solve
%         z^{k+1} = argmin_z ||z||_1 + (rho/2)||x^{k+1} - z + u^k||_2^2
%         This is a standard l1 minimization problem solved via soft-thresholding:
%         z^{k+1} = S_{1/rho}(x^{k+1} + u^k)
%
%  3. u-update:
%         u^{k+1} = u^k + x^{k+1} - z^{k+1}
%
% Input:
%   A      - m x n coefficient matrix (m < n for underdetermined systems)
%   b      - m x 1 vector of observations
%   rho    - augmented Lagrangian parameter
%   alpha  - over-relaxation parameter (1.0 <= alpha <= 1.8 typical)
%
% Output:
%   x       - n x 1 solution vector
%   history - structure containing:
%               - objective value at each iteration
%               - primal and dual residuals
%               - tolerances for convergence%% Constants and initializationt_start = tic;               % Start timerMAX_ITER = 1000;             % Maximum number of iterationsABSTOL   = 1e-4;             % Absolute tolerance for convergenceRELTOL   = 1e-2;             % Relative tolerance for convergence[~, n] = size(A);            % Dimensions of the problem% Initialize variablesz = zeros(n, 1);             % Auxiliary variable zu = zeros(n, 1);             % Dual variable (scaled Lagrange multiplier)% precompute static variables for x-update (projection on to Ax=b)AAt = A * A';P = eye(n) - A' * (AAt \ A);q = A' * (AAt \ b);% Print table header for iteration logsfprintf('%-5s %-12s %-12s %-12s %-12s %-12s\n',...'Iter', 'r_norm', 'eps_pri', 's_norm', 'eps_dual', 'Objective');for k = 1:MAX_ITER% x-update: Projection onto {x | Ax = b}x = P * (z - u) + q;% z-update: Soft-thresholdingzold = z;x_hat = alpha * x + (1 - alpha) * zold;z = shrinkage(x_hat + u, 1/rho);% u-update: Dual variableu = u + (x_hat - z);% Diagnostics and reportinghistory.objval(k)  = norm(z, 1);  % Objective function valuehistory.r_norm(k)  = norm(x - z);  % Primal residual normhistory.s_norm(k)  = norm(-rho * (z - zold));  % Dual residual norm% Convergence toleranceshistory.eps_pri(k) = sqrt(n) * ABSTOL + RELTOL * max([norm(x), norm(-z)]);history.eps_dual(k) = sqrt(n) * ABSTOL + RELTOL * norm(rho * u);% Print iteration detailsfprintf('%-5d %-12.4f %-12.4f %-12.4f %-12.4f %-12.2f\n', k, ...history.r_norm(k), history.eps_pri(k), ...history.s_norm(k), history.eps_dual(k), history.objval(k));% Check for convergenceif (history.r_norm(k) < history.eps_pri(k) && ...history.s_norm(k) < history.eps_dual(k))break;endend% Print total computation timefprintf('Total time: %.2f seconds\n', toc(t_start));
endfunction y = shrinkage(a, kappa)
% Soft-thresholding operatory = sign(a) .* max(abs(a) - kappa, 0);
end

一个调用的例子

clear; close all; clc;% Example of ADMM for Basis Pursuit (BP)
rng(0);  % Set random seed for reproducibility% Problem dimensions
n = 30;  % Number of variables
m = 30;  % Number of equations% Generate random data
A = randn(m, n);               % Measurement matrix
x = sprandn(n, 1, 0.1 * n);    % Sparse true signal
b = A * x;                     % Observed measurementsxtrue = double(x);  % Save true signal for comparison% Solve Basis Pursuit problem using ADMM
rho = 1.0;   % Augmented Lagrangian parameter
alpha = 1.0; % Over-relaxation parameter
[x, history] = ADMM_bp(A, b, rho, alpha);% Plot results
figure('Position', [100, 100, 800, 600]);% Plot original signal and reconstructed signal
subplot(2, 1, 1);
stem(1:n, xtrue, 'b', 'LineWidth', 1.5); hold on;
stem(1:n, x, 'r--', 'LineWidth', 1.5);
xlabel('Index', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Value', 'FontSize', 12, 'FontWeight', 'bold');
title('Original Signal vs. Reconstructed Signal', 'FontSize', 14, 'FontWeight', 'bold');
legend({'True Signal', 'Reconstructed Signal'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);% Plot convergence history
subplot(2, 1, 2);
semilogy(1:length(history.objval), history.objval, 'k-', 'LineWidth', 1.5); hold on;
semilogy(1:length(history.r_norm), history.r_norm, 'r-', 'LineWidth', 1.5);
semilogy(1:length(history.s_norm), history.s_norm, 'b-', 'LineWidth', 1.5);
xlabel('Iteration', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Log Scale', 'FontSize', 12, 'FontWeight', 'bold');
title('Convergence History', 'FontSize', 14, 'FontWeight', 'bold');
legend({'Objective Value', 'Primal Residual', 'Dual Residual'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);% Display results
fprintf('True Signal with none zero element:\n');
disp(xtrue(1:10));fprintf('Reconstructed Signal (first 10 values):\n');
disp(x(1:10));

运行结果

2.3. ADMM求解Lasso问题

2.3.1 数学形式

Lasso 是一种带有 ${\ell }_1$​ 正则化的线性回归问题，其目标是平衡数据拟合误差和模型的稀疏性。数学形式为：
$$\underset{x}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\lambda \Vert x{\Vert }_1,$$
其中：

• $x\in {\mathbf{R}}^n$；
• $A\in {\mathbf{R}}^{m\times n},b\in {\mathbf{R}}^m$：数据矩阵和观测向量；
• $\lambda >0$​：正则化参数，用于控制稀疏性。

2.3.2 将问题转换为ADMM形式

Lasso 问题可以通过引入辅助变量 $z$​ 转化为：
$$\underset{x,z}{\min }f(x)+g(z),\text{subject to }x-z=0,$$
其中：

• $f(x)=\frac{1}{2}\Vert Ax-b{\Vert }_2^2$：用于拟合数据；
• $g(z)=\lambda \Vert z{\Vert }_1$：用于实现稀疏性约束。

增广拉格朗日函数为：
$$L(x,z,u)=f(x)+g(z)+\frac{\rho }{2}\Vert x-z+u{\Vert }_2^2,$$
其中：

• $u$：是缩放的拉格朗日乘子；
• $\rho >0$：是增广拉格朗日参数。

2.3.3 ADMM 的迭代步骤

1. $x$-update：
   $$x^{k+1}=\arg \underset{x}{\min }f(x)+\frac{\rho }{2}\Vert x-z^k+u^k{\Vert }_2^2.$$

目标函数为：
$$\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\frac{\rho }{2}\Vert x-z^k+u^k{\Vert }_2^2.$$
这是一个标准的带正则化的最小二乘问题（Ridge Regression），其显式解为：
$$x^{k+1}=(A^{T}A+\rho I{)}^{-1}(A^{T}b+\rho (z^k-u^k)),$$

• $A^{T}A+\rho I$ 总是可逆的（$\rho >0$）。

在实现时，可以通过对 $A^{T}A+\rho I$ 进行一次 Cholesky 分解 【见附录1】来加速后续迭代。当$m<n$时，如果仍然去直接分解 $(A^{\mathsf{T}}A+\rho I)\in {\mathbb{R}}^{n\times n}$，其规模会较大（因为 $n$大于 $m$​），可能效率并不理想。此时更常见、更高效的做法是转而分解
$$(A{A}^{\mathsf{T}}+\rho I)\in {\mathbb{R}}^{m\times m},$$
这样做的原因是，当 $m<n$ 时，$(A{A}^{\mathsf{T}}+\rho I)$ 的维度$m\times m$更小，分解规模相应更小，计算会更高效。由$A^{T}A+\rho I$得到($A{A}^{\mathsf{T}}+\rho I$)可以借用著名的 Woodbury 恒等式(亦称 Sherman–Morrison–Woodbury 公式)【见附录3】：
$$(\rho I+A^{\mathsf{T}}A{)}^{-1}=\frac{1}{\rho }I-\frac{1}{\rho }{A}^{\mathsf{T}}(I+\frac{1}{\rho }A{A}^{\mathsf{T}}{)}^{-1}\frac{1}{\rho }A.$$
将它稍作整理，可以写成
$$(\rho I+A^{\mathsf{T}}A{)}^{-1}=\frac{1}{\rho }I-\frac{1}{{\rho }^2}{A}^{\mathsf{T}}(I+\frac{1}{\rho }A{A}^{\mathsf{T}}{)}^{-1}A.$$
因此，当我们要解
$$(A^{\mathsf{T}}A+\rho I)x=q$$
时（$q=A^{T}b+\rho (z^k-u^k)$），可以直接乘上它的逆得到
$$x=(A^{\mathsf{T}}A+\rho I{)}^{-1}q=[\frac{1}{\rho }I-\frac{1}{{\rho }^2}{A}^{\mathsf{T}}(I+\frac{1}{\rho }A{A}^{\mathsf{T}}{)}^{-1}A]q.$$

---

2. $z$-update：
   $$z^{k+1}=\arg \underset{z}{\min }g(z)+\frac{\rho }{2}\Vert x^{k+1}-z+u^k{\Vert }_2^2.$$

这是一个**${\ell }_1$​-范数最小化问题**，可以通过 软阈值算子（Soft-Thresholding Operator） 得到解析解：
$$z^{k+1}=S_{\lambda /\rho }(x^{k+1}+u^k),$$
其中：
$$S_{\kappa }(a)=\text{sign}(a)\cdot \max (|a|-\kappa ,0).$$
软阈值操作会对小于$\lambda /\rho$的系数进行稀疏化。

---

3. $u$-update：
   $$u^{k+1}=u^k+x^{k+1}-z^{k+1}.$$

2.3.4 Matlab程序和例子

ADMM_lasso函数

function [x, history] = ADMM_lasso(A, b, lambda, rho, alpha)
% ADMM_lasso   Solves the Lasso problem using ADMM
%
% This function solves the Lasso problem, which is defined as:
%   minimize (1/2)||Ax - b||_2^2 + lambda * ||x||_1,
%
% where:
%   - A: m x n data matrix
%   - b: m x 1 observation vector
%   - lambda: regularization parameter controlling sparsity
%
% The Lasso problem seeks to balance the data fitting term (||Ax - b||_2^2)
% with the sparsity-inducing regularization term (||x||_1).
%
% In ADMM form, the problem is reformulated as:
%   minimize f(x) + g(z),
%   subject to x - z = 0,
%
% where:
%   - f(x) = (1/2)||Ax - b||_2^2,
%   - g(z) = lambda * ||z||_1.
%
% The augmented Lagrangian is:
%   L(x, z, u) = f(x) + g(z) + (rho/2)||x - z + u||_2^2,
%
% ADMM proceeds with the following updates:
%   1. x-update: Solve a ridge regression problem:
%        x^{k+1} = argmin_x f(x) + (rho/2)||x - z^k + u^k||_2^2.
%      This step involves solving:
%        x^{k+1} = (A^T A + rho * I)^{-1}(A^T b + rho * (z^k - u^k)).
%
%   2. z-update: Apply the soft-thresholding operator:
%        z^{k+1} = S_{lambda/rho}(x^{k+1} + u^k).
%      Soft-thresholding is defined as:
%        S_kappa(a) = sign(a) * max(|a| - kappa, 0).
%
%   3. u-update: Update the scaled dual variable:
%        u^{k+1} = u^k + x^{k+1} - z^{k+1}.
%
% Input:
%   - A: m x n matrix of predictors
%   - b: m x 1 vector of responses
%   - lambda: regularization parameter
%   - rho: augmented Lagrangian parameter
%   - alpha: over-relaxation parameter (1.0 <= alpha <= 1.8 is typical)
%
% Output:
%   - z: solution vector (n x 1)
%   - history: structure containing convergence metrics:
%       * objval: objective value at each iteration
%       * r_norm: primal residual norm
%       * s_norm: dual residual norm
%       * eps_pri: primal feasibility tolerance
%       * eps_dual: dual feasibility tolerance%% Constants and initializationt_start = tic;               % Start timerMAX_ITER = 1000;             % Maximum number of iterationsABSTOL   = 1e-4;             % Absolute tolerance for convergenceRELTOL   = 1e-2;             % Relative tolerance for convergence[m, n] = size(A);            % Dimensions of the problem% Save a matrix-vector multiplyAtb = A' * b;z = zeros(n,1);              % Initialize auxiliary variable zu = zeros(n,1);              % Initialize scaled dual variable u% Cache the factorization[L, U] = factor(A, rho);% Print table header for iteration logsfprintf('%-5s %-12s %-12s %-12s %-12s %-12s\n',...'Iter', 'r_norm', 'eps_pri', 's_norm', 'eps_dual', 'Objective');for k = 1:MAX_ITER% x-update: Solve ridge regression problem% if m >= n, conduct Cholesky factorization directly;% if m < n, use the Sherman–Morrison–Woodbury equation:% (rho I + ATA)^-1 = 1/rho I - 1/rho^2 AT(I + 1/rho AAT)^-1Aq = Atb + rho*(z - u);    % Temporary valueif( m >= n )    % If skinnyx = U \ (L \ q);else            % If fatx = q/rho - (A'*(U \ ( L \ (A*q) )))/rho^2;end% z-update with relaxationzold = z;x_hat = alpha*x + (1 - alpha)*zold;z = shrinkage(x_hat + u, lambda/rho);% u-update: Update scaled dual variableu = u + (x_hat - z);% Diagnostics, reporting, termination checkshistory.objval(k)  = objective(A, b, lambda, x, z);history.r_norm(k)  = norm(x - z);history.s_norm(k)  = norm(-rho*(z - zold));history.eps_pri(k) = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z));history.eps_dual(k)= sqrt(n)*ABSTOL + RELTOL*norm(rho*u);% Print iteration detailsfprintf('%-5d %-12.4f %-12.4f %-12.4f %-12.4f %-12.2f\n', k, ...history.r_norm(k), history.eps_pri(k), ...history.s_norm(k), history.eps_dual(k), history.objval(k));if (history.r_norm(k) < history.eps_pri(k) && ...history.s_norm(k) < history.eps_dual(k))break;endend% Print total computation timefprintf('Total time: %.2f seconds\n', toc(t_start));
endfunction p = objective(A, b, lambda, x, z)% Computes the objective value:%   p = (1/2)||Ax - b||_2^2 + lambda * ||z||_1p = ( 1/2*sum((A*x - b).^2) + lambda*norm(z,1) );
endfunction y = shrinkage(a, kappa)% Soft-thresholding operator:%   y = sign(a) .* max(abs(a) - kappa, 0)y = sign(a) .* max(abs(a) - kappa, 0);
endfunction [L, U] = factor(A, rho)
% FACTOR Precomputes Cholesky factorization for x-update in ADMM
% This function prepares a precomputed matrix factorization to speed up 
% the x-update step in the ADMM iterations. Depending on the size and shape 
% of A, it efficiently handles the matrix to solve the linear system.
% 
% INPUT:
%   A    - m x n matrix (data matrix in the optimization problem)
%   rho  - augmented Lagrangian parameter
% 
% OUTPUT:
%   L    - Lower triangular matrix from Cholesky decomposition
%   U    - Upper triangular matrix (transpose of L)
% 
% BACKGROUND:
% In the ADMM x-update step, we solve a linear system of the form:
%   (A^T A + rho I) x = q (if m >= n, "skinny" case), or
%   x = A^T y with (I + 1/rho A A^T) y = b (if m < n, "fat" case).
% This function computes the Cholesky factorization of the respective 
% matrices for efficient iterative solving.% Get dimensions of the input matrix A[m, n] = size(A);% Determine the matrix to factorize based on the shape of Aif (m >= n)  % "Skinny" case: More rows than columns% Factorize (A^T A + rho * I)% This is a (n x n) symmetric positive definite matrixL = chol(A' * A + rho * speye(n), 'lower');else  % "Fat" case: More columns than rows% Factorize (I + 1/rho * A A^T)% This is a (m x m) symmetric positive definite matrixL = chol(speye(m) + (1/rho) * (A * A'), 'lower');end% Convert the result to sparse format for memory efficiency and speed% This ensures MATLAB recognizes the triangular structureL = sparse(L);U = sparse(L');
end

调用示例

clear; close all; clc;% Example of ADMM for Lasso
rng(0);  % Set random seed for reproducibility%% ============ 1. 生成随机数据 ============
m = 500;           % 样本数(观测维度)
n = 1000;           % 特征数(变量维度)
k = 50;            % 稀疏度(真值 x_true 中非零元素个数)
A = randn(m, n);  % 随机生成设计矩阵% 生成稀疏真值 x_true
x_true = zeros(n, 1);
supp = randperm(n, k);        % 在 n 个位置中随机选 k 个作为非零
x_true(supp) = randn(k, 1)*3; % 放大一点非零量级% 生成观测 b (带一些噪声)
noise = 0.01*randn(m, 1);     % 噪声
b = A * x_true + noise;% 正则化系数(可调)
lambda = 0.1;[x, history] = ADMM_lasso(A, b, lambda, 1.0, 1.0);% Plot original vs. reconstructed signals
figure('Position', [100, 100, 800, 600]);
subplot(2, 1, 1);
stem(1:n, x_true, 'b', 'LineWidth', 1.5); hold on;
stem(1:n, x, 'r--', 'LineWidth', 1.5);
xlabel('Index', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Value', 'FontSize', 12, 'FontWeight', 'bold');
title('Original vs. Reconstructed Signal (First 100 Components)', 'FontSize', 14, 'FontWeight', 'bold');
legend({'Original Signal', 'Reconstructed Signal'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);% Plot convergence history
subplot(2, 1, 2);
semilogy(1:length(history.objval), history.objval, 'k-', 'LineWidth', 1.5); hold on;
semilogy(1:length(history.r_norm), history.r_norm, 'r-', 'LineWidth', 1.5);
semilogy(1:length(history.s_norm), history.s_norm, 'b-', 'LineWidth', 1.5);
xlabel('Iteration', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Log Scale', 'FontSize', 12, 'FontWeight', 'bold');
title('Convergence History', 'FontSize', 14, 'FontWeight', 'bold');
legend({'Objective Value', 'Primal Residual Norm', 'Dual Residual Norm'}, 'FontSize', 12);
grid on;
set(gca, 'FontSize', 12, 'LineWidth', 1.2);

结果

附录1. Cholesky 分解

Cholesky 分解是一种针对对称正定矩阵的分解方法；如果有一个对称正定矩阵 $M$​，我们可以把它分解为
$$M=R^{T}R$$
其中 $R$ 是一个上三角矩阵（或有些实现中会取下三角矩阵）。在实际计算中，它常被用来快速解方程组、提高数值稳定性以及节省运算量。

比如在最小二乘问题中，如果我们要解
$$(A^{T}A)x=A^{T}b,$$
一个直接的做法是先对 $A^{T}A$ 做 Cholesky 分解 $A^{T}A=R^{T}R$，然后我们只需要两次三角回代就能快速得到 $x$​。这是因为解方程组
$$R^{T}Rx=A^{T}b$$
可以分两步走：

1. 先解 $R^{T}y=A^{T}b$ 得到 $y$；
2. 再解 $Rx=y$ 得到 $x$。

这相对于一般的高斯消去法，能显著降低运算量并且更稳定。而且，当我们在迭代算法中要反复求解相似的线性方程组（例如 ADMM 的子问题），我们可以在最开始就对矩阵做 Cholesky 分解，然后重复使用这个分解，从而大幅加速后续的求解过程。

根据这个原理，我们先分解得到$R$, 然后根据$R$计算$x$，即
$$x=R^{-1}(R^T{)}^{-1}{A}^{T}b$$

附录2. 基追踪里$x$更新

对于
$$x^{k+1}:=\arg \underset{x}{\min }\frac{\rho }{2}\Vert x-z^k+u^k{\Vert }_2^2\text{subject to }Ax=b.$$
展开目标函数 $\Vert x-z^k+u^k{\Vert }_2^2$​：
$$\Vert x-z^k+u^k{\Vert }_2^2=(x-z^k+u^k{)}^T(x-z^k+u^k).$$
这可以写成：
$$x^{T}x-2x^T(z^k-u^k)+(z^k-u^k{)}^T(z^k-u^k).$$
因此目标函数变为：
$$\frac{\rho }{2}\left(x^{T}x-2x^T(z^k-u^k)+\text{constant}\right).$$
其中“constant”项与 $x$ 无关，可以忽略，优化问题化为：
$$\underset{x}{\min }\frac{\rho }{2}{x}^{T}x-\rho x^T(z^k-u^k)\text{subject to }Ax=b.$$
为引入约束 $Ax=b$，使用拉格朗日乘子 $\lambda$​，构造拉格朗日函数：
$$\mathcal{L}(x,\lambda )=\frac{\rho }{2}{x}^{T}x-\rho x^T(z^k-u^k)+{\lambda }^T(Ax-b).$$
对 $x$ 和$\lambda$ 求偏导，设置导数为 0，得到以下方程组：

对 $x$​ 的导数：
$$\frac{\partial \mathcal{L}}{\partial x}=\rho x-\rho (z^k-u^k)+A^T\lambda =0.$$
整理得到：
$$x=z^k-u^k-\frac{1}{\rho }{A}^T\lambda .$$
对 $\lambda$​ 的导数：
$$\frac{\partial \mathcal{L}}{\partial \lambda }=Ax-b=0.$$
即：
$$Ax=b.$$

---

联立求解
$$A\left(z^k-u^k-\frac{1}{\rho }{A}^T\lambda \right)=b.$$
展开整理：
$$A{z}^k-A{u}^k-\frac{1}{\rho }A{A}^T\lambda =b.$$

$$A{A}^T\lambda =\rho (A{z}^k-A{u}^k-b).$$

从而得到$\lambda$​ 的解：
$$\lambda =(A{A}^T{)}^{-1}\rho (A{z}^k-A{u}^k-b).$$
将$\lambda$ 的解代回
$$x=z^k-u^k-\frac{1}{\rho }{A}^T\lambda .$$
代入 $\lambda =(A{A}^T{)}^{-1}\rho (A{z}^k-A{u}^k-b)$​：
$$x=z^k-u^k-A^T(A{A}^T{)}^{-1}(A{z}^k-A{u}^k-b).$$
进一步整理为：
$$x^{k+1}=\left(I-A^T(A{A}^T{)}^{-1}A\right)(z^k-u^k)+A^T(A{A}^T{)}^{-1}b.$$

附录3. Woodbury 恒等式（Sherman–Morrison–Woodbury 公式）

**Woodbury 恒等式（Sherman–Morrison–Woodbury 公式）**是一条在数值分析和应用数学中非常著名、且常常被用来对矩阵进行「低秩修正」（low-rank update）求逆的公式。它能够将一个大规模矩阵的逆问题，转化为对一个更小规模矩阵做运算，从而显著节省计算开销。

给定可逆矩阵 $B\in {\mathbb{R}}^{n\times n}$，以及矩阵 $U\in {\mathbb{R}}^{n\times k}$、$C\in {\mathbb{R}}^{k\times k}$、$V\in {\mathbb{R}}^{k\times n}$（其中 $C$​ 也要求可逆），那么 Woodbury 恒等式表述如下：
$$(B+UCV{)}^{-1}=B^{-1}-B^{-1}U(C^{-1}+V{B}^{-1}U{)}^{-1}V{B}^{-1}.$$
当 $U$ 和 $V$ 的列秩（或行秩）很小，且 $C$ 维度小（例如$k\ll n$），就可以用 Woodbury 公式来避开对一个 $n\times n$ 大矩阵的求逆或分解，而改去操作一个$k\times k$ 的矩阵，大大降低了计算量和内存消耗。

对于$(A^{\mathsf{T}}A+\rho I{)}^{-1}$，套用公式有$B=\rho I,U=A^T,C=I,V=A$
$$(A^{\mathsf{T}}A+\rho I{)}^{-1}=\frac{1}{\rho }-\frac{1}{\rho }{A}^T(I+\frac{1}{\rho }A{A}^T{)}^{-1}A\frac{1}{\rho }$$

参考文献

Boyd S, Parikh N, Chu E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends® in Machine learning, 2011, 3(1): 1-122.