Back To ADMM Home.

Table of Contents

1. Huber Fitting
2. Basis pursuit problem
3. L1 Regulaized Loss
4. Lasso
5. Sparse Inverse Covariance Selection
6. Total Variation Reconstruction

6.1 Huber Fitting

The huber fitting problem is to fit the data with huber loss (which is a loss more robust to outliers compared to l2 norm). But as it get quadratic when getting closer to zero, it won’t ouput a sparse result.

Equivalent to the following problem:

The corresponding ADMM updates are :

Where the update of x is solving a least square problem. We will solve it using factorization:

  [L U] = factor(A);
  q = Atb + A'*(z - u);
  x = U \ (L \ q);

the update of z is a little more complicate. But we can solve it exactly by considering the case $|z|<1$ and $|z|>1$, lead to the result:

With matlab code:

  tmp = A*x - b + u;
  z = rho/(1 + rho)*tmp + 1/(1 + rho)*shrinkage(tmp, 1 + 1/rho);

Result:

Code and Test Script.

6.2 Basis pursuit problem

The problem is a l1 minizmization in a linear constrained set :

With m < n, the problem is to find a sparse solution (sparsity corresponding to l1 norm), with underdetermined system of linear equations.

Reform the problem with two functions : \(f(x) = \|x\|_{1}\) and \(g(z) = I_{C}(z)\). g is the indicate function of the set C, where \(C = \{ z \| Az = b \}\). And adding the consensus constraint $x= z$, we will have the following equivalent problem:

Using the scaled ADMM udpates we have:

From the evaluation of proximal operators, we know the proximal operator of indicator function is equivalent to the projection operator. And the proximal oprator of the l1 norm is a shrinkage funtion. As a result, we have the following updates:

6.2.2 Updates

The shrinkage function is:

  function y = shrinkage(a, kappa)
      y = max(0, a-kappa) - max(0, -a-kappa);
  end

The projection into C is :

This solution could be found by solving the following problem (finding the projection of z into the set C):

For the lagrangian:

The gradient with respect to x is zero:

So the dual function is :

Subsititute into the equation $x^{*} = z - A^{T}\lambda/2$, we have:

6.2.3 x Relaxation

However in the test, I found this relaxation did show positive effect.

6.2.4 Result

Code and Test Script.

6.3 L1 Regulaized Loss

Having a l1 loss term for the variable, will produce sparity in the variable pattern. The general problem is:

Transform into consensus form:

The ADMM updates are :

6.3.1 Logistic function

Consider the linear system $Aw+bv=z$ and $x = [v,w]$

The x update is a unconstrainted optimization problem:

Which can be solved using Netwon’s method:

Where:

The Newton’s step will be:

We can also use LBFGs for solving it (when the amount of variables is huge), Code.

6.3.4 Result

Code and Test Script.

6.4 Lasso

See in Proximal Algorithms pages.

Using LSQR we can handle huge problems. Or in Code and Code using LSQR. see their results here and here with LSQR.

The least squares (LSQR) algorithm is an adaptation of the conjugate gradients (CG) method for rectangular matrices. Analytically, LSQR for Ax = b produces the same residuals as CG for the normal equations A’Ax = A’b, but LSQR possesses more favorable numeric properties and is thus generally more reliable. LSQR.

6.5 Sparse Inverse Covariance Selection

Given a dataset consisting of samples from a zero mean Gaussian distribution in $\mathbf{R}^{n}$ :

Consider the task of estimating the covariance matrix under the prior assumption that $\Sigma^{-1}$ is sparse. Since $(\Sigma^{-1})_{i,j} = 0$ means that the ith and jth components of the random variable are conditionally independent, this problem is equivalent to the structure learning problem of estimating the topology of the undirected graphical model representation of the Gaussian. It is also called the covariance selection problem.

For an example, in SLAM problem if we want to determine the corresponding relationship of the feature observations (determine whether it is a match inlier or a outlier). We can apply the covariance selection problem upon the system’s covariance.

If S is the empirical covariance matrix $(1/N)\sum_{i=1}^{N}a_{i}a_{i}^{T}$, then the estimation problem can be written as :

With variable $X \in S_{+}^{n}$, where $|\cdot|_{1}$ is defined elementwise.

6.5.1 ADMM

The upper problem is equivalent to the following problem:

The ADMM algorithm for sparse inverse covariance selection is :

6.5.2 Z update

All the calculation in Z update are done in element-wise, as a result the update of Z is a elementwise soft thresholding:

6.5.3 X update

X update could be further simplify using the following two properties:

Proof of the first derivative:

Using the upper gradient result, and using the first-order optimality condition is that the gradient should vanish:

As $X \succ 0$, we can rewrite the equation and take the orthogonal eigen-value decomposition of the righthand side.

Then we will have :

Now we can construct a diagonal solution of this equation (as we only need a solution for ADMM update) :

Which are always positive since $\rho > 0$. It shows that the computational effort of the X-update is that of an eigenvalue decomposition of a symmetric matrix. Matlab code:

  % x-update
  [Q,L] = eig(rho*(Z - U) - S);
  es = diag(L);
  xi = (es + sqrt(es.^2 + 4*rho))./(2*rho);
  X = Q*diag(xi)*Q';

6.6 Total Variation Reconstruction

Consider a total variation reconstruction of a 1d data(vector date).

Construct a D matrix for the calculation of derivatives, we have the problem:

The upper norm $\mid \mid_{F}$ means elementwised l1 sum. And apply the consensus form (ADMM):

Apply the ADMM updates :

We can easily get the update of x using the first order optimal condition, result in the matlab code :

  x = (I + rho*DtD) \ (b + rho*D'*(z-u));

And the update of z :

  z = shrinkage(D*x + u, lambda/rho);

We will get the convergence result.

We further compare the ADMM method with the CVX SDP method.

  cvx_begin quiet
      variable xtv(n,1);
      xtv_diff = xtv(2:end,1) - xtv(1:end-1,1); % x (horiz) differences
      minimize(norm(xtv_diff, 1) + gamma_l1*norm(xtv-b,2)); %tv roughness measure
  cvx_end

We saw that the CVX SDP method result highly depends on the parameter gamma_l1, while the ADMM results are much more stable and fast.

The running time results are :

Back To ADMM Home.