一篇论文复现的整体思路和复现记录(三，基础实现篇)

工欲善其事，必先利其器，要写出好的MATLAB代码，先从最基础的代码开始做起。

1 ADMM原理

以下内容整理自Standford University的Boyd老师的课件和论文。
ADMM问题的基本形式：
最优化问题形式包括两组可分离自变量和线性等式约束：
$$\begin{array}{rl}\underset{x,z}{\min }& f(\mathbf{x})+g(\text{z})\\ s.t.& \mathbf{Ax}+\mathbf{Bz}=\mathbf{c}\end{array}$$
写出该问题对应的拉格朗日函数式：
$$L_{\rho }(\mathbf{x},\mathbf{z},\lambda )=f(\mathbf{x})+g(\text{z})+{\mathbf{y}}^T(\mathbf{Ax}+\mathbf{Bz}-\mathbf{c})+\frac{\rho }{2}\Vert \mathbf{Ax}+\mathbf{Bz}-\mathbf{c}{\Vert }_2^2$$
按如下步骤，按照Gauss-Seidel方法更新迭代：
$$\begin{array}{rl}{x}^{k+1}& =\underset{x}{\mathrm{argmin}}{L}_{\rho }(\mathbf{x},{\mathbf{z}}^k,{\mathbf{y}}^k)\\ z^{k+1}& =\underset{z}{\mathrm{argmin}}{L}_{\rho }({\mathbf{x}}^{k+1},\mathbf{z},{\mathbf{y}}^k)\\ y^{k+1}& =y^k+\rho ({\mathbf{Ax}}^{k+1}+{\mathbf{Bz}}^{k+1}-c)\end{array}$$
进一步，如在拉格朗日函数中定义${\mathbf{u}}^k=(1/\rho ){\mathbf{y}}^k$，放缩对偶变量，拉格朗日函数变为：
$$L_{\rho }(\mathbf{x},\mathbf{z},\lambda )=f(\mathbf{x})+g(\text{z})+\frac{\rho }{2}\Vert \mathbf{Ax}+\mathbf{Bz}-\mathbf{c}+\mathbf{u}{\Vert }_2^2+const$$
相对应的迭代式子变为：
$$\begin{array}{rl}{\mathbf{x}}^{k+1}& =\underset{x}{\mathrm{argmin}}f(\mathbf{x})+\frac{\rho }{2}\Vert \mathbf{Ax}+{\mathbf{Bz}}^k-\mathbf{c}+{\mathbf{u}}^k{\Vert }_2^2\\ {\mathbf{z}}^{k+1}& =\underset{z}{\mathrm{argmin}}g(\mathbf{z})+\frac{\rho }{2}\Vert {\mathbf{Ax}}^{k+1}+\mathbf{Bz}-\mathbf{c}+{\mathbf{u}}^k{\Vert }_2^2\\ {\mathbf{u}}^{k+1}& ={\mathbf{u}}^k+{\mathbf{Ax}}^{k+1}+{\mathbf{Bz}}^{k+1}-c\end{array}$$
到此为止，上述仍然是一些抽象的概念，与实际问题并无什么联系。好在Boyd老师在网站上面挂出来了一些案例，但是这些案例难懂，我费了一些功夫终于看明白其中一二。

2 部分特殊表示

有些符号看得不是很清楚，总结于以下

邻近算子（$\text{Proximity Operator}$）

前面已经知道$x$的更新值为
$$x^{+}=\underset{x}{\mathrm{argmin}}f(x)+\frac{\rho }{2}\Vert Ax+{Bz}^k-c+u^k{\Vert }_2^2$$
令$v=-Bz+c-u$，
$$x^{+}=\underset{x}{\mathrm{argmin}}f(x)+\frac{\rho }{2}\Vert Ax-v{\Vert }_2^2$$
又令$A=I$：
$$\begin{array}{c}{x}^{+}=\underset{x}{\mathrm{argmin}}f(x)+\frac{\rho }{2}\Vert x-v{\Vert }_2^2\end{array}\text{\hspace{1em}(*)}$$
该式子的右端项可以用${\text{prox}}_{f,\rho }(v)$来表示，我暂时将其翻译为：函数$f$带惩罚项$\rho$的邻近算子。所以，后面有一些式子会用式子(*)来简略表示。

向集合的投影（$\text{Projection}$）

当函数$f$足够简单，$x$的更新值成为前面所说到的邻近算子的形式，可以用解析的办法分析。其中一个例子是，假如说$f$是一个非空闭凸集的指示函数（$\text{indicator function}$），那么也可以将x的更新值表示为：
$$x^{+}=\underset{x}{\mathrm{argmin}}f(x)+\frac{\rho }{2}\Vert x-v{\Vert }_2^2={\Pi }_{\mathcal{C}}(v)$$
其中，${\Pi }_{\mathcal{C}}$表示向$\mathcal{C}$上的欧式范数的投影.

软阈值（$\text{Soft Thresholding}$）

感谢前辈的文章1，2，基本了解软阈值的情况。
软阈值问题的形式为：
$$S_{\kappa }(a)=\{\begin{array}{ll}a-\kappa & a>\kappa \\ 0& |a|\le \kappa \\ a+\kappa & a<-\kappa \end{array}=(a-\kappa {)}_{+}-(-a-\kappa {)}_{+}$$
是在求解优化问题形如：
$$\underset{x}{\mathrm{argmin}}\Vert x-B{\Vert }_2^2+\lambda \Vert x{\Vert }_1$$
的时候作用的，用于标记$x$更新值的取值。其中$B=[b_1,b_2,\dots ,b_n]$。由于$\Vert x{\Vert }_1$并不可微，所以通过分类讨论的结果：
$$S_{\lambda /2}(b)=\{\begin{array}{ll}b-\frac{\lambda }{2}& b>\frac{\lambda }{2}\\ 0& |b|<\frac{\lambda }{2}\\ b+\frac{\lambda }{2}& b<-\frac{\lambda }{2}\end{array}=(b-\frac{\lambda }{2}{)}_{+}-(-b-\frac{\lambda }{2}{)}_{+}$$

例如，在Boyd老师的论文中就提到，$x_i$更新值：
$$x_i^{+}=\underset{x_i}{\mathrm{argmin}}(\lambda |x_i|+(\rho /2)(x_i-v_i{)}^2).$$
可以得到相应更新值：
$$x_i^{+}:=S_{\lambda /\rho }(v_i)$$

3 部分代码解析

以下选择Lasso问题的代码进行分析。代码摘自网站，此处只是复刻编程的流程，整理一下思路。

3.1 代码架构

代码分为lasso.m, objective.m, shrinkage.m, factor.m等，在运行中实际起到作用分别是：

1. function lasso.m是整个ADMM的算法执行流程，包括数据记录，数据迭代，数据输出，迭代终点判断等内容；
2. objective.m是整个优化函数的目标函数；
3. shrinkage.m表示整个优化函数的目标函数；
4. factor.m则是根据A的形状的不同，进行的分解。
5. 通过实际案例检验所写的代码是否正确反映了算法。

3.2 代码解构

3.2.1 lasso.m

lasso问题的目标函数为：
$$minimize\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\lambda \Vert x{\Vert }_1$$
写成ADMM方法所能够求解的格式：
$$\begin{array}{rl}minimize& \frac{1}{2}\Vert Ax-b{\Vert }_2^2+\lambda \Vert z{\Vert }_1\\ \text{subject to}& x-z=0\end{array}$$
迭代的表达式为：
$$\begin{array}{rl}{x}^{k+1}& :=(A^{T}A+\rho I{)}^{-1}(A^{T}b+\rho (z^k-u^k))\\ z^{k+1}& :=S_{\lambda /\rho }(x^{k+1}+u^k)\\ u^{k+1}& :=u^k+x^{k+1}-z^{k+1}\end{array}$$

function [z, history] = lasso(A, b, lambda, rho, alpha)
% lasso  Solve lasso problem via ADMM
% [z, history] = lasso(A, b, lambda, rho, alpha);
% Solves the following problem via ADMM:
%   minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1
% The solution is returned in the vector x.
% history is a structure that contains the objective value, the primal and
% dual residual norms, and the tolerances for the primal and dual residual
% norms at each iteration.
% rho is the augmented Lagrangian parameter.
% alpha is the over-relaxation parameter (typical values for alpha are
% between 1.0 and 1.8).
% More information can be found in the paper linked at
%:http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
%

t_start = tic;
Global constants and defaults
QUIET    = 0;
MAX_ITER = 1000;
ABSTOL   = 1e-4; 
RELTOL   = 1e-2;
Data preprocessing
[m, n] = size(A);

% save a matrix-vector multiply
Atb = A'*b;
ADMM solver
x = zeros(n,1);
z = zeros(n,1);
u = zeros(n,1);

% cache the factorization
[L U] = factor(A, rho);

if ~QUIET
    fprintf('%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
      'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
end

for k = 1:MAX_ITER

    % x-update
    q = Atb + rho*(z - u);    % temporary value
    if( m >= n )    % if skinny
       x = U \ (L \ q);
    else            % if fat
       x = q/rho - (A'*(U \ ( L \ (A*q) )))/rho^2;
    end

    % z-update with relaxation
    zold = z;
    x_hat = alpha*x + (1 - alpha)*zold;
    z = shrinkage(x_hat + u, lambda/rho);

    % u-update
    u = u + (x_hat - z);

    % diagnostics, reporting, termination checks
    history.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);

    if ~QUIET
        fprintf('%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
            history.r_norm(k), history.eps_pri(k), ...
            history.s_norm(k), history.eps_dual(k), history.objval(k));
    end

    if (history.r_norm(k) < history.eps_pri(k) && ...
       history.s_norm(k) < history.eps_dual(k))
         break;
    end

end

if ~QUIET
    toc(t_start);
end
end

3.2.2 objective.m

function p = objective(A, b, lambda, x, z)
    p = ( 1/2*sum((A*x - b).^2) + lambda*norm(z,1) );
end

3.2.3 shrinkage.m

此处表达的是式子：
$$z^{k+1}:=S_{\lambda /\rho }(x^{k+1}+u^k)=(x^{k+1}+u^k-\lambda /\rho {)}_{+}-(-x^{k+1}-u^k-\lambda /\rho )$$

function z = shrinkage(x, kappa)
    z = max( 0, x - kappa ) - max( 0, -x - kappa );
end

3.2.4 factor.m

$A$是一个$m\times n$的矩阵。当$m<n$，由于在$x$的更新式中，$A^{T}A+\rho I$是一个$n\times n$的矩阵，而提一个$1/\rho$后，变为$I+(1/\rho )A{A}^T$，是一个$m\times m$大小的矩阵。由于矩阵求逆的复杂度为$O(n^3)$，因此，后面这种做法更便于求解。结合稀疏矢量技术完成文中代码。

function [L U] = factor(A, rho)
    [m, n] = size(A);
    if ( m >= n )    % if skinny
       L = chol( A'*A + rho*speye(n), 'lower' );
    else            % if fat
       L = chol( speye(m) + 1/rho*(A*A'), 'lower' );
    end

    % force matlab to recognize the upper / lower triangular structure
    L = sparse(L);
    U = sparse(L');
end

3.2.5 难点和实施

难点主要在于要矩阵表达式的运算和推导，容易出错，因此给复现带来了困难。实际案例（Example）见于链接，为上面函数的应用。