L1 Minimization With CVXOPT

From corrada.com/articles

Recovery of precision errors for a collection of models requires that a \(\ell_1\)-minimization with equality constraints problem be solved. The original problem is to recover the minimum \(\ell_1\) vector x subject to it solving the linear equation system \(Ax=b\). The L1-magic package has an Appendix A in its documentation stating that the problem is equivalent to the following minimization problem

$$\begin{align} \underset{x,u}{\min} \sum_i u_i \,\, \text{subject to} \,\, x_i - u_i & \leq 0, \\ -x_i - u_i & \leq 0, \\ Ax & = b \end{align}$$

A dummy vector u of the same dimensions as the wanted x vector has been introduced to enforce the minimization constraint

\(\min \sum_{i=1}^{m} |x_i|\)

A solution with Mathematica

The solution I currently have working involves Mathematica. The online documentation for Mathematica's LinearProgramming function has the following solution:

L1Minimization[A_,b_] := Module[
  {B, c, Aall, ball, x, lb, AT},
  {m, n} = Dimensions[A];
  AT = Transpose[A];
  B = SparseArray[{{i_,j_} -> 1},{m,n}];
  Aall = Join[Transpose[Join[B,-AT], Transpose[Join[B,AT]]];
  ball = Join[-b,b];
  c = Join[Table[1,{m}], Table[0,{n}]];
  lb = Table[-Infinity, {m+n}];
  x = LinearProgramming[c, Aall, ball, lb];
  x = Drop[x,m]
]

But Mathematica is not without its warts and this solution does not work on my Windows machine running version 6.0.1.0 because LinearProgramming does not behave right with a SparseArray matrix (eventhough the documentation says it can handle objects of this type!). This is the solution I got to work:

L1Minimization[A_, b_] := Module[
  {B, c, Aall, ball, x, lb, AT, m, n, i, ownSparse},
  {m, n} = Dimensions[A];
  AT = Transpose[A];
  B = SparseArray[{{i_, j_} -> 0}, {n, m}];
  Aall = Transpose[Join[B, AT]];
  For[i = 1, i < n + 1, i++,
   ownSparse = Table[If[j == i || j == (n + i), 1, 0], {j, 2*n}];
   AppendTo[Aall, ownSparse];
   ];
  For[i = 1, i < n + 1, i++,
   ownSparse = Table[If[j == i || j == (n + i), 1, 0], {j, 2*n}];
   ownSparse[[n + i]] = -1;
   AppendTo[Aall, ownSparse];
   ]; 
  ball = Join[Table[{b[[i]], 0}, {i, Length[b]}], Table[{0.0, 1}, {2*n}]];
  c = Join[Table[1, {n}], Table[0, {n}]];
  lb = Table[-Infinity, {n + n}];
  x = LinearProgramming[c, Aall, ball, Method -> "InteriorPoint"];
  x = Drop[x, n]]

The solution explained

Let us break down the solution part by part to see how it works.

  {m, n} = Dimensions[A];

The minimization procedure is applied to under-determined linear systems so \(n > m\). Let us a simple example for A to illustrate the transformations that are needed so take \(m=2, \, n=3\). For example (the numbers do not mean anything but are just a guide),

\(A =\begin{pmatrix} 0.1 & 0.2 & 0.3 \\ 0.4 & 0.5& 0.6 \end{pmatrix}\)

  AT = Transpose[A];

This produces the matrix

\(AT =\begin{pmatrix} 0.1 & 0.4 \\ 0.2 & 0.5 \\ 0.3 & 0.6 \end{pmatrix}\)

B = SparseArray[{{i_, j_} -> 0}, {n, m}];
  Aall = Transpose[Join[B, AT]];

The matrix Aall is now of the form

\(Aall =\begin{pmatrix} 0 & 0 & 0 & 0.1 & 0.2 & 0.3 \\ 0 & 0 & 0 & 0.4 & 0.5 & 0.6 \end{pmatrix}\)

  For[i = 1, i < n + 1, i++,
   ownSparse = Table[If[j == i || j == (n + i), 1, 0], {j, 2*n}];
   AppendTo[Aall, ownSparse];
   ];
  For[i = 1, i < n + 1, i++,
   ownSparse = Table[If[j == i || j == (n + i), 1, 0], {j, 2*n}];
   ownSparse[[n + i]] = -1;
   AppendTo[Aall, ownSparse];
   ]; 

The first For loop appends the vector \(\begin{pmatrix}1,0,0,1,0,0\end{pmatrix}\) and its two cyclic permutations to the matrix Aall. The second For loop appends the vector \(\begin{pmatrix}1,0,0,-1,0,0\end{pmatrix}\) and its two cyclic permutations. So the final form of the Aall matrix is

\(Aall =\begin{pmatrix} 0 & 0 & 0 & 0.1 & 0.2 & 0.3 \\ 0 & 0 & 0 & 0.4 & 0.5 & 0.6 \\ 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & -1 & 0 & 0 \\ 0 & 1 & 0 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & 0 & -1 \end{pmatrix}\)

  ball = Join[Table[{b[[i]], 0}, {i, Length[b]}], Table[{0.0, 1}, {2*n}]];

This creates the matrix

\(ball =\begin{pmatrix} b_1 & 0 \\ b_2 & 0 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \end{pmatrix}\)

Mathematica is told by this matrix what kind of constraint to apply. The second entry on each row specifies the type of constraint. The 0 value means that an equality is required. The first two entries correspond to the constraint \(Ax = b\). The 1 value imposes a greater than or equal constraint. So the last six entries correspond to the constraint equations between the original vector x we are trying to recover and the dummy vector u created to guarantee that we are getting the \(\ell_1\)-minimized vector.

  c = Join[Table[1, {n}], Table[0, {n}]];

This enforces the minimization of the sum of the components of dummy u vector.

Putting it all together, the final lines of the code are requiring that the following problem be solved (going down the rows of the linear program defined by the code:

\(\begin{align} \underset{x,u}{\min} \sum_i u_i \,\, \text{subject to} \,\, Ax & = b, \\ u_i + x_i & \geq 0, \\ u_i - x_i & \geq 0 \end{align}\)

This is equivalent to the linear program at the top of the page.

A solution with cvxopt

How would this Mathematica code be implemented in cvxopt?