ADMM算法

ADMM

对偶下降

$$\begin{gathered} minimizef(x) \\ subjecttoAx=b \end{gathered}$$

其中，$x\in R^n$，$A\in R^{m\times n}$，$f:R^n\to R$

该问题的拉格朗日式如下
$$L(x,y)=f(x)+y^T(Ax-b)$$
对偶方程如下：
$$g(y)=*{inf}_{x}L(x,y)=-*{sup}_{x\in domf}(-y^{T}Ax-f(x))-b^{T}y=-f^{\ast }(-A^{T}y)-b^{T}y$$
为什么这样做呢？
$$\begin{gathered} \underset{x}{\min }L(x,y)=d^{\ast }\le p^{\ast }=\underset{x}{\min }f(x)+I(Ax-b) \\ whereI(u)=\{\begin{array}{ll}0& \text{u=0}\\ \infty & \text{u}̸=0\end{array} \end{gathered}$$
对偶问题最优化
$$maximizeg(y)$$
令$y\ast =*{argmax}_{y}g(y)$，$x\ast =argmi{n}_{x}L(x,y\ast )$

如果强对偶成立，则$d\ast =p\ast$。$x\ast$是原问题的最优解

对偶下降算法步骤如下：

$$\begin{gathered} x^{k+1}:=argmi{n}_{x}L(x,y^k) \\ {y}^{k+1}:=y^k+{\alpha }^k(A{x}^{k+1}-b) \end{gathered}$$

先最小化$L(x,y^k)$，后最大化$L(x^{k+1},y)$

特点：

• ${\alpha }^k$选择合适，且满足一定条件下，$x^k$，$y^k$收敛到最优点
• 然而，对许多例子不成立。比如$f(x)=Bx$，在更新$x$时不存在最小值

对偶分解

每次迭代需要广播和集合两个步骤：$x$更新可以独立进行，$y$更新要集中进行

增广拉格朗日乘子法

$$L_{\rho }(x,y)=f(x)+y^T(Ax-b)+(\rho /2){\Vert Ax-b\Vert }_2^2$$

这等价于如下问题原拉格朗日：
$$\begin{gathered} minimizef(x)+(\rho /2){\Vert Ax-b\Vert }_2^2 \\ subjecttoAx=b \end{gathered}$$

算法步骤：

$$\begin{gathered} x^{k+1}:=argmi{n}_x{L}_{\rho }(x,y^k) \\ {y}^{k+1}:=y^k+\rho (A{x}^{k+1}-b) \end{gathered}$$

为什么选择$\rho$作为迭代步长？

最优条件必须满足
$$\begin{gathered} A{x}^{\ast }-b=0 \\ \nabla f(x^{\ast })+A^T{y}^{\ast }=0 \end{gathered}$$
根据定义，$x^{k+1}$最小化$L_{\rho }(x,y^k)$，因此
$$\begin{gathered} 0={\nabla }_x{L}_{\rho }(x^{k+1},y^k) \\ ={\nabla }_{x}f(x^{k+1})+A^T(y^k+\rho (A{x}^{k+1}-b)) \\ ={\nabla }_{x}f(x^{k+1})+A^T{y}^{k+1} \end{gathered}$$
随着迭代的进行，只需要原条件的残差收敛到0，即可

评价：

• 大大增加了收敛性
• 不能进行分解

ADMM算法

$$\begin{gathered} minimizef(x)+g(z) \\ subjecttoAx+Bz=c \end{gathered}$$

最优解表示如下：
p ∗ = i n f { f ( x ) + g ( z )   ∣   A x + B z = c } p^*=inf\{ f(x)+g(z)\ |\ Ax+Bz=c \} p∗=inf{f(x)+g(z) ∣ Ax+Bz=c}
增广拉格朗日：
$$L_{\rho }(x,y)=f(x)+g(z)+y^T(Ax+Bz-c)+(\rho /2){\Vert Ax+Bz-c\Vert }_2^2$$
算法步骤：
$$\begin{gathered} x^{k+1}:=*{argmin}_x{L}_{\rho }(x,z^k,y^k) \\ {z}^{k+1}:=*{argmin}_z{L}_{\rho }(x^{k+1},z,y^k) \\ {y}^{k+1}:=y^k+\rho (A{x}^{k+1}+B{z}^{k+1}-c) \end{gathered}$$

归一化形式：

令$r=Ax+Bz-c$，则
$$y^{T}r+(\rho /2){\Vert r\Vert }_2^2=(\rho /2){\Vert r+u\Vert }_2^2-(\rho /2){\Vert u\Vert }_2^2$$
其中，$u=(1/\rho )y$是归一化对偶变量

归一化算法步骤如下：
$$\begin{gathered} x^{k+1}=*{argmin}_x(f(x)+(\rho /2){\Vert Ax+B{z}^k-c+u^k\Vert }_2^2) \\ {z}^{k+1}=*{argmin}_z(g(z)+(\rho /2){\Vert A{x}^{k+1}+Bz-c+u^k\Vert }_2^2) \\ {u}^{k+1}=u^k+A{x}^{k+1}+B{z}^{k+1}-c \end{gathered}$$
可见，
$$\begin{gathered} u^k=u^0+\sum _{j=i}^k{r}^j \\ {r}^k=A{x}^k+B{z}^k-c \end{gathered}$$

收敛性

• 假设1： e p i   f = { ( x , t ) ∈ R n × R   ∣   f ( x ) ≤ t } epi \ f=\{(x,t)\in R^n×R \ |\ f(x)\le t \} epi f={(x,t)∈Rn×R ∣ f(x)≤t}是闭的且非空的凸集
• 假设2：存在$(x^{\ast },z^{\ast },y^{\ast })$使得下式成立

$$L_0(x^{\ast },z^{\ast },y)\le L_0(x^{\ast },z^{\ast },y^{\ast })\le L_0(x,z,y^{\ast })$$

满足上述条件下，有

• 残差收敛：$r^k\to 0ask\to \infty$
• 目标函数收敛：$f(x^k)+g(z^k)\to p^{\ast }ask\to \infty$。其中，$p^{\ast }$是最优值
• 对偶变量收敛：$y^k\to y^{\ast }ask\to \infty$。其中，$y^{\ast }$是对偶最优点

注意，$x^k,y^k$不一定收敛到最优值，仅在上述条件满足下。

最优条件

• 原问题的可行性

$$A{x}^{\ast }+B{z}^{\ast }-c=0$$

• 对偶可行性

$$0=\nabla f(x^{\ast })+A^T{y}^{\ast }$$

$$0=\nabla g(z^{\ast })+B^T{y}^{\ast }$$

因为$z^{k+1}$最小化$L_{\rho }(x^{k+1},z,y^k)$，所以有
$$\begin{gathered} 0=\nabla g(z^{k+1})+B^T{y}^k+\rho B^T(A{x}^{k+1}+B{z}^{k+1}-c) \\ =\nabla g(z^{k+1}++B^T{y}^k+\rho B^T{r}^{k+1}) \\ =\nabla g(z^{k+1})+B^T{y}^{k+1} \end{gathered}$$
所以，式(25)总是成立的。

因为$x^{k+1}$最小化$L_{\rho }(x,z^k,y^k)$，所以有
$$\begin{gathered} 0=\nabla f(x^{k+1})+A^T{y}^k+\rho A^T(A{x}^{k+1}+B{z}^k-c) \\ =\nabla f(x^{k+1}++A^T(y^k+\rho r^{k+1}+\rho B(z^k-z^{k+1}))) \\ =\nabla f(x^{k+1})+A^T{y}^{k+1}+\rho A^{T}B(z^k-z^{k+1}) \end{gathered}$$
因此，$s^{k+1}=\rho A^{T}B(z^{k+1}-z^k)$可以看做式(24)的残差

停止条件

收敛条件得出：
$$f(x^k)+g(z^k)-p^{\ast }\le -(y^k{)}^T{r}^k+(x^k-x^{\ast }{)}^T{s}^k$$
假设$||x^k-x^{\ast }||\le d$，则
$$f(x^k)+g(z^k)-p^{\ast }\le -(y^k{)}^T{r}^k+d||s^k|{|}_2\le ||y^k|{|}_2||r^k|{|}_2+d||s^k|{|}_2$$
因此，终止条件设为：
$$||r^k|{|}_2\le {ϵ}^{pri}and||s^k|{|}_2\le {ϵ}^{dual}$$
其中，
ϵ p r i = p ϵ a b s + ϵ r e l m a x { ∣ ∣ A x k ∣ ∣ 2 , ∣ ∣ B z k ∣ ∣ 2 , ∣ ∣ c ∣ ∣ 2 } ϵ d u a l = n ϵ a b s + ϵ r e l ∣ ∣ A T y k ∣ ∣ 2 \epsilon^{pri} = \sqrt{p}\epsilon^{abs}+\epsilon^{rel}max\{ ||Ax^k||_2,||Bz^k||_2 ,||c||_2\} \\ \epsilon^{dual} = \sqrt{n}\epsilon^{abs}+\epsilon^{rel}||A^Ty^k||_2 ϵpri=p ​ϵabs+ϵrelmax{∣∣Axk∣∣2​,∣∣Bzk∣∣2​,∣∣c∣∣2​}ϵdual=n ​ϵabs+ϵrel∣∣ATyk∣∣2​
其中，$p,n$是各自向量的维度

ADMM变种

变化的惩罚参数

$${\rho }^{k+1}:=\{\begin{array}{ll}{\tau }^{incr}{\rho }^k,& \text{if}||r^k|{|}_2>\mu ||s^k|{|}_2\\ {\rho }^k/{\tau }^{decr},& \text{if}||s^k|{|}_2>\mu ||r^k|{|}_2\\ {\rho }^k,& otherwise\end{array}$$

其中，一般$\mu =10>1,{\tau }^{incr}=2>1,{\tau }^{decr}=2>1$

直观地，$\rho$变大，会使得$r^k$变小

广义增广

$$(\rho /2)||r|{|}_2^2\to (1/2)r^{T}Pr$$

其中，P是正定的。这可以看做一个新的标准ADMM。不过，$r=0\to Fr=0$。其中，$F^{T}F=P$

超松弛

在更新z和y时，
$$A{x}^{k+1}\to {\alpha }^{k}A{x}^{k+1}-(1-{\alpha }^k)(B{z}^k-c)$$
一般，选择超松弛参数${\alpha }^k\in [1.5,1.8]$

常见形式

下面只考虑$x$更新步骤，其他类似
$$\begin{gathered} x^{+}=argmi{n}_x(f(x)+\rho /2||Ax-v|{|}_2^2) \\ wherev=-Bz+c-u \end{gathered}$$

二次目标函数

$$f(x)=(1/2)x^{T}Px+q^{T}x+r$$

则
$$x^{+}=(P+\rho A^{T}A{)}^{-1}(\rho A^{T}v-q)$$
式(43)可以利用稀疏性、保存重复的分解等来加快速度

元素软阈值

当$f(x)=\lambda ||x|{|}_1，A=I$时。对于每个标量$x_i$的更新为
$$x_i^{+}:=argmi{n}_{x_i}(\lambda |x_i|+(\rho /2)(x_i-v_i{)}^2)$$
尽管，$f(x)$不可导。但是存在简单的解析解：
$$x_i^{+}:=S_{\lambda /\rho }(v_i)$$
软阈值操作符$S$：
$$S_k(a)=\{\begin{array}{ll}a-k& a>k\\ 0& |a|\le k\\ a+k& a<-k\end{array}$$
或者表示为
$$S_k(a)=(a-k{)}_{+}-(-a-k{)}_{+}$$
这也是shrinkage operator，如下
$$S_k(a)=(1-k/|a|{)}_{+}a,fora̸=0$$

一阶范数问题

最小化$||Ax-b|{|}_1$，转化为ADMM形式
$$\begin{gathered} minimize||z|{|}_1 \\ subjecttoAx-z=b \end{gathered}$$
上面式子中，$f=0,g=||\ast |{|}_1$。假设$A^{T}A$可逆，则ADMM算法步骤如下：
$$\begin{gathered} x^{k+1}=(A^{T}A{)}^{-1}{A}^T(b+z^k-u^k) \\ {z}^{k+1}=S_{1/\rho }(A{x}^{k+1}-b+u^k) \\ {u}^{k+1}=u^k+A{x}^{k+1}-z^{k+1}-b \end{gathered}$$

openMVG代码

eigen线性代数

template <typename Linear_SolverT>
inline void Compute
(
  const Eigen::SparseMatrix<double>& spd_mat,
  Linear_SolverT * linear_solver
)
{
  linear_solver->compute(spd_mat);
}

template <typename Linear_SolverT>
inline void Compute
(
  const Eigen::MatrixXd& spd_mat,
  Linear_SolverT * linear_solver
)
{
  linear_solver->compute(spd_mat.sparseView());
}

迭代参数

struct Options {
    int max_num_iterations = 1000;
    // Rho is the augmented Lagrangian parameter.
    double rho = 1.0;
    // Alpha is the over-relaxation parameter (typically between 1.0 and 1.8).
    double alpha = 1.0;

    double absolute_tolerance = 1e-4;
    double relative_tolerance = 1e-2;
};

私有成员

  Options options_;

  // Matrix A where || Ax - b ||_1 is the problem we are solving.
  MatrixType a_;

  // Cholesky linear solver.
#ifdef EIGEN_MPL2_ONLY
  using Linear_Solver_T = Eigen::SparseLU<Eigen::SparseMatrix<double>>;
#else
  // Since our linear system will be a SPD matrix we can
  // utilize the Cholesky factorization.
  using Linear_Solver_T = Eigen::SimplicialLLT<Eigen::SparseMatrix<double>>;
#endif
  Linear_Solver_T linear_solver_;

  Eigen::VectorXd Shrinkage
  (
    const Eigen::VectorXd& vec, const double kappa
  ) const
  {
    Eigen::ArrayXd zero_vec(vec.size());
    zero_vec.setZero();
    return zero_vec.max( vec.array() - kappa) -
           zero_vec.max(-vec.array() - kappa);
  }

初始化

  L1Solver
  (
    const Options& options,
    const MatrixType& mat
  )
  : options_(options), a_(mat)
  {
    // Analyze the sparsity pattern once. Only the values of the entries will be
    // changed with each iteration.
    const MatrixType spd_mat = a_.transpose() * a_;
    l1_solver_internal::Compute(spd_mat, &linear_solver_);
  }

核心

  bool Solve
  (
    const Eigen::VectorXd& rhs,  // b = rhs
    Eigen::VectorXd* solution
  )
  {
    // Since constructor was called before we check Compute status
    if (linear_solver_.info() != Eigen::Success)
    {
      std::cerr << "Cannot compute the matrix factorization" << std::endl;
      return false;
    }

    Eigen::VectorXd& x = *solution;
    Eigen::VectorXd z(a_.rows()), u(a_.rows());
    z.setZero();
    u.setZero();

    Eigen::VectorXd a_times_x(a_.rows()), z_old(z.size()), ax_hat(a_.rows());
    // Precompute some convergence terms.
    const double rhs_norm = rhs.norm();
    const double primal_abs_tolerance_eps =
      std::sqrt(a_.rows()) * options_.absolute_tolerance;
    const double dual_abs_tolerance_eps =
      std::sqrt(a_.cols()) * options_.absolute_tolerance;

    for (int i = 0; i < options_.max_num_iterations; ++i)
    {
      // Update x.
      x.noalias() = linear_solver_.solve(a_.transpose() * (rhs + z - u));
      a_times_x.noalias() = a_ * x;
      ax_hat.noalias() = options_.alpha * a_times_x;
      ax_hat.noalias() += (1.0 - options_.alpha) * (z + rhs);

      // Update z and set z_old.
      std::swap(z, z_old);
      z.noalias() = Shrinkage(ax_hat - rhs + u, 1.0 / options_.rho);

      // Update u.
      u.noalias() += ax_hat - z - rhs;

      // Compute the convergence terms.
      const double r_norm = (a_times_x - z - rhs).norm();
      const double s_norm =
        (-options_.rho * a_.transpose() * (z - z_old)).norm();
      const double max_norm =
        std::max({a_times_x.norm(), z.norm(), rhs_norm});
      const double primal_eps =
        primal_abs_tolerance_eps + options_.relative_tolerance * max_norm;
      const double dual_eps =
        dual_abs_tolerance_eps +
        options_.relative_tolerance *
          (options_.rho * a_.transpose() * u).norm();

      // Log the result to the screen.
      // std::ostringstream os;
      // os << "Iteration: " << i << "\n"
      //   << "R norm: " << r_norm << "\n"
      //   << "S norm: " << s_norm << "\n"
      //   << "Primal eps: " << primal_eps << "\n"
      //   << "Dual eps: " << dual_eps << std::endl;
      // std::cout << os.str() << std::endl;

      // Determine if the minimizer has converged.
      if (r_norm < primal_eps && s_norm < dual_eps)
      {
        return true;
      }
    }
    return false;
    }