rank-deficient NNLS

Question

I want to find the minimum-norm solution to a rank-deficient least-squares problem, subject to positivity constraints, e.g. $$\min_x\ \|x\|^2 \quad s.t.\quad Ax = b,\ x \geq 0$$ where $A$ is large, sparse, rectangular, and rank deficient.

What is the most efficient way of robustly solving this problem? Is there recommended code out there for doing so? If $A$ were full-rank, I could use NNLS; without the inequality constraints, I could use a sparse QR decomposition. The above optimization problem is obviously a QP, so worst-case I assume an interior point method should work well, although in the past I have had trouble getting e.g. IpOpt to work with rank-deficient equality constraint matrices.

EDIT: A typical $A$ is of size 50,000 $\times$ 70,000 with 1,000,000 nonzeroes. I would prefer C code, but FORTRAN, python, etc. is fine too. For performance reasons I would like to avoid Matlab if possible.

Answer 1

We solve almost the same problem using the Split-Bregman algorithm (see this for implementation details).

Essentially you convert the problem into $\min_x \|x\|^2+\lambda \|Ax-b^*\|+\gamma\| u -x \|$ where $u$ is $x$ with all negative values removed, and $b^*$ is $b$ with the residual added. This is then solved repeatedly, with $u$ and $b^*$ changing every iteration.

I can write up the entire algorithm for you if you'd like, but it may be easier to simply look at the pseudo-code in the references.

Answer 2

Others have already supplied the two most likely answers to this question, but I'll add a bit of comparison and a way to help decide between the two approaches. I'd suggest either

1. A primal-dual interior point method for quadratic programming. For this approach you'd want to make use of a library that someone else has developed. A lot depends on whether you want/need open source software or whether something that is free for research (but not open source) or commercial and not free for research would do. Some codes to consider include CPLEX, Gurobi, LOQO, and Clp.

or

1. A first order method (Augmented Lagrangian, ADMM, split Bregman, etc.) These are possible to implement yourself without needing to use a packaged library.

The big advantage of the first order methods is that they can solve very large problems that are too big to solve by the primal-dual interior point method but these methods don't produce very accurate solutions and can converge very slowly.

Your problem is of a size that should be reasonably possible to solve using a primal dual interior point method (assuming that you've got a computer with sufficient memory and CPU.) If so, you'll probably be happier with the quality of solutions produced by the interior point method.

You could also consider a more general purpose LP/SOCP/SDP solver such as Mosek, SeDuMi, or SDPT3, but since your problem actually is a QP it doesn't make a lot of sense to reformulate your problem as an SOCP.

Here is a quick test that may help you in determining whether a primal-dual interior point method would be usable for this problem. Compute $H=AA^{T}+I$ (using e.g. MATLAB) then compute the Cholesky factorization of $H$ (using a suitable sparse Cholesky routine- MATLAB has this built in.) If the Cholesky factorization can reasonably be computed on your machine, then a primal-dual interior point code should have enough memory to solve the QP. If $H$ is just way too big and dense for this to be practical, then using a primal-dual method is out of the question and you need to look at first order methods.

Answer 3

You might be able to form a (sparse) $QR$ decomposition of $A$, and then use the $R$ matrix and the system $Rx = Q^{T}b$ to form an equivalent full rank linear system.

Independent of that, you could try other QP solvers, not just IPOPT. CPLEX and Gurobi have free academic licenses and contain QP solvers. Something like CVX, or SDP solvers like SeDuMi and SDPT3 would work as well.

Split augmented Lagrangian approaches like LKlevin is proposing have worked really well for specific problems (Boyd has a good paper for $\ell_{1}$ regression using the alternating direction method of multipliers, for instance), and they can be faster than interior-point solvers by taking advantage of problem structure.