8
$\begingroup$

What numerical methods are suitable to solve the following problem

$$\min_x \tfrac{1}{2}x^T A x + b^Tx + \lambda ||x||_1$$

where $x,b\in\mathbf{R}^n$, and $A\in \mathbf{R}^{n\times n}$ is positive definite, and $\lambda\in\mathbf{R}$ is positive?

Typically this problem arises in L1-regularized regression problems (e.g. LASSO) but the common method to solve that (LARS) requires that you can compute regression residuals, which you can't in this case.

Are there standard tools (open source or commercial) that can solve problems like this?

$\endgroup$
  • $\begingroup$ Boyd has a recent monograph called Proximal Algorithms that explains useful methods for problems like this. $\endgroup$ – littleO Jul 17 '14 at 3:58
  • 2
    $\begingroup$ To get something working quickly, you could very easily (a few lines of code in Matlab) implement the proximal gradient method based on the second paragraph of Michael Grant's answer. More information about the proximal gradient method can be found in Vandenberghe's 236c notes here. Then FISTA is only a slight modification of the proximal gradient method. Evaluating the prox operator of the L1 norm is discussed here. $\endgroup$ – littleO Jul 17 '14 at 12:29
  • 1
    $\begingroup$ @littleO, thanks for adding these. If you put them in an answer I would vote them up. One thing I would add---indeed, I think I will edit my answer to do so---is that many of the theoretical discussions use a fixed step size based on a Lipschitz continuity measure; this includes Vandenberghe's notes. That is extremely conservative in practice. Step size adaptation is essential for good performance. $\endgroup$ – Michael Grant Jul 17 '14 at 14:31
  • 1
    $\begingroup$ I'll add one more buzz word: look up "split bregman", or equivalently alternating direction method of multipliers. Works wonders. $\endgroup$ – icurays1 Jul 17 '14 at 15:03
  • $\begingroup$ I would argue that split Bregman and ADMM are overkill for problems of this specific type, due to the strong convexity of the smooth term and the simple structure of the proximal term. But I don't feel strongly about that. And you helped me remember that my co-author has a rather impressive list of first-order methods on his old CalTech page, part of his Sparse and Low Rank Approximation Wiki. $\endgroup$ – Michael Grant Jul 17 '14 at 15:41
7
$\begingroup$

Oh, absolutely, there are lots of tools. Do a search for proximal gradient methods in your favorite search engine; these are generalizations of projected gradient methods. For specific tools, search for TFOCS (disclosure below), FISTA, NESTA, SPGL1, GPSR, SpaRSA, L1LS... and the bibliographies for these will lead to even more options. Even better, see this rather extensive list of first-order methods compiled by Stephen Becker as part of his larger sparse- and low-rank approximation wiki.

A basic proximal gradient method considers your function as the sum of a smooth term $f(x)=\tfrac{1}{2}x^TAx+b^Tx$ and a non-smooth but "prox-capable" term $g(x)=\lambda\|x\|_1$. At each iteration, it uses the gradient of the smooth term $v=\nabla f(x)=Ax+b$ and computes the next step using a proximal minimization $$x^{+} = \mathop{\textrm{arg min}}_z ~ \lambda \|z\|_1 + \tfrac{1}{2t} \|x-tv-z\|_2^2$$ where $t$ is a step size. Many presentations of these methods assume a fixed step size determined by a Lipschitz continuity measure, but such step sizes are always very conservative. So a simple step size adaptation strategy is essential for good performance.

More advanced methods use the same gradient and proximal minimization computations, but mix that information together in different ways to achieve a provable improvement in worst-case performance. But it turns out that for your function is strongly convex, due to your claim that $A$ is positive definite. So a standard proximal gradient algorithm structure is likely a very good choice.

I co-wrote TFOCS with Stephen Becker and Emmanuel Candes. TFOCS is a MATLAB-based toolbox that lets you build custom first-order solvers using a variety of algorithms, smooth functions, linear operators, and prox functions. We wrote an accompanying journal article for it. Our goal was to structure the code in such a way that it is easy to try a variety of common first-order methods on your models, and to make the code that implements the iterations themselves easy to read. But really, it's just one of many software packages out there, and of course if you don't want to use MATLAB, you'll have to try something else!

If you're willing and able to to factorize $A=R^TR$, then you can indeed use standard LASSO tools: $$\min_x \tfrac{1}{2} \| R x - \bar{b} \|_2^2 + \lambda \|x\|_1 - \tfrac{1}{2} \bar{b}^T\bar{b}, \qquad \bar{b}\triangleq - R^{-T} b$$ Don't limit yourself to LARS; there are a lot of alternative approaches out there, and many of the packages above specifically support the LASSO. But the nice thing about many (but not all) of the tools above, including TFOCS, is that you do not necessarily need to do the factorization if you prefer not to.

Finally, this is of course a convex optimization problem, so any general-purpose convex optimization framework will be able to handle this with aplomb, albeit less efficiently than these more specialized methods. Typically they will split $x$ into positive and negative parts: $$x=x^{+} - x^{-}, \quad x^{+}, x^{-} \succeq 0$$ This enables the substitution $\|x\|_1=\vec{1}^T(x^{+}+x^{-})$, in which case the problem becomes a standard QP.

Major edit history:

  • Clarified that the OP's function is strongly convex due to the positive definiteness of $A$.
  • Added a note about the importance of step size adaptation in practice; thanks littleO.
  • Added links to Stephen Becker's sparse approximation wiki; thanks icurays1.
$\endgroup$
  • $\begingroup$ Thanks Michael, this is really helpful - I'm exploring several different options you suggested. I never have to solve a particularly tricky problem (mostly quadratic programming with an $L_1$ term and sometimes linear inequality constraints) but I need to solve a lot of them, so speed is a factor for me. Thanks for your help! $\endgroup$ – Chris Taylor Jul 18 '14 at 17:35
  • $\begingroup$ You're welcome. But do take AJ's suggestion seriously! At least give the general-purpose solvers a try first. Who knows, maybe you'll find it easier to just farm all of your problems out to a cluster than building a highly optimized solver :-) $\endgroup$ – Michael Grant Jul 18 '14 at 17:46
  • $\begingroup$ I'll certainly give it a shot. I've already tried several general-purpose solvers (open source and commercial) which have been okay, but not as fast as a couple of hand-rolled special-purpose solvers. I've not tried CVX yet though... $\endgroup$ – Chris Taylor Jul 18 '14 at 21:38
  • $\begingroup$ Ah, well if you're already that far along in your work that you're looking for speed, then I would not try something like CVX, which is certainly not designed for speed. $\endgroup$ – Michael Grant Jul 19 '14 at 3:18
  • $\begingroup$ @MichaelGrant, Could you explain the intuition behind the Accelerated Methods (Such as FISTA)? It seems they add knowledge from the previous iteration and gain speed because of that. Is there a way to add information for $ l $ iteration backwards and add even more speed? $\endgroup$ – Royi May 31 '16 at 17:53
3
$\begingroup$

CVX (MATLAB) or CVXPY (Python) would allow you to solve that really easily, as long as your problem isn't too large.

The CVX code would be just:

cvx_begin
    variable x(n)
    minimize( 1/2*quad_form(x,A) + b'*x + lambda*norm(x,1) )
cvx_end
$\endgroup$
  • $\begingroup$ Indeed, this really should be the first thing anyone does, before even considering a more specialized algorithm. There's no sense in worrying about complexity or speed until you have reason to believe you'll be satisfied with the results. $\endgroup$ – Michael Grant Jul 18 '14 at 15:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.