SNOPT — Sparse Sequential Quadratic Programming

1 · Motivation — the problem SNOPT exists to solve

SNOPT ("Sparse Nonlinear OPTimizer") solves the single most common shape of decision problem in engineering and science:

\[ \min_{x}\ f(x) \qquad \text{subject to} \qquad l \le c(x) \le u \]

Here \(x\in\mathbb{R}^n\) are the decision variables we choose, \(f:\mathbb{R}^n\to\mathbb{R}\) is the objective (cost, energy, error), and \(c:\mathbb{R}^n\to\mathbb{R}^m\) is a vector of constraints — each one held between a lower bound \(l_i\) and an upper bound \(u_i\) (either may be \(\pm\infty\)).

Why it matters. Real problems — aircraft trajectory optimization, optimal power flow, chemical-plant scheduling, computable general-equilibrium economics, model-predictive control — routinely have thousands to millions of variables and constraints. The functions involved are often expensive (a single evaluation may run a physics simulation), so the solver's number of function evaluations is the real currency, not raw arithmetic speed. SNOPT is engineered around two facts: (1) its SQP method needs few function evaluations, and (2) most large problems have sparse structure — each constraint touches only a handful of variables — which sparse linear algebra can exploit. That combination is why it has survived four decades of industrial use.

Takeaway. SNOPT = constrained minimization where "large" and "expensive" are the hard parts, solved with few evaluations and sparse linear algebra.

2 · The problem form — one notation to rule them all

The "range" form \(l \le c(x) \le u\) is a deliberately uniform bookkeeping device. It absorbs three things an engineer thinks of separately:

Variable bounds

\(x_j^{\min} \le x_j \le x_j^{\max}\)

a "constraint" row that just returns \(x_j\)

Linear constraints

\(l_i \le a_i^{\top} x \le u_i\)

rows of a sparse matrix \(A\)

Nonlinear constraints

\(l_i \le g_i(x) \le u_i\)

smooth functions \(g_i\)

▼

\[ l \le c(x) \le u,\qquad c(x)=\big(\,x;\ Ax;\ g(x)\,\big) \]

Setting \(l_i=u_i\) gives an equality; \(u_i=+\infty\) gives a one-sided inequality; \(l_i=-\infty,\ u_i=+\infty\) gives an inactive (ignored) constraint.

Why the uniform range form? It lets one algorithm handle bounds, linear constraints, and nonlinear constraints identically — every constraint is just "a function that must land in an interval." The solver never needs to know which kind it is; it only needs each constraint's value \(c_i(x)\) and gradient \(\nabla c_i(x)\).

Two qualifications, stated up front because they define SNOPT's scope. First, the objective \(f\) may be linear or nonlinear, so SNOPT also covers linear programs (LP) and quadratic programs (QP) as special cases. Second, SNOPT finds local optima and is most reliable when \(f\) and \(c\) are smooth; it assumes you can supply first derivatives (or tolerates finite-difference estimates when you cannot).

Takeaway. Bounds, linear, and nonlinear constraints collapse into one ranged constraint vector \(l\le c(x)\le u\) — the solver only needs values and gradients.

3 · What "optimal" means — the KKT conditions

To build an algorithm, we need a precise, checkable statement of optimality. For the equality-constrained case \(c(x)=0\) (take \(l=u=0\) for now), the right object is the Lagrangian:

\[ L(x,\lambda) \;=\; f(x) + \lambda^{\top} c(x), \qquad \lambda\in\mathbb{R}^m \]

Why introduce \(\lambda\)? Unconstrained optimality says \(\nabla f(x^\star)=0\). With constraints, we can only move along the constraint surface, so the correct requirement is weaker: the objective's gradient must be balanced by the constraint gradients — i.e. \(\nabla f\) must be a combination of the \(\nabla c_i\). The multipliers \(\lambda_i\) are exactly the weights in that combination.

The first-order necessary (KKT) conditions at a local minimizer \(x^\star\) (assuming a standard constraint qualification holds) are:

\[ \nabla f(x^\star) + J(x^\star)^{\top}\lambda^\star = 0, \qquad c(x^\star)=0, \]

where \(J(x)\in\mathbb{R}^{m\times n}\) is the Jacobian of \(c\) (its \(i\)-th row is \(\nabla c_i(x)^{\top}\)). The first equation is stationarity; the second is feasibility. For inequalities \(c_i(x)\le 0\), add complementarity and sign conditions \(\lambda_i \ge 0,\ \lambda_i c_i(x)=0\): a multiplier can be nonzero only if its constraint is active.

\(x^\star\)

\(\nabla c\)

\(\nabla f\)

\(f(x)=f^\star\)

\(c(x)=0\)

Figure 1 — KKT geometry. The optimum \(x^\star\) is where the objective level set (blue dashed circle) is tangent to the constraint (solid circle). At that point \(\nabla f\) (blue) and \(\nabla c\) (light blue) are collinear and opposite: \(\nabla f + \lambda \nabla c = 0\) with \(\lambda > 0\).

Takeaway. Optimality is a balance of gradients — \(\nabla f\) must equal a combination of constraint gradients, with weights \(\lambda\) — plus feasibility (and, for inequalities, complementarity).

4 · The core idea — SQP is Newton's method on the KKT system

We have a system of \(n+m\) equations in \(n+m\) unknowns \((x,\lambda)\). Write it as \(F(x,\lambda)=0\):

\[ F(x,\lambda) \;=\; \begin{bmatrix} \nabla f(x) + J(x)^{\top}\lambda \\ c(x) \end{bmatrix} \;=\; 0. \]

Why Newton? Newton's method converges quadratically once close to a solution — it roughly doubles correct digits per iteration. That speed, not the raw cost of one step, is what keeps the function-evaluation count small, which (recall §1) is the whole point for expensive functions.

Newton's iteration solves \(F'(z_k)\,\Delta z = -F(z_k)\). Differentiating \(F\) with respect to \((x,\lambda)\) gives a saddle-point (KKT) matrix:

\[ \underbrace{\begin{bmatrix} H_k & J_k^{\top} \\ J_k & 0 \end{bmatrix}}_{F'(z_k)} \begin{bmatrix} \Delta x \\ \Delta \lambda \end{bmatrix} \;=\; -\begin{bmatrix} \nabla f_k + J_k^{\top}\lambda_k \\ c_k \end{bmatrix}, \]

where \(H_k = \nabla^2_{xx} L(x_k,\lambda_k) = \nabla^2 f_k + \sum_i (\lambda_k)_i \nabla^2 c_i(k)\) is the Hessian of the Lagrangian, and \(c_k = c(x_k)\). The trick that makes this an optimization method rather than a root-finder:

Why solve a QP instead of that linear system directly? Two reasons. (1) Inequalities. The Newton system above only handles equalities; a quadratic program can carry inequality constraints \(l\le J_k d + c_k \le u\) and its active-set solver decides which constraints bind at each step — that is exactly the hard combinatorial part of constrained optimization, delegated to a QP solver. (2) Equivalence. The KKT conditions of a QP are a linear system of the very same shape, so solving a QP and taking a Newton step are one and the same operation.

The QP subproblem at iteration \(k\) is:

\[ \min_{d}\ \ \tfrac12 d^{\top} H_k\, d + \nabla f_k^{\top} d \qquad \text{s.t.} \qquad l \le J_k\, d + c_k \le u . \]

Its solution \(d\) is the search direction, and its multiplier \(\mu\) is the new multiplier estimate (in the equality case, \(\mu = \lambda_k + \Delta\lambda\)). To see the equivalence, write the QP's stationarity condition \(H_k d + \nabla f_k + J_k^{\top}\mu = 0\) and feasibility \(J_k d + c_k = 0\); substituting \(d=\Delta x\) and \(\mu=\lambda_k+\Delta\lambda\) reproduces the Newton system exactly.

Start

\((x_k,\lambda_k)\)

▼

Evaluate

\(f,\ c,\ \nabla f,\ J\)

→

Build QP

\(\tfrac12 d^{\top}H d + \nabla f^{\top} d\) s.t. \(l \le J d + c \le u\)

▼

Solve QP

direction \(d\), multiplier \(\mu\)

→

Line search

step \(\alpha\) on merit function

▼

Update

\(x_{k+1}=x_k+\alpha d,\ \lambda_{k+1}=\mu\)

▼

Check KKT

converged? → stop, else loop

Figure 2 — One SQP iteration. Form a quadratic model of the Lagrangian with linearized constraints, solve it, move, repeat. Every ingredient here needs a careful choice (§5–§7); this diagram is just the skeleton.

Takeaway. SQP = Newton's method on the KKT conditions, implemented as: build a QP that quadratically models the Lagrangian under linearized constraints, solve it, step.

5 · Globalization — a merit function so it converges from anywhere

Raw Newton steps are only reliable near a solution: far away, a full step can move toward the constraint too aggressively or even increase the objective. We need a single scalar that tells us "is this new point better?", so we can shrink the step until it is. SNOPT uses an augmented-Lagrangian merit function. For equalities, its essential form is:

\[ M(x;\lambda,\rho) \;=\; f(x) + \lambda^{\top} c(x) + \tfrac{\rho}{2}\|c(x)\|_2^2 , \]

Why this form? It fuses two objectives that are often in tension. The term \(f+\lambda^{\top}c\) is the Lagrangian, which encodes "improve the objective as judged by the current multiplier." The term \(\tfrac{\rho}{2}\|c\|^2\) is a penalty that screams "stay feasible." The parameter \(\rho\) dials the balance; as \(\rho\) grows, minimizing \(M\) enforces feasibility harder. Crucially, \(M\)'s stationary point and gradient match the KKT conditions as \(c\to 0\), so the merit function and the true problem agree at a solution.

A line search then works along the QP direction \(d\): try \(\alpha=1\) (the full Newton step), and backtrack \(\alpha \leftarrow \alpha/2\) until the merit function decreases sufficiently:

\[ M(x_k + \alpha d) \;<\; M(x_k) + \eta\,\alpha\, M'(x_k; d), \qquad 0<\eta<1. \]

\(M(x_k)\)

accept: \(M(x_k+\alpha d)<M(x_k)\)

\(\alpha=1\) rejected

step length \(\alpha\) →

Figure 3 — A line search. The merit function \(M(x_k+\alpha d)\) is traced as the step length \(\alpha\) grows. The full step \(\alpha=1\) raises \(M\) (it overshoots the constraint), so we backtrack to the green point where \(M\) actually decreases.

For inequalities, SNOPT adds slack variables \(s\) and uses the augmented Lagrangian \(M = f - \lambda^{\top}(c-s) + \tfrac{\rho}{2}\|c-s\|^2\) (with \(s\) bounded by the constraint range), which handles one-sided constraints gracefully. This is the machinery behind the documentation's promise: "an augmented Lagrangian merit function ensures convergence from an arbitrary point."

Takeaway. Newton gives the direction; a merit function + line search gives the guarantee that every accepted step makes genuine progress, from any starting point.

6 · The Hessian — quasi-Newton, and limited-memory for scale

The QP needs \(H_k\), the Hessian of the Lagrangian. Two problems: computing second derivatives \(\nabla^2 f\) and \(\nabla^2 c_i\) is expensive and error-prone, and \(H_k\) may be indefinite (its QP then has no clean minimum). SNOPT sidesteps both by approximating \(H_k\) with a quasi-Newton matrix \(B_k\), built from gradient differences only.

Why gradient differences? The secant equation \(B_{k+1} s = y\) — with \(s = x_{k+1}-x_k\) and \(y = \nabla L_{k+1}-\nabla L_k\) — says "whatever \(B\) is, it must be consistent with the curvature the gradients actually observed along the last step." The BFGS update is the minimal change to \(B_k\) that satisfies the secant equation and preserves symmetry & positive definiteness:

\[ B_{k+1} = B_k \;-\; \frac{B_k s s^{\top} B_k}{s^{\top} B_k s} \;+\; \frac{y y^{\top}}{y^{\top} s}. \]

The payoff: the user supplies only first derivatives, and \(B_k\) gradually builds a curvature model of the Lagrangian as iterations proceed — typically enough for superlinear (near-quadratic) convergence.

For genuinely large \(n\), even storing the dense \(n\times n\) matrix \(B_k\) is impossible. The limited-memory BFGS (L-BFGS) variant keeps only the last \(r\) pairs \((s_j, y_j)\) and represents \(B_k\) implicitly, so storage is \(O(nr)\) instead of \(O(n^2)\). This is what lets the "sparse" solver also run when the Hessian is effectively dense.

Dense \(B\)

\(O(n^2)\) storage

exact-in-spirit curvature, but memory explodes for large \(n\)

Limited-memory \(B\)

\(O(n r)\) storage

last \(r\) curvature pairs, applied implicitly — scales to huge \(n\)

Exact \(H\) (optional)

user-supplied \(\nabla^2\)

best convergence if you can provide true second derivatives

Figure 4 — Three ways to get the Hessian term in the QP. SNOPT defaults to limited-memory quasi-Newton, but also accepts an exact Hessian (and, for problems where the Lagrangian Hessian is naturally sparse, a structured one).

Takeaway. The second-derivative term is approximated by BFGS from gradients alone — and L-BFGS keeps that approximation's memory cost linear, which is what makes the method "large-scale."

7 · Why "Sparse" is in the name

In a large model, each nonlinear constraint \(g_i\) typically involves only a few of the \(n\) variables — an aircraft's aerodynamics at one node depends on a handful of neighboring state variables, not all of them. The Jacobian \(J\) is therefore mostly zeros:

Sparse Jacobian

nonzeros \(\propto n\), not \(n^2\)

Dense (naive) Jacobian

nonzeros \(\propto n^2\) — infeasible at scale

Figure 5 — Sparsity is the difference between feasible and impossible. The KKT systems solved inside each QP inherit this structure; factorizing them with a sparse \(LU\) (SNOPT uses the LUSOL package) keeps work proportional to the number of nonzeros rather than \(n^3\).

Why does this matter to the algorithm? The QP subproblem's KKT system is solved by factorization. A dense factorization of an \(n\)-dimensional system costs \(O(n^3)\) — useless beyond a few thousand variables. A sparse factorization costs roughly \(O(\text{nonzeros})\), which is what turns "thousands of variables" into "millions."

Sparsity is also why SNOPT is most efficient in two specific regimes the authors call out: when only some variables enter nonlinearly, or when there are relatively few degrees of freedom at a solution (many constraints active). Both keep the effective factorizations small.

Takeaway. The "Sparse" in SNOPT is not cosmetic — exploiting Jacobian sparsity in the KKT factorizations is what separates a toy solver from an industrial one.

8 · Robustness machinery — elastic mode, scaling, warm start

Three practical features round out SNOPT from "the math works" to "it works in production."

Elastic mode — what if no feasible point exists?

A model with a typo in a bound is infeasible; a hard solver just fails and tells you nothing useful. SNOPT instead adds elastic variables \(v \ge 0\) that allow each constraint to be violated by \(v_i\), then minimizes the total violation. The result is the closest feasible point and a clear report of which constraints are broken — an invaluable debugging signal.

\[ \min_{x,v}\ f(x) + \rho \sum_i v_i \qquad \text{s.t.} \qquad l - v \le c(x) \le u + v,\ \ v \ge 0 . \]

Scaling

If one variable is measured in meters and another in nanometers, the raw problem is numerically hostile (the KKT matrices become ill-conditioned and Newton stalls). SNOPT applies automatic scaling of variables and constraints to keep the numerics tame — a small thing that silently saves a great many runs.

Cold vs. warm start

A cold start begins with no prior information. A warm start feeds the solver a previous solution's active set, basis, and multipliers — so re-solving a slightly perturbed problem (as in model-predictive control, or a sequence of design iterations) costs a fraction of a cold solve. Warm starting is a major reason SQP/active-set methods persist in time-critical control loops.

Takeaway. Elastic mode, scaling, and warm starts are the engineering layer that turns a correct algorithm into a dependable tool — especially for infeasible models and repeated solves.

9 · Worked example — the whole machine on a 2-variable problem

To make every moving part concrete, take a toy problem with a clean, hand-checkable answer:

\[ \min_{x}\ \ (x_1-2)^2 + (x_2-2)^2 \qquad \text{s.t.} \qquad x_1^2 + x_2^2 = 1. \]

Geometrically: find the point on the unit circle closest to \((2,2)\) — obviously the diagonal point \((1/\sqrt2,\,1/\sqrt2)\), but let's derive it and then watch SQP find it numerically.

Step 1 — KKT by hand. With \(L = f + \lambda(x_1^2+x_2^2-1)\), stationarity gives

\[ 2(x_1-2) + 2\lambda x_1 = 0, \qquad 2(x_2-2) + 2\lambda x_2 = 0 . \]

Both equations say \(x_j(1+\lambda) = 2\), so \(x_1 = x_2 = \dfrac{2}{1+\lambda}\). The constraint \(x_1^2+x_2^2=1\) then forces \(2\left(\dfrac{2}{1+\lambda}\right)^2 = 1\), i.e. \((1+\lambda)^2 = 8\). Taking the positive branch (the closest point):

\[ \lambda^\star = 2\sqrt2 - 1 \approx 1.828, \qquad x^\star = \Big(\tfrac{1}{\sqrt2},\ \tfrac{1}{\sqrt2}\Big) \approx (0.7071,\ 0.7071). \]

The objective there is \(f^\star = 9 - 4\sqrt2 \approx 3.343\). (The other branch \(\lambda = -1-2\sqrt2\) is the far point \((-1/\sqrt2,-1/\sqrt2)\) — also a KKT point, but a local maximum of \(f\) on the circle, a reminder that SQP finds local KKT points, not the global one.)

Step 2 — watch it iterate. The demo below runs exactly the SQP loop of §4–§5 on this problem. Click the plot to set a starting point (or use a preset), then Step or Run. The Hessian of the Lagrangian is computed exactly here (\(H = 2(1+\lambda) I\)) so you can watch pure SQP, unclouded by quasi-Newton noise.

Iteration	0
\(x_1,\ x_2\)	2.000, 2.000
\(\lambda\)	0.000
\(f(x)\)	0.000
\(c(x)=x_1^2+x_2^2-1\)	7.000
step \(\alpha\)	—
‖KKT residual‖	—

The path converges to \(x^\star\approx(0.7071,0.7071)\), \(\lambda^\star\approx1.828\), matching the hand solution. This is the algorithm illustrating itself — not a proof of anything, just the iterations made visible.

Takeaway. A handful of SQP iterations (each one QP-solve + line search) lands on the KKT point to several digits — the quadratic convergence promised by Newton.

10 · Where SNOPT fits — active-set SQP vs. interior-point, and the family tree

Modern large-scale NLP gives you two families of solvers, and they make different tradeoffs:

Active-set SQP

SNOPT, NPSOL, filterSQP

• few function evaluations
• excellent warm starts
• decides the active set explicitly
• strong with many active constraints

Interior-point

IPOPT, KNITRO, LOQO

• polynomial complexity bound
• very robust on huge, degenerate problems
• less natural warm-starting
• great when many inequalities are active

Figure 6 — Two ways to handle inequality-constrained NLP. Neither dominates; the choice turns on whether you need few evaluations/warm starts (SQP) or worst-case scale and robustness (interior-point).

SNOPT descends from a celebrated lineage at Stanford/UCSD built by Philip Gill, Walter Murray, and Michael Saunders: MINOS (reduced-gradient, 1970s) → NPSOL (dense SQP) → SNOPT (sparse SQP) → SQOPT (sparse QP). The definitive reference is their 2005 SIAM Review paper, "SNOPT: An SQP algorithm for large-scale constrained optimization."

Why reach for SNOPT specifically? The documentation's own summary is the best guidance: SNOPT is "especially effective if the objective or constraint functions (and their gradients) are expensive to evaluate," when "only some variables enter nonlinearly," or when there are "relatively few degrees of freedom at a solution." In those regimes its low evaluation count and sparse factorization shine; for a dense, huge, highly-nonlinear, many-active-inequality problem, an interior-point solver may be the safer default.

Takeaway. SNOPT is the sparse-SQP specialist: few evaluations, superb warm starts, and sparsity-aware linear algebra — ideal for expensive, structured, constraint-rich problems.

11 · Summary

SNOPT solves \(\min_x f(x)\) s.t. \(l \le c(x) \le u\) for large, sparse, smooth problems. Its engine is SQP: repeatedly solve a quadratic model of the Lagrangian under linearized constraints, which is exactly a Newton step on the KKT conditions. A quasi-Newton (L-BFGS) approximation supplies the Hessian from gradients alone; an augmented-Lagrangian merit function with a line search guarantees progress from any start; sparse \(LU\) factorizations make it scale. Elastic mode, scaling, and warm starts make it dependable in practice.

The one-line mental model. SQP = "quadratically approximate the Lagrangian, linearize the constraints, solve, step, repeat" — with sparse linear algebra doing the heavy lifting.