Why is SQP better than Augmented Lagrangian for nonlinear programming?

In the technical report on Galahad[1], the authors state, in the context of general nonlinear programming problems,

To our minds, there had never really been much doubt that SQP [sequential quadratic programming] methods would be more successful [than Augmented Lagrangian methods] in the long term...

What could be the basis for that belief? I.e., are there any theoretical results that suggest SQP methods should be faster/more reliable than Augmented Lagrangian methods?

[1] Galahad, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization, by Gould, Orban, and Toint

SQP methods require that the objective is twice differentiable (cf https://en.m.wikipedia.org/wiki/Sequential_quadratic_programming) while Augmented Lagrangians work even when the objective is nondifferentiable (hence their recent resurgence in the image processing community cf ftp://arachne.math.ucla.edu/pub/camreport/cam09-05.pdf)

I don't know about the galahad software, but if is supposed to solve differentiable optimization problems it will probably do much better by using a method that is allowed to differentiate the objective function.

• It isn't true that SQP requires twice differentiable objective functions. You may simply get a method that has a smaller rate of convergence if the objective function has less differentiability, but that is exactly the same as with augmented Lagrangian methods. – Wolfgang Bangerth May 28 '14 at 0:23

In terms of outer iterations, SQP should win because it includes second derivative information, whereas augmented lagrangian methods such as ADMM do not.

However, one thing to keep in mind is that each iteration for these methods involves solving a linear system, so to do a fair comparison you have to take into account how easy these systems are to solve.

For augmented lagrangian (alternating) methods, each iteration you are solving something like, $$(A^TA + \rho I)x = b,$$ where $A$ is a forward operator straight from the objective function that is known and usually easier to deal with or precondition, and $\rho$ is the penalty parameter. (eg, your problem is $\min_x ||Ax-b||^2$ subject to some regularization and constraints).

For SQP methods you are solving something like $$Hx = g,$$ where $H$ is the Hessian (or approximation thereof) which is usually only available implicitly in terms of it's action on vectors, and $g$ is the gradient. The Hessian contains not just $A$, but also a combination of other matrices and matrix inverses coming from linearizing the constraints and regularization.

Preconditioning Hessians is a pretty tricky business and is much less studied than preconditioning forward problems. A standard method is to approximate the Hessian inverse with L-BFGS, but this is of limited effectiveness when the Hessian inverse is high-rank. Another popular method is to approximate the Hessian as a sum of a low-rank matrix plus an easy to invert matrix, but this also has limited effectiveness for hard problems. Other popular Hessian estimation techniques are based on sparse approximations, but continuum problems often have Hessians that have poor sparse approximations.

• +1, although I'd like to caution against blanket statements (by which I don't mean specifically this answer). For example, in PDE-constrained optimization, applying $A$ often involves solving a nonlinear PDE, while $H$ can be applied by solving two linear PDEs -- which can be significantly cheaper (and easier to precondition) if the original PDE is nasty. – Christian Clason May 27 '14 at 22:54
• So, $H$ can be applied by solving two PDEs, but to apply $H^{-1}$ you need to solve 2 PDEs per kryolv iteration in your solver. On the other hand $A$ is a forward operator so it usually doesn't involve any PDE solves at all. Typically one actually knows the matrix $A$ explicitly, eg., a 5 point finite difference stencil on a mesh. Preconditioners for $A$ can be used to build preconditioners for $A^TA + \rho I$, but it is harder to use them to precondition $H$. – Nick Alger May 27 '14 at 23:20
• If $A$ is a linear forward operator (which is not the case in nonlinear PDE-constrained optimization), then you are of course correct. Otherwise, applying $A$ requires a linear PDE solve per Newton iteration (or fixed point iteration), followed by another for $A^T$ (which is always linear). Which of the two methods requires fewer total work (say, by number of linear PDE solves) depends very much on the specific problem. Different tools for different jobs, is all I'm saying. – Christian Clason May 28 '14 at 6:45
• I agree about different tools for different jobs. The Gauss-Newton Hessian for the PDE constrained optimization problem I have in mind - $\min_{q,u} \frac{1}{2}||Cu - y||^2 + \frac{\alpha}{2}||Rq||^2$ such that $Au=q$ - is $H = A^{-T}C^TCA^{-1} + \alpha R^T R$, and the full Hessian is this plus other terms. So here $H$ contains two inverses and $H^{-1}$ contains two inverses within an inverse. – Nick Alger May 28 '14 at 14:52
• And I had the the constraint $S(q) = u$ in mind (e.g., $S$ maps $q$ to the solution $u$ of $-\nabla\cdot(q\nabla u) = f$, which appears in parameter identification or topology optimization). – Christian Clason May 28 '14 at 15:52