Sequential Quadratic Programming

Reading: Biegler (2010), Nonlinear Programming, §6.2, pp. 135–151; Nocedal & Wright (2006), Numerical Optimization, Chapter 18, pp. 529–562.

This guide supplements Lecture 23 with complete QP certificates, elastic repairs, line-search and trust-region algorithms, reduced-space details and local-rate conditions. The implementation laboratory then exposes each dense QP/KKT solve and preserves successful runs and useful failure cases.

Start with the QP model, elastic repairs, trust-region families, or local rates. Constrained Globalization supplies the general acceptance proofs and Maratos correction. Equality Newton supplies the equality-QP bridge; Quasi-Newton Methods supplies the general BFGS derivation.

The prototype implements scalar-relaxed line-search SQP. The trust-region, SLQP and sparse methods are fully described here but require different computational machinery; they are not implemented by the prototype.

Formulating the inequality QP¶

On an open $D\subset\mathbb R^n$, take $f:D\to\mathbb R$, $\mathbf h:D\to\mathbb R^{m_h}$ and $\mathbf g:D\to\mathbb R^{m_g}$ of class $C^2$. Row Jacobians have dimensions $m_h\times n$ and $m_g\times n$; the curvature model $\mathbf B$ is symmetric $n\times n$. Use the course convention $L=f+\mathbf v^{\mathsf T}\mathbf h+\mathbf u^{\mathsf T}\mathbf g$, with equalities $\mathbf h=\mathbf0$, inequalities $\mathbf g\le\mathbf0$, and $\mathbf u\ge\mathbf0$. At iteration $k$, put $\mathbf a_k=\nabla f(\mathbf x^k)$ and $\mathbf B_k=\nabla_{\mathbf{xx}}^2L(\mathbf x^k,\mathbf v^k,\mathbf u^k)$ or a chosen approximation. An IQP subproblem is

This uses a quadratic Lagrangian model and affine constraint models. Linearizing a nonlinear inequality does not certify feasibility of the nonlinear trial. The QP supplies a step and multiplier estimates; an acceptance procedure from Constrained Globalization selects the new point. Equality Newton owns the exact equality KKT/QP bridge and its local assumptions.

In the bound formulation used by the code, introduce nonnegative slacks with $\mathbf g(\mathbf x)+\mathbf s=\mathbf0$ and relabel the enlarged primal vector and equality vector. The equality Jacobian below is $\mathbf J$. Inequalities are now represented by these equalities and simple bounds. The subproblem is

Here the code’s argument g means the objective gradient, its J is a row Jacobian, and its lam means the course equality multiplier $\mathbf v$. Its bound-QP solver returns a signed gradient zeta: positive at a lower bound, negative at an upper bound. At a solved QP, $\boldsymbol\zeta=\mathbf u_L-\mathbf u_U$ and stationarity is $\mathbf a+\mathbf B\mathbf d+\mathbf J^{\mathsf T}\bar{\mathbf v}-\mathbf u_L+\mathbf u_U=\mathbf0$. Do not read this signed quantity as a nonnegative inequality multiplier.

Biegler §6.2, pp.135–136; Nocedal–Wright (18.11), p.533.

Function signs and multiplier signs¶

Nocedal–Wright writes $L_{NW}=f-\boldsymbol\lambda_E^{\mathsf T}\mathbf h-\boldsymbol\lambda_I^{\mathsf T}\mathbf c_I$ with $\mathbf c_I\ge\mathbf0$. For the same equality function, $\mathbf v=-\boldsymbol\lambda_E$. When $\mathbf g=-\mathbf c_I$, however, $\mathbf u=\boldsymbol\lambda_I$. Thus equality multipliers reverse sign but inequality multipliers do not under this function conversion. For a lower bound, $c_L=x-x_L$ and $g_L=x_L-x$; for an upper bound, $c_U=x_U-x$ and $g_U=x-x_U$. Both nonnegative bound multipliers retain their values. State the constraint function and Lagrangian together when translating.

If a bound set were known, matrices $\mathbf E_L,\mathbf E_U$ collect the corresponding unit columns. The QP would impose $\mathbf E_L^{\mathsf T}(\mathbf x^k+\mathbf d)=\mathbf E_L^{\mathsf T}\mathbf x_L$ and its upper-bound analogue. The general inequality QP keeps the bounds as inequalities and chooses which are active at its solution. Active bounds at the current point and at the QP endpoint need not be the same.

A correctly scoped zero-step certificate¶

Flat QP minimizers span more than the selected endpoint face.

If a QP KKT solution has $\mathbf d=\mathbf0$, its affine residuals equal the original constraint values. QP stationarity becomes NLP stationarity, and dual feasibility/complementarity apply at the same point. Therefore a zero QP KKT step gives an NLP KKT point, without needing a separate uniqueness argument. It does not alone prove second-order optimality or a global minimum. In computation, verify the QP’s KKT residuals and solve status rather than assuming its returned step is exact.

A useful converse holds with sufficient strict convexity on the equality tangent. Suppose the current point is NLP KKT and $\mathbf B$ is positive definite on $\ker\mathbf J_h$. Zero is a feasible QP step. For any other feasible $\mathbf d$, stationarity gives

The equality term is zero. Positive multipliers only occur at active inequalities; there $\mathbf J_{g,i}\mathbf d\le0$, so each remaining term is nonnegative. Since $\mathbf J_h\mathbf d=0$ and $\mathbf d\ne0$, $\frac12\mathbf d^{\mathsf T}\mathbf B\mathbf d>0$. Hence the QP objective at every nonzero feasible step exceeds its value at zero, making zero the unique global QP minimum. Full positive definiteness is a convenient sufficient special case. This converse imposes more than curvature on one selected all-active face.

The distinction matters in Biegler Theorem 6.1, p.136. Consider $f=x_2^2$, no equalities, $0\le x_1\le1$, current $(0,0)$, and QP Hessian $\operatorname{diag}(0,2)$. The QP may return $(1,0)$ with zero upper-bound multiplier. At that endpoint the active row $[1,0]$ is independent, the null curvature is2, and its KKT matrix is nonsingular. Yet all $(d_1,0)$ with $0\le d_1\le1$ minimize the QP. The original current point is NLP KKT with the same zero multipliers, despite the selected nonzero step. Positivity on the all-active null space misses the weak bound’s critical directions. In the book’s fixed-face proof, zero does not satisfy that endpoint’s upper-bound equality. Use the certificate/converse above rather than promoting active-face nonsingularity to global primal uniqueness. Second Order Conditions supplies the critical-cone distinction.

A sharper fixed-face converse can instead assume strict complementarity at the selected QP endpoint, the same multipliers satisfying NLP KKT at the current point, and positive curvature on that selected active null space. Every positively weighted bound/inequality is active at both points, so its affine row gives $\mathbf J_A\mathbf d=\mathbf0$. Subtracting the two stationarity equations gives $\mathbf B\mathbf d=\mathbf0$. Therefore $\mathbf d^{\mathsf T}\mathbf B\mathbf d=0$ and active-null positivity forces $\mathbf d=\mathbf0$. The flat example fails this added strict-complementarity assumption. This is a different sufficient converse, not a reason to omit weakly active directions from the general theorem.

Active-set identification and well-posedness¶

Strictly and weakly weighted bounds respond differently to center perturbations.

At a feasible $\mathbf x^*$, define $\mathcal A^*=\{i:g_i(\mathbf x^*)=0\}$. The active Jacobian stacks $\mathbf J_h$ and the rows $\nabla g_i^{\mathsf T}$ for $i\in\mathcal A^*$. LICQ means these rows are linearly independent; strict complementarity means $u_i^*>0$ for every $i\in\mathcal A^*$. Under strict complementarity, SOSC uses the true solution Hessian $\mathbf W^*=\nabla_{\mathbf{xx}}^2L(\mathbf x^*,\mathbf v^*,\mathbf u^*)$ and requires $\mathbf p^{\mathsf T}\mathbf W^*\mathbf p>0$ for every nonzero direction in the null space of this active Jacobian. L18/L21 supply the general critical-cone and second-order statements.

For the exact SQP curvature $\mathbf B=\nabla_{\mathbf{xx}}^2L(\mathbf x,\mathbf v,\mathbf u)$, Nocedal–Wright Theorem 18.1, p.533, assumes a local KKT solution, LICQ, strict complementarity and SOSC. There is a neighborhood of the solution and its relevant multiplier vector where a local QP solution branch has the same active set. A general Hessian approximation needs its own closeness/regularity conditions; the exact-Hessian theorem does not apply to an arbitrary symmetric model. Biegler Theorem 6.5, pp.141–142, gives the analogous eventual-set property for a converging SQP sequence. This stabilizes the active-set algebra; it does not prove global QP uniqueness far away or remove the need for accepted full steps and local smoothness/Hessian assumptions.

Strict complementarity is sufficient for this identification theorem, not necessary for every example to identify its set. For $\min x^2$ subject to $x\ge0$, the solution has zero multiplier, but the exact QP step from a positive start reaches zero immediately and its active set remains identifiable. Without strict complementarity, the theorem no longer excludes changing weakly active sets; it does not assert that all such problems fail.

Likewise, LICQ is a sufficient KKT factorization condition, not a necessary condition for unique primal optimization: $\min x^2$ with duplicate inequalities $-x\le0,-2x\le0$ has a unique minimum. A solver may select an independent working subset. It must still distinguish truly redundant rows from dependent gradients with inconsistent residuals: the equations $d=0$ and $d=-1$ have the same gradient but cannot both be satisfied. Discarding the latter merely because its gradient is dependent changes the problem. Preserve primal/dual residual checks when rank safeguards or least-squares solves are used.

Three separate questions govern the subproblem: is its feasible set nonempty, is its objective bounded below/attained, and which primal/dual solution is unique? Full indefiniteness need not make an equality QP unbounded when reduced curvature is positive. A compact inequality feasible set bounds a continuous quadratic objective even if it is indefinite. Independent active rows with positive curvature on that face make its KKT system nonsingular; they need not settle global minimization across other faces. These are distinct issues.

IQP, EQP and SLQP¶

An IQP outer method solves an inequality QP and reads a working-set estimate from its solution. An EQP outer method estimates a working set first, then computes an equality-QP step on that set. An active-set QP subsolver can be used inside IQP; choosing it does not convert the outer algorithm to EQP. Nocedal–Wright §18.2, pp.533–534.

Warm starts reuse a previous QP working set or factorization when successive subproblems are related. They can reduce work near a stable solution; they do not guarantee that the next set is identical. Whether IQP or EQP is cheaper depends on derivative evaluation cost, sparsity, conditioning, degeneracy and the quality of the set estimate. Neither family requires exhaustive enumeration of every inequality subset.

SLQP uses an LP phase to estimate a working set, followed by an EQP phase. A hard-constraint LP can itself be infeasible, so the practical scheme uses a penalty LP on an infinity-norm trust box. Its objective is the linear objective model plus the piecewise-linear constraint violation. A basic optimal solution can identify a working set, but not every LP optimum returned by every solver is a vertex. Select independent active rows and retain feasibility/restoration safeguards.

The fuller Nocedal–Wright construction, pp.551–552, uses a Cauchy point $\mathbf p^C=\alpha_{LP}\mathbf p^{LP}$ giving sufficient decrease in the merit model. It selects an independent subset of constraints active in the LP model. In the EQP objective, violated constraints contribute the derivative of their current penalty branch: equality residuals contribute $\rho_{\rm pen}\operatorname{sign}(h_i^{lin})\nabla h_i$, and violated course inequalities $g_i^{lin}>0$ contribute $\rho_{\rm pen}\nabla g_i$. Working-set rows are enforced at equality inside a Euclidean trust ball with its own radius. Let $\mathbf p^Q$ be that candidate. The total step is

chosen to improve the merit model. This is not simply discarding the LP point and taking an unconstrained equality-QP solution. LP and EQP radii serve different roles. The penalty update below also applies, with a zero-Hessian LP model.

Inconsistent linearizations: complete worked example¶

A feasible nonlinear point can coexist with contradictory affine step requirements.

The constraints $x\le1,x^2\ge4$ are feasible at $x=-2$, but at $x=1$ their linearizations require $d\le0$ and $d\ge3/2$. A feasible NLP can have an empty local affine model inside its physical bounds.

Biegler Example 6.3, p.140, uses

with $0\le x_2\le2$ and $x_3,x_4\ge0$. At $(0,.1,0,0)$ both residuals are .99, not the printed first value 1.01. The affine rows are

with $-.1\le d_2\le1.9$ and $d_3,d_4\ge0$. Adding the rows gives $.4d_2=1.98+d_3+d_4\ge1.98$, hence $d_2\ge4.95$, contradicting its upper bound. The held errata restores the omitted starting-point components but does not repair the 1.01 arithmetic.

The global solution is $(-3,2,6,0)$, with objective -1. A complete verification eliminates the equalities: nonnegative slacks imply $1-x_2^2\le x_1\le x_2^2-1$, so $1\le x_2\le2$. At a fixed $x_2$, minimize by $x_1=1-x_2^2$, giving $f=1-x_2^2+x_2$. Its derivative $1-2x_2<0$ on $[1,2]$, so the unique minimizing $x_2$ is2 and the remaining variables follow. Course equality multipliers $(0,1)$, upper $x_2$ multiplier3 and lower $x_4$ multiplier1 satisfy KKT. This global calculation uses this example’s special structure; a generic SQP status is not a global certificate.

Two elastic repairs and what relaxation means¶

Residual slacks and their multiplier cap.

The $\ell_1$ elastic QP introduces $\mathbf s,\mathbf t\ge\mathbf0$:

For a bound-feasible current point, zero step and component positive/negative parts of $\mathbf h$ give a feasible lifted point. More generally, a nonempty physical box allows a bound-feasible step and appropriate residual slacks. Strict positive primal curvature gives a finite unique primal minimum; mere feasibility with an indefinite unbounded objective does not. At optimality with positive weight, the slacks represent the residual’s actual positive and negative parts. Their stationarity implies $-\rho_{\rm pen}\le\bar v_j\le\rho_{\rm pen}$; a positive $s_j$ forces $\bar v_j=\rho_{\rm pen}$ and a positive $t_j$ forces $\bar v_j=-\rho_{\rm pen}$. This multiplier cap belongs to the l1 formulation, not the scalar formulation. Biegler6.13p138; Nocedal–Wright18.12p536.

The scalar repair implemented below minimizes

with the physical step bounds and $0\le\xi\le1$. The book’s6.16only prints $\xi\ge0$; the notebook adds the upper bound. A zero step with $\xi=1$ is feasible when the current point satisfies the box. For $\mathbf B\succ0,M>0$, the lifted objective Hessian $\operatorname{diag}(\mathbf B,M)$ is positive definite, giving a unique minimizer on the convex feasible set. The $\xi^2$ term alone does not make flat primal directions unique.

At $\xi=0$ the chosen step satisfies the original affine equalities. For $0<\xi<1$ it contracts their residual by $1-\xi$. A positive $\xi$ does not prove that exact affine feasibility was impossible: for $h=-1,J=1,a=0,B=1000,M=1$, and physical $x\in[0,2]$ at $x=0$, the exact step1is feasible but the weighted optimum chooses $\xi=999/1001$ and $d=2/1001$. A finite weight can prefer partial relaxation. With only the book’s lower bound on $\xi$, the data $a=100,B=M=1,h=-1,J=1$ and broad physical bounds produce $\xi=50,d=-49$. Thus a contraction interpretation needs an explicit domain/safeguard.

At $\xi=1$ the step is tangent to the affine equalities; it need not be zero. For $h=x_1-x_1^2+1$, current $(0,0)$, lower bound $x_1\ge0$, $\mathbf a=(0,-1)$ and $\mathbf B=\mathbf I$, the guarded scalar QP forces $\xi=1,d_1=0$ but selects $d_2=1$. A feasible nonlinear point exists at $x_1=(1+\sqrt5)/2$. Diagnose lost contraction/restoration needs; do not infer global infeasibility or zero tangent movement.

For Example 6.3, the relaxed row sum gives $1.98(1-\xi)\le.76$, so the minimum feasible relaxation is $1-.76/1.98=.616161\ldots$. The actual solver below computes its selected relaxation separately. At the all-zero start, $\xi=1$ is forced by the two rows and nonnegative slacks, but there remain feasible $d_2$ values in $[0,2]$. For this example’s gradient and positive Hessian, the minimizing primal step is zero. Preserve that specific failure; do not turn it into a universal statement about $\xi=1$.

Elastic merit descent: complete sign argument¶

The elastic merit-descent argument uses the actual contraction and bound signs.

Assume a bound-feasible current point and an exact QP KKT solution, including nonnegative bound multipliers and endpoint complementarity. For the actual scalar step, $\mathbf J\mathbf d=-(1-\xi)\mathbf h$. The forward norm expansion from Constrained Globalization gives

At the QP endpoint, bound complementarity and a bound-feasible current point imply $\bar{\mathbf u}_L^{\mathsf T}\mathbf d\le0$ and $\bar{\mathbf u}_U^{\mathsf T}\mathbf d\ge0$. Dot QP stationarity with $\mathbf d$ and substitute its equality row:

For $0\le\xi<1$, $1\le p\le\infty$, dual exponents $1/p+1/q=1$ with $1/\infty=0$, and $\rho_{\rm pen}\ge\|\bar{\mathbf v}\|_q+\epsilon_{\rm pen}$ with $\epsilon_{\rm pen}>0$ give

With $\mathbf B\succ0$, a nonzero step is strictly descending; at zero residual/zero step it is stationary, not a strict-descent exception. Biegler6.21p140prints a positive quadratic sign followed by a strict-negative claim; the preceding stationarity calculation requires the negative sign above. Full positive curvature is used here. The weaker reduced-curvature variants require the appropriate coupling-aware analysis from L22, not this unchanged multiplier rule.

A safeguarded line-search algorithm¶

A complete scalar-relaxed line-search outline is:

1. Start inside a consistent physical box and the differentiability domain $D$; evaluate every trial inside $D$. Evaluate $f,\mathbf a,\mathbf h,\mathbf J$ and the chosen curvature model.

2. Solve the guarded scalar QP, checking its primal/dual residuals and inner solve status. If it produces no adequate contraction or progress, diagnose restoration or a stall rather than declaring the entire NLP infeasible.

3. Compute the actual forward merit derivative. For the full-positive model and $\xi<1$, use the multiplier bound above; alternatively use a properly safeguarded direct slope rule. Compare all trials with the same current weight.

4. Backtrack from a full step with $0<\eta<1/2$ until $\phi_p(\mathbf x+\alpha\mathbf d)\le\phi_p(\mathbf x)+\eta\alpha D\phi_p$. Use L22’s correction safeguard near the solution. If a correction fails, discard it; if the line search exhausts its limit, report failure/restoration rather than unconditionally accepting its last trial.

5. Update the primal and specified dual estimates. Re-evaluate actual new-point feasibility, stationarity, dual signs and complementarity; a small step alone is not convergence. Then form the fixed-multiplier secant pair when updating curvature.

Biegler Algorithm 6.1p141takes full QP multiplier updates independently of primal $\alpha$. Nocedal–Wright Algorithm 18.3p545damps both: $\boldsymbol\lambda^{k+1}=\boldsymbol\lambda^k+\alpha(\hat{\boldsymbol\lambda}-\boldsymbol\lambda^k)$. They are different algorithm choices. Neither implies intermediate complementarity automatically.

Nocedal–Wright’s model-based rule uses $V(\mathbf d)=\|\mathbf h+\mathbf J\mathbf d\|_1$ and $q_\rho=f+\mathbf a^{\mathsf T}\mathbf d+(\sigma/2)\mathbf d^{\mathsf T}\mathbf B\mathbf d+\rho V(\mathbf d)$. Require $q_\rho(\mathbf0)-q_\rho(\mathbf d)\ge\kappa\rho[V(\mathbf0)-V(\mathbf d)]$, $0<\kappa<1$. If the violation decrease $\Delta V>0$, this gives

For an exact affine-feasibility step $\Delta V=\|\mathbf h\|_1$, this is18.36p542. For a relaxed step use its actual $\Delta V$; do not divide by zero at an already feasible tangent step. Algorithm 18.3uses $\sigma=1$. This rule measures predicted progress without relying directly on the multipliers; its behavior depends on the problem and model, not an unconditional superiority claim.

The existing code adds its own hybrid weight safeguard when relaxation is large. It is not literally either textbook algorithm in every branch. Its original stopping/failure behavior remains below, with separate diagnostic interpretation.

What the global argument must establish¶

Uniform decrease, summability and limit transfer are separate proof steps.

The term global convergence refers to approach to first-order stationary/KKT points from a larger class of starts, under algorithmic assumptions. It does not mean the globally best solution of a nonconvex NLP.

A transparent sufficient argument assumes bound-feasible primal iterates in a compact subset of $D$, QP solutions satisfying the scalar model’s KKT equations, bounded QP multipliers (including bound multipliers), an eventually fixed weight satisfying $\rho_{\rm pen}\ge\|\bar{\mathbf v}^k\|_q+\epsilon_{\rm pen}$ with one fixed $\epsilon_{\rm pen}>0$, accepted step lengths $\alpha_k\ge\underline\alpha>0$, $\mathbf B_k\succeq b\mathbf I$ with $b>0$ and uniformly bounded norm, and contraction $\xi_k\le1-\gamma$ with $\gamma>0$. Continuity on that compact set gives a merit lower bound by Weierstrass. The derivative bound above and accepted Armijo steps with one fixed $\eta>0$ yield

Summing against the merit lower bound gives finite sums of squared steps and violations, hence $\mathbf d_k\to\mathbf0$ and $\mathbf h_k\to\mathbf0$. Bounded curvature makes the QP stationarity defect $-\mathbf B_k\mathbf d_k$ vanish. Compactness and bounded duals allow subsequences with convergent primal/dual values; continuity gives stationarity and dual signs. QP complementarity at $\mathbf x^k+\mathbf d_k$ transfers to $\mathbf x^k$ because the difference in each product is a bounded multiplier times a vanishing step. Thus the limit is NLP KKT. This proof does not establish one convergent whole sequence, SOSC or a global minimum.

Biegler Theorem 6.4p141states a related KKT limit-point result for compact primal/trial sequences, uniform continuity, positive bounded curvature, LICQ, solvable QPs with $\xi<1$ and a sufficiently large weight. Its precedingp140discussion separately emphasizes contraction bounded away from1. Pointwise $\xi<1$ or pointwise positivity is not a uniform constant; an implementation must establish the safeguards its proof uses. Nocedal–Wright18.3p557explicitlyincludes bounded steps/dual values, a fixed sufficiently large penalty and its particular model-decrease line search. The prototype’s status string does not verify those global hypotheses. The finite-update issue for a bounded monotone weight is explained fully at the L22 reference.

Fixed-multiplier curvature and damping¶

Use $\mathbf s_k=\mathbf x^{k+1}-\mathbf x^k$ and

Both gradients fix the same updated multipliers. For linear physical bounds, the bound-gradient contributions cancel; the code therefore uses objective and equality-gradient terms. Changing the multiplier in the second term instead introduces an unwanted dual-step contribution. The general BFGS and damping derivations remain at Quasi-Newton Methods; this section applies them to constrained curvature.

For $\mathbf B\succ0,\mathbf s\ne\mathbf0$, put $a=\mathbf s^{\mathsf T}\mathbf B\mathbf s>0$ and $b=\mathbf s^{\mathsf T}\mathbf y$. Powell’s rule takes $\theta=1$ if $b\ge.2a$, otherwise $\theta=.8a/(a-b)$, and $\mathbf r=\theta\mathbf y+(1-\theta)\mathbf B\mathbf s$. The actual rule chooses $0<\theta\le1$; zero is a limiting conceptual case, not a finite selected value. If it selects1, the positive-curvature test has passed, so that actual ordinary update preserves definiteness. Forcing1while ignoring the test can produce an indefinite update. The convex combination is in secant vectors; updated matrices are not generally their affine mixture.

A constrained family makes negative measured curvature concrete. Let $f=-2x_1^2+\frac12(x_2-2)^2$, $h=x_1$, with solution $(0,2)$. Its true full Hessian is $\operatorname{diag}(-4,1)$ but its tangent curvature is1. At $(-\epsilon,2-3\epsilon)$ with $\mathbf B=\operatorname{diag}(2,3)$, the equality-QP step is $(\epsilon,\epsilon)$ and the fixed-multiplier secant is $(-4\epsilon,\epsilon)$. Thus $\mathbf s^{\mathsf T}\mathbf y=-3\epsilon^2$ even arbitrarily close to the strict constrained minimum. At $\epsilon=1$, $a=5,b=-3,\theta=.5,\mathbf r=(-1,2)$ and the actual existing damping routine yields $\left[\begin{smallmatrix}2.2&-3.2\\-3.2&5.2\end{smallmatrix}\right]$ with positive eigenvalues. A full positive approximation cannot converge to this specifically indefinite true Hessian; SOSC does not imply every problem has an indefinite full Hessian.

Reduced space and physical bounds¶

For full-row-rank $\mathbf J\in\mathbb R^{m\times n}$ and a null basis $\mathbf Z\in\mathbb R^{n\times(n-m)}$, use $\mathbf d=\mathbf n+\mathbf Z\mathbf q$, with $\mathbf J\mathbf n=-\mathbf h$ and $\mathbf J\mathbf Z=\mathbf0$. The unrelaxed reduced objective is

where $\mathbf H_R=\mathbf Z^{\mathsf T}\mathbf B\mathbf Z$. Its bounds are $\mathbf x_L\le\mathbf x+\mathbf n+\mathbf Z\mathbf q\le\mathbf x_U$. An arbitrary null basis does not turn these into independent coordinate bounds on $\mathbf q$. With $\mathbf J=[1,1]$, $\mathbf Z=(1,-1)$ and current $(.5,.5)$ inside $[0,1]^2$, the true tangent coordinate interval is $[-.5,.5]$; $q=1$ violates both physical bounds.

For scalar relaxation, the reconstructed step must be $\mathbf d=(1-\xi)\mathbf n+\mathbf Z\mathbf q$. Exact substitution into the full lifted objective retains

The normal and coupling terms depend on $\xi$ and cannot be dropped as constants. Biegler6.33p145is a reduced approximate model with frozen/dropped terms, not this exact substitution. Its subsequent Algorithm 6.3Step 4prints an unscaled normal reconstruction. Use the actual relaxed normal in a bound-feasible implementation and test its actual merit derivative.

The difference can survive line search: for $h=x_1+20x_1^2-1$, current $(0,0)$ and $0\le x_1\le.25$, the normal is $(1,0)$. A reduced elastic candidate $q=0,\xi=.75$ has valid scaled step $(.25,0)$. The unscaled $(1,0)$ violates the upper bound, although the NLP has feasible $x_1=.2$. A backtracked $\alpha=.275$ trial has violation .7875, below the weight1Armijo bound .9725 with $\eta=.1$, but still violates the physical bound by .025. Merit acceptance cannot substitute for preserving hard bounds. The original full lifted notebook solver already returns its actual primal step and is kept unchanged.

The exact multiplier from a full QP obeys

Dropping $\mathbf B\mathbf d$ produces a first-order/least-squares estimate, not the exact QP multiplier. With $\mathbf Y=\mathbf J^{\mathsf T}$ it uses $-(\mathbf J\mathbf J^{\mathsf T})^{-1}\mathbf J(\mathbf a-\mathbf u_L+\mathbf u_U)$. Nocedal–Wright pp. 539–540 also analyzes dropping the normal/tangent coupling term near a solution; it is an approximation justified by normal-step behavior and the rate conditions below, not an algebraic identity. Reduced BFGS uses $\mathbf s_R=\alpha\mathbf q$ and the projected fixed-multiplier gradient change. If the basis changes, its curvature coordinates need consistent transport/interpretation. L21’s equality reference owns the full normal/coupling derivation.

Large sparse factorization and active-set updates¶

Reduced-space methods can be economical when $n-m$ is small, but basis construction, dense fill, the number of active bounds and conditioning also matter. Many degrees of freedom can favor sparse full KKT/Hessian operations or limited-memory updates rather than a dense full BFGS matrix. Biegler pp. 144–147 uses examples and approximate dimension thresholds; these are design heuristics, not universal cutoffs or a claim that every QP implementation is dense.

An active-set solver can reuse a factorization when adding a row. Let $\mathbf K\mathbf z=\mathbf r$ be the current nonsingular KKT system and append $\mathbf e^{\mathsf T}\mathbf z=b$ with multiplier $u$. Solve $\mathbf K\mathbf z_0=\mathbf r$ and $\mathbf K\mathbf w=\mathbf e$. If $\mathbf e^{\mathsf T}\mathbf w\ne0$,

This is the Schur update in6.36p147; it uses solves, not explicitly forming an inverse. A zero or small denominator warns about dependence/conditioning and requires a consistency/rank safeguard. To drop an active row, permit a slack in that row and impose zero on its multiplier, yielding the analogous augmented6.37system. These operations update the working system, not a proof of global primal uniqueness. The prototype below deliberately recomputes dense systems rather than implementing these sparse updates.

Bound-aware trust-region families¶

Independent coordinate bounds need not enforce the total physical trust radius.

For $\Delta>0$, define $\|\mathbf d\|_\infty=\max_{1\le i\le n}|d_i|$. Its closed radius-$\Delta$ ball is the coordinate box $[-\Delta,\Delta]^n$. A hard SQP trust model adds $\|\mathbf d\|_\infty\le\Delta$ to the affine constraints and physical bounds. Its feasible set can be empty even for a feasible NLP. For the line $d_1+d_2=2.04$, its distance is 1.442 in Euclidean norm and 1.02 in infinity norm; both miss radius 1, but the ball and box are different sets. State the chosen norm when interpreting a circle illustration. A finite trust radius makes a feasible quadratic model bounded, including indefinite curvature; it does not automatically cure rank deficiency, infeasibility or inaccurate inner solves.

Composite relaxation. The generic normal/tangent proof and acceptance mechanism remain at Constrained Globalization. For bounds, retain $\mathbf x_L\le\mathbf x+\mathbf n+\mathbf Z\mathbf q\le\mathbf x_U$ in the tangent model and coordinate the total trust radius. Biegler6.39–6.45p149–150restricts the problem to dependent variables $\mathbf y$ without bounds and independent variables $\mathbf z$ with simple bounds. With $\mathbf J=[\mathbf J_y,\mathbf J_z]$ and invertible $\mathbf J_y$, it uses $\mathbf Y=[\mathbf I;\mathbf0]$ and $\mathbf Z=[-\mathbf J_y^{-1}\mathbf J_z;\mathbf I]$. Independent coordinate bounds on $\mathbf z+\mathbf q$ then have their stated meaning. These bases are not generally orthogonal, so the independent coordinate radius must be coordinated with $\|\mathbf Y\mathbf p_Y+\mathbf Z\mathbf q\|\le\Delta$.

For $h=y+2z-1$, current $(1.3,0)$, normal $(-.3,0)$ and $\mathbf Z=(-2,1)$, $q=.2$ gives total step $(-.7,.2)$ with length .7280. It satisfies an independent $z$ bound $|q|\le.25$ but violates a total step radius .6. The true squared length is $.09+1.2q+5q^2$, including coupling. This special-coordinate lesson extends the L22 composite model; it does not repeat its generic proof. Projected/truncated CG can exploit negative curvature in a bounded tangent subproblem, but step-quality and constraint safeguards still apply.

S-l1QP. Set $V(\mathbf d)=\sum_j|h_j+\mathbf J_{h,j}\mathbf d|+\sum_i\max(g_i+\mathbf J_{g,i}\mathbf d,0)$. Minimize $\mathbf a^{\mathsf T}\mathbf d+\frac12\mathbf d^{\mathsf T}\mathbf B\mathbf d+\rho_{\rm pen}V(\mathbf d)$ on the trust box. If physical bounds are hard, retain their intersection with that box; otherwise include their violations consistently in $V$. A bound-feasible current point makes zero belong to the hard domain. Softened constraint rows cannot make it empty. The nonempty hard domain is closed and bounded in finite-dimensional space, hence compact. By Weierstrass, a continuous model attains a global minimum there even if curvature is indefinite; uniqueness and cheap solution do not follow. Slack lifting makes the piecewise formulation a smooth QP with linear constraints, but it can still be nonconvex. Nocedal–Wright18.49–18.50pp549–550.

Here the penalty changes the computed step, not only its acceptance. The complete18.5update procedure onpp553–554first solves with the previous weight. If the selected step is exactly affine-feasible, retain the weight subject to the model-decrease safeguard. Otherwise solve a minimum-violation model on the same trust domain to obtain $\mathbf p_\infty$. If that reference reaches zero violation, increase the weight and resolve until the penalty-model step does so; if it cannot, require a prescribed fraction of its achievable violation reduction. Finally increase/resolve as needed so model reduction exceeds a fixed fraction of weight times the actual predicted violation reduction. If that reduction is zero, use the appropriate stationary/restoration branch; do not divide by it. With SLQP, the analogous LP stage uses zero quadratic curvature. These extra solves trade inner work against better penalty selection; they are not a universal efficiency guarantee.

Filters. A hard filter trust method solves the bounded affine QP when feasible and invokes restoration when its model cannot be solved/accepted. Historical acceptance and near-feasibility objective switching remain necessary. Biegler p. 148 separately describes an elastic filter line-search extension with projected-curvature, bounded-dual and sufficient-consistency assumptions; the next page gives the hard filter trust model. An elastic problem’s lifted feasibility does not guarantee that its residual is sufficiently small or that its numerical solve succeeds. Fletcher–Leyffer–Toint Theorem7contains method-specific restoration/KKT/regularity alternatives; see the complete L22 filter reference. Do not assign that theorem unchanged to every penalty/elastic/filter variant.

Local rates: what the hypotheses actually guarantee¶

A constructed sequence has two-step superlinearity and persistent one-step ratios.

Exact equality/SQP Newton steps have a quadratic local primal-dual rate under LICQ, SOSC, locally Lipschitz second derivatives, close primal/dual starts and accepted full steps. Stable active sets allow the inequality problem to use that local model; they do not supply every hypothesis by themselves. L21 owns the equality proof and L22 owns the Maratos safeguard.

Let $\mathbf P_k$ be an orthogonal projector onto the locally fixed active tangent and $\mathbf s_k=\mathbf x^{k+1}-\mathbf x^k$. A one-sided projected Dennis–Moré condition controls $\|\mathbf P_k(\mathbf B_k-\mathbf W^*)\mathbf s_k\|/\|\mathbf s_k\|\to0$. A two-sided condition controls only $\|\mathbf P_k(\mathbf B_k-\mathbf W^*)\mathbf P_k\mathbf s_k\|/\|\mathbf s_k\|\to0$. The missing normal/tangent coupling is a real difference: with $\mathbf P=\operatorname{diag}(0,1)$, $\mathbf B=2\mathbf I$, $\mathbf W=\left[\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\right]$ and $\mathbf s=(1,1)$, the two-sided numerator is0but the one-sided numerator is1.

Biegler6.24–6.27p142summarizes exact, one-sided superlinear and two-sided two-step conditions. Nocedal–Wright18.6p559gives ordinary full BFGS Q-superlinearity with strong full-Hessian positivity, close initial data and its regularity assumptions. Its printed parenthetical $r_k=s_k$ is a typo: ordinary BFGS uses $r_k=y_k$. The same page distinguishes damped BFGS R-superlinearity from one-step Q-superlinearity. Theorem 18.7assumes the sequence already converges, bounded approximations, its smoothness hypothesis and the two-sided projected condition to establish two-step superlinearity. The latter assumes $\mathbf x^k\to\mathbf x^*$ and nonzero errors $e_k=\|\mathbf x^k-\mathbf x^*\|$ and means: for every $\varepsilon>0$, there is $K$ such that $e_{k+2}/e_k<\varepsilon$ for all $k\ge K$. Equivalently,

It need not force every one-step ratio to zero. Let $a_0=.5$, $a_{j+1}=a_j^2$, $e_{2j}=a_j$ and $e_{2j+1}=a_j/2$. Both two-step ratios equal $a_j\to0$, while every even-to-odd ratio is .5. This is an error-sequence illustration, not a claimed trajectory of the prototype. Finite run tables can illustrate a rate but cannot prove asymptotic conditions. Nor is “quadratic if exact, two-step otherwise” a universal classification of every SQP method.

Interpreting the implementation laboratory¶

The original computational cells below are retained. They implement a dense active-set bound-QP solver, scalar-relaxed line-search SQP, damped BFGS and an optional reoptimized correction. They do not implement the trust-region/SLQP/sparse variants explained above.

• hessian='exact' uses the actual Lagrangian Hessian followed by full-space convexification. It does not necessarily take the unmodified exact Newton step. hessian='bfgs' uses the safeguarded update.

• The inner solver can fall back to least squares on singular KKT systems and has a finite iteration cap; its return does not itself prove an exact feasible QP optimum. It does not warm-start working sets across outer iterations.

• The scalar upper bound $\xi\le1$ and hybrid penalty logic are explicit prototype safeguards. A relaxation near1causes the original code’s stall status, not a proof of global infeasibility.

• The stopping test uses small primal step and equality violation rather than every original-model KKT residual. The original outer result does not return bound multipliers. Diagnostics identify any reconstructed final-bound multipliers separately from solver-returned equality values.

• The line search has a finite retry loop, and its current prototype path does not return a separate failure if that loop exhausts. The generic correction is equality-based and clips to bounds; the bounded example’s acceptance is empirical, and a general local-rate proof needs a bound-aware correction. The Maratos circle has no active physical bounds.

Fresh actual runs give near-roundoff original-point feasibility/stationarity for the circle and active-bound examples with both modes. On Biegler6.3from $(0,.1,0,0)$, BFGS reaches $(-3,2,6,0)$; the convexified-exact mode hits the100-iteration limit with violation $2.5148\times10^{-5}$ and stationarity defect .21557. From the all-zero start both keep the $\xi=1$ failure with violation1. Those outcomes are useful teaching evidence. They remain visible rather than being removed or “repaired” to make a comparison look successful. Decisions between full/reduced curvature, line search/trust methods and alternative solvers should use actual derivative cost, sparsity, residuals and failure traces rather than categorical method rankings.

Implementation laboratory: problem and KKT solve¶

The code uses $L=f+\mathbf v^{\mathsf T}\mathbf h$ with a row Jacobian and hard bounds $\mathbf x_L\le\mathbf x\le\mathbf x_U$. The bound QP and sign table define the model. lam stores the course equality multiplier $\mathbf v$; the objective-gradient argument g is not a vector of inequality functions. Nocedal–Wright equality multipliers reverse sign for the same equality function, while their nonnegative inequality multipliers retain their values when $\mathbf g=-\mathbf c_I$.

For equality-only subproblems, solve_kkt solves

A singular-system least-squares fallback is a numerical safeguard; inspect residuals before treating the returned values as a QP KKT solution.

Step 2: bounds, by an active set method¶

The QP chooses the active set, and re-chooses it every iteration -- that is what we bought by relaxing the known-active-set problem (6.9)--(6.10). The price is that a QP with inequalities is not a linear solve. We pay it with a textbook primal active set method: fix a working set of bounds, solve the resulting equality QP with solve_kkt, then add a blocking bound or release one whose multiplier has the wrong sign.

Step 3: scalar elastic relaxation¶

The original function below implements the guarded scalar QP, with $0\le\xi\le1$. A bound-feasible current point makes $(\mathbf d,\xi)=(\mathbf0,1)$ feasible. Positive full primal curvature and $M>0$ make its convex objective strictly convex. The upper bound on $\xi$ is an explicit addition to the book’s printed domain.

Positive relaxation can be preferred even when the original affine model is feasible. The original function’s short docstring should be read with these precise distinctions. Stack $\mathbf z=[\mathbf d;\xi]$ to obtain one more variable, Jacobian column $-\mathbf h$, and Hessian diagonal entry $M$.

Does the QP solver actually work?¶

Before trusting it inside SQP, check it against scipy.optimize.minimize with method='SLSQP' on 200 random, strictly convex, bound constrained QPs.

Hessian options¶

hessian="exact" evaluates the true Lagrangian Hessian, then applies a full-space convexification. Its resulting step can differ from exact SQP Newton. hessian="bfgs" uses damped BFGS and numerical skip safeguards. The constrained damping example explains why tangent SOSC can coexist with negative measured full-space curvature. The dense inner solver uses a positive model rather than exploiting indefinite trust-region curvature.

Line-search prototype¶

The full safeguarded algorithm and penalty comparison distinguish Biegler’s full QP multiplier update, Nocedal–Wright’s damped primal-dual update, and this prototype’s hybrid penalty logic. For its actual scalar step,

The code uses a multiplier/direct-slope maximum when relaxation is small and a direct-slope safeguard when relaxation is large. This avoids permanently importing a large elastic shadow price into the merit weight. It is a prototype choice, not a literal implementation of every textbook branch.

The original stopping test checks a small step and equality violation. The final diagnostics below also evaluate original-point stationarity and physical-bound complementarity. The finite backtracking loop lacks a separate exhaustion failure return; do not read its status string as verifying every globalization hypothesis.

Example 1: a problem with an analytic solution¶

By symmetry $\mathbf{x}^* = (-1, -1)^T$ and $f^* = -2$. Stationarity $\nabla f + \mathbf{J}_{\mathbf{h}}^T \lambda = \mathbf{0}$ reads $(1,1)^T + (-2,-2)^T \lambda = \mathbf{0}$, so $\lambda^* = 1/2$, and $\nabla^2_{\mathbf{xx}}\mathcal{L} = \lambda^* \cdot 2 \mathbf{I} = \mathbf{I} \succ 0$.

Observed error history¶

The plot shows $\|\mathbf x^k-\mathbf x^*\|$ for these particular runs. Local-rate conditions distinguish exact full-step Newton, Q-superlinear, R-superlinear and two-step guarantees. The convexified-exact mode may modify the true Hessian, and damped BFGS is not unconditionally classified as only two-step superlinear. A finite plot illustrates behavior; it does not prove an asymptotic theorem.

Example 2: an active bound¶

Without the bound the solution is $(1, 1)^T$ by symmetry. The bound $x_2 \leq 0.8$ cuts that off, so at the solution $x_2 = 0.8$ and $x_1 = \sqrt{2 - 0.64}$. The multipliers follow from stationarity. This exercises the active set machinery: the QP has to find that bound.

Example 3: Biegler’s Example 6.3 -- an inconsistent linearization¶

The feasible set is not empty -- $\mathbf{x}^* = [-3, 2, 6, 0]^T$ is the global solution -- but the linearization at $\mathbf{x}^0 = [0, 0.1, 0, 0]^T$ is inconsistent, and plain SQP stops dead.

Minimum feasible relaxation and the selected QP step¶

Adding the affine rows forces $d_2\ge4.95$, whereas its upper bound is 1.9. The relaxed rows require $1.98(1-\xi)\le.4(1.9)$, so the minimum feasible relaxation is $1-.76/1.98=.616161\ldots$. Compare that lower threshold with the actual selected value. This particular example has an inconsistent original affine model; positive selected relaxation alone does not prove inconsistency in general.

What the original algorithm does¶

From $(0,.1,0,0)$, the original BFGS prototype reaches $(-3,2,6,0)$ in the run below. The independent elimination argument, rather than the generic solver status, proves this example’s global minimum. A separate convexified-exact run below reaches its iteration cap with non-negligible stationarity error.

At the all-zero start, adding the relaxed rows gives $2(1-\xi)=-d_3-d_4\le0$, forcing $\xi=1$ and $d_3=d_4=0$. There are still feasible tangent values $0\le d_2\le2$; this example’s gradient and positive model choose the zero step. The original prototype retains its stall status. Restoration is a useful next mechanism, not a guarantee that every restored method succeeds.

Maratos computation¶

Globalization for Constrained Optimization supplies the complete rejection/correction proofs and acceptance safeguards. For this experiment, use $h=x_1^2+x_2^2-1$, $f=2h-x_1$, the solution $(1,0)$, multiplier $v^*=-3/2$, and Hessian $\mathbf W=\mathbf I$. Nocedal–Wright’s opposite multiplier sign comes from their $f-\lambda^{\mathsf T}c$ convention. The following cells keep the original model and algorithm computations.

Check the equality-QP step¶

At $\mathbf x^k=(\cos t,\sin t)$, compare the KKT solve with $\mathbf d=(\sin^2t,-\sin t\cos t)$. Its error is half the old error squared, but the full trial increases both objective and violation when $\sin t\ne0$.

The next table checks several positive angles near the solution. This behavior occurs arbitrarily close to it; a strict increase is not claimed at every nonzero angle (multiples of $\pi$ are exceptions).

Reoptimized replacement step¶

Use $\bar{\mathbf h}=\mathbf h(\mathbf x^k+\mathbf d)-\mathbf J\mathbf d$ as the new QP residual. The solve returns the whole replacement $\hat{\mathbf d}$, so test $\mathbf x^k+\hat{\mathbf d}$. The increment is $\hat{\mathbf d}-\mathbf d$. For linear constraints, this QP returns the original step and adding both would double it. The complete correction reference derives the constraint Taylor model and distinguishes this reoptimized variant from a minimum-normal correction: they coincide for this identity-Hessian example but can differ with Hessian coupling.

Shared step and acceptance illustration¶

The following archive and figure use the model and solver above. The handout and the comprehensive reference use this same exported figure.

With and without the correction¶

Compare the reported step lengths and residuals for the starts below. Full-step acceptance is a local, safeguarded property; a correction’s mere merit decrease does not establish Armijo acceptance at arbitrary starts.

Shared Lecture 23 computations¶

These figures use the actual solver functions above. The first separates the smallest feasible relaxation from the chosen weighted optimum; the second applies damping to a constrained step with positive tangent curvature. The final table checks original-point residuals, including preserved failure runs.

Damping at a constrained minimum¶

Use $f=-2x_1^2+\frac12(x_2-2)^2$, $h=x_1$, current $(-1,-1)$ and $\mathbf B=\operatorname{diag}(2,3)$. The exact full Hessian is indefinite, but the tangent direction $(0,1)$ has curvature 1. solve_kkt computes the step; damped_bfgs_update supplies the update. Scaling the same construction toward $(0,2)$ keeps $\mathbf s^{\mathsf T}\mathbf y<0$.

Original-point diagnostics across both Hessian modes¶

All runs below use eps_d=eps_h=1e-12 and max_iter=100. Equality multipliers come from the original result. The prototype does not return bound multipliers: the table reconstructs nonnegative values from the final gradient on bounds active within 1e-7, then checks stationarity and complementarity. This is a diagnostic reconstruction, not an inner-solver certificate. Missing bounds use zero multipliers. The plotting floor 1e-16 only makes zero residuals visible on the log scale.

Algorithms beyond the prototype¶

The native guide above includes the trust-region families, SLQP, reduced-space bounds, and sparse factorization updates. The code uses dense systems and bound-form inequalities, with no outer working-set warm start or feasibility restoration implementation. The complete globalization reference describes restoration and filter safeguards. Reaching a useful nonlinear point from a poor start remains an algorithmic task; these omissions are not guarantees of failure or success.