What is the fastest way to solve Ax=b (subject to constraints and an absolute term)

Question

I am trying to solve/optimize $Ax=b$ in the least squares sense subject to

1. box constraints;
2. a few (less than 5) equality/inequality constraints; and
3. an absolute function penalty (or some other piecewise linear convex function).

That is, min $f(\mathbf{x})=|A\mathbf{x}-b|^2+\mathbf{c}\cdot|\mathbf{x}|$ where $\mathbf{c}\cdot|\mathbf{x}|=c_1|x_1|+...+c_n|x_n|$

In particular, $A\mathbf{x}=b$ is massively underdetermined (i.e. $A$ has several thousand columns but hundreds of rows) and any valid solution will do if many exist. This code is performance critical. I've exhausted everything I can think of. Is there any obvious way to solve this efficiently (either a clever formulation or appropriate choice of algorithm)?

Language agnostic solutions would be appreciated.

---

tl;dr - Skip the rest. Here's what I have tried so far

Approach 1 - Keep it simple

Use Bounded BFGS and evaluate $f(\mathbf{x})$ directly. I am trying to improve on this as it is too slow.

Approach 2 - QP:

$f(\mathbf{x})=\mathbf{x^T}A^TA\mathbf{x}-2\mathbf{x^T}A^Tb+\mathbf{c}\cdot|\mathbf{x}|=\frac{1}{2}\mathbf{x^T}Q\mathbf{x}+\mathbf{d}\cdot\mathbf{x}+\mathbf{c}\cdot|\mathbf{x}|$

So this is basically standard quadratic form but $Q$ is now enormous and positive semi-definite making it difficult to work with. Slightly slower.

Approach 2 - Use QR decomposition:

The idea is to reduce the size of the quadratic matrix using QR decomposition and a co-ordinate transformation.

$A^T=QR$ and $\mathbf{x^T}Q=\mathbf{y^T}$

$f(\mathbf{x})=\mathbf{x^T}A^TA\mathbf{x}-2\mathbf{x^T}A^Tb+\mathbf{c}\cdot|\mathbf{x}|=\mathbf{y^T}RR^T\mathbf{y}+\mathbf{d}\cdot Q\mathbf{y}+\mathbf{c}\cdot|Q\mathbf{y}|$

This makes $A^TA$ very small and fast to execute. However, the transformation forces the thousands of (quick) box constraints to become linear inequality constraints and this hammers the performance. A lot slower.

Any other ideas? Thanks.

Answer 1

This problem can be formulated as a standard positive semidefinite QP with bounded variables.

First, deal with the absolute value terms in the objective by letting $x=u-v$, where $u\geq 0$ and $v \geq 0$. Then write the problem as

$\min \;\; (A(u-v)-b)^{T}(A(u-v)-b)+c^{T}u+c^{T}v$

subject to

$B(u-v)=d$ (equality constraints)

Upper bounds on u and v.

$u,v \geq 0$.

Then use an active set method for QP to solve this QP. In this approach, most of the variables are fixed at their lower or upper bounds, which dramatically reduces the size of the QP. You solve a simple linear programming feasibility problem with the constraints $B(u-v)=d$ and the bounds constraints to find an initial feasible solution and an initial active set of variables.

In each iteration of the active set method you solve the reduced size QP over the current set of active variables, and then check optimality conditions to see if any of the fixed variables should be released from their bounds and whether any of the free variables should be pinned to their upper or lower bounds.

Active set methods for QP are widely available- don't reinvent the wheel.