Continuous Optimisation HT26, Inequality constrained optimisation problems

Problem (iCP):

where $f : \mathbb R^n \to \mathbb R$, $c = (c _ 1, \ldots, c _ p) : \mathbb R^n \to \mathbb R^p$ are smooth. The aim is to find KKT points (or local minimisers of this problem.

Strictly feasible set:

i.e. the interior of the feasible region.

Assumption: $\Omega^0 \ne \emptyset$, i.e. the ∆scq (Slater) condition.

For each $\mu > 0$, the logarithmic barrier function is:

The associated logarithmic barrier subproblem is

This is essentially unconstrained. The constraint $c(x) > 0$ is enforced implicitly by the barrier term, since $c _ i(x) \to 0^+ \implies -\log c _ i(x) \to +\infty$.

Role of $\mu$:

• $\mu \to \infty$: the barrier term dominates, so a minimiser $x(\mu)$ is pushed deep into $\Omega^o$, far from the boundary. Strict feasibility is well-enforced.
• $\mu \to 0$: the barrier term vanishes pointwise on $\Omega^o$, so $x(\mu)$ approaches the constrained minimiser. This is where ill-conditioning can occur.

Gradient:

where $J(x)$ is the constraint Jacobian (∆constraint-jacobian-matrices), and $c(x)^{-1} = (1/c _ 1(x), \ldots, 1/c _ p(x))$.

Stationarity for $\text{iCP} _ \mu$: We have that $x(\mu)$ is a minimiser iff $\nabla f _ \mu(x(\mu)) = 0$, or equivalently

with $\lambda _ i(\mu) > 0$ on $\Omega^o$.

KKT for $\text{iCP}$: We require

and $\lambda^\ast \ge 0$, $\lambda _ i^\ast c _ i(x^\ast) = 0$.

Comparison: These two expressions have the same algebraic shape:

This identifies $\lambda _ i(\mu) := \mu / c _ i(x(\mu))$ as a natural Lagrange multiplier estimate, and motivates the central question: under what conditions do $x(\mu) \to x^\ast$ and $\lambda _ i(\mu) \to \lambda _ i^\ast$ as $\mu \to 0$?

The central path is the set

where $x(\mu)$ is a (global) minimiser of $f _ \mu$ on $\Omega^o$. As $\mu$ decreases, $x(\mu)$ traces a curve deep in $\Omega^o$ toward the boundary of $\Omega$, and under suitable conditions (see ∆central-path-existence) converges to a constrained minimiser of $\text{iCP}$ as $\mu \to 0$.

Suppose $\Omega^o \ne \emptyset$, and $x^\ast$ is a local minimiser of $\text{iCP}$ with KKT multipliers $\lambda^\ast$, satisfying:

1. Strict complementarity: $\lambda _ i^\ast > 0$ whenever $c _ i(x^\ast) = 0$.
2. ∆licq: $\nabla c _ i(x^\ast), i \in \mathcal A := \{i : c _ i(x^\ast) = 0\}$, are linearly independent.
3. Second-order sufficient condition: there exists $\alpha > 0$ such that

where $\mathcal L$ is the ∆lagrangian-of-constrained-optimisation-problem.

Then there is a unique $\mathcal C^1$ path $\mu \mapsto x(\mu)$ of (global) minimisers of $f _ \mu$ defined for all $\mu > 0$ sufficiently small, and $x(\mu) \to x^\ast$ as $\mu \to 0$.

Convex:

Nonconvex:

Difference:

In the nonconvex case, there are two local minimisers. Hence there are two central paths starting within the feasible region (for $\mu$ sufficiently small), both of which converge to each local minimiser.

Inputs: initial $\mu^0 > 0$, starting iterate $x^0 \in \Omega^o$.

Algorithm: Let $k = 0$. Until convergence:

• Choose $0 < \mu^{k+1} < \mu^k$ (commonly $\mu^{k+1} = 0.1 \mu^k$ or $\mu^{k+1} = (\mu^k)^2$).
• Choose a starting point $x_0^k \in \Omega^o$ for the inner minimisation (commonly $x_0^k := x^k$).
• Inner minimisation: from $x_0^k$, use an unconstrained minimiser to find an approximate minimiser $x^{k+1}$ of $f_{\mu^{k+1}}$.
• $k := k+1$

Convergence requirement: $\mu^k \to 0$ as $k \to \infty$.

Inner-minimisation algorithms: $f_\mu$ has log singularities at the boundary of $\Omega$, so we may use:

• Linesearch methods: use special linesearch to cope with the singularity of the log.
• Trust-region methods: shape the trust region to cope with the log contours, reject any candidate with $c(x^k + s^k) \not> 0$.

Suppose:

• $f, c \in \mathcal C^1$
• $\|\nabla f _ {\mu^k}(x^k) \| \le \epsilon^k$, where $\epsilon^k \to 0$ as $k \to \infty$
• $\mu^k \to 0$ as $k \to \infty$
• $x^k \to x^\ast$ where $\nabla c _ i(x^\ast), i \in \mathcal A := \{i : c _ i(x^\ast) = 0\}$, are linear independent (∆licq at $x^\ast$)

Then $x^\ast$ is a KKT point of $\text{iCP}$, and $\lambda^k \to \lambda^\ast$, where $\lambda^\ast$ is the vector of Lagrange multipliers at $x^\ast$.

(compare to ∆quadratic-penalty-method-convergence, ∆aug-lagrangian-global-convergence).

By ∆licq, $J _ \mathcal A(x^\ast)$ has full row rank, so the pseudo-inverse $J _ {\mathcal A}(x^\ast)^+ = (J _ \mathcal A J _ \mathcal A^\top)^{-1} J _ {\mathcal A}$ is well-defined. By continuity, $J _ {\mathcal A}(x^k)^+$ exists, is bounded, and converges to $J _ {\mathcal A}(x^\ast)^+$ for $k$ large.

Let $\mathcal I = \{1, \ldots, p\} \setminus \mathcal A$. From $c(x^k) > 0$ and $x^k \to x^\ast$: $c(x^\ast) \ge 0$ with $c _ {\mathcal A}(x^\ast) = 0$ and $c _ {\mathcal I}(x^\ast) > 0$, so $x^\ast$ is feasible (primal feasibility).

Define $\lambda^\ast = (\lambda _ {\mathcal A}^\ast, \lambda^\ast _ {\mathcal I})$ via

Complementary slackness $\lambda _ i^\ast c _ i(x^\ast) = 0$ is immediate. $\lambda _ {\mathcal I}^\ast = 0$ means inactive $c _ {\mathcal A}(x^\ast) = 0$ kills active ones.

It remains to show that:

1. $\lambda^k _ {\mathcal I} \to 0$
2. $\lambda^k _ {\mathcal A} \to \lambda^\ast _ {\mathcal A} \ge 0$
3. Lagrangian stationarity $\nabla f(x^\ast) = J _ {\mathcal A}(x^\ast)^\top \lambda^\ast _ {\mathcal A}$.

For (1), note that for $i \in \mathcal I$, $\lambda^k _ i = \mu^k / c _ i(x^k) \to / c _ i(x^\ast) = 0$, since $\mu^k \to 0$ and $c _ i(x^\ast) > 0$.

For (2), write $r^k := \nabla f(x^k) - J _ {\mathcal A}(x^k)^\top \lambda^k _ {\mathcal A}$ for the active-block residual. The multiplier estimate (∆barrier-multiplier-estimates) gives $\nabla f _ {\mu^k}(x^k) = \nabla f(x^k) - J(x^k)^\top \lambda^k$, and splitting $J^\top \lambda^k = J _ {\mathcal A}^\top \lambda^k _ {\mathcal A} + J _ {\mathcal I}^\top \lambda^k _ {\mathcal I}$ rearranges this to $r^k = \nabla f _ {\mu^k}(x^k) + J _ {\mathcal I}(x^k)^\top \lambda^k _ {\mathcal I}$. Hence

where:

• $(\star 1)$: substituting the identity $r^k = \nabla f _ {\mu^k}(x^k) + J _ {\mathcal I}(x^k)^\top \lambda^k _ {\mathcal I}$ derived above.
• $(\star 2)$: triangle inequality, then submultiplicativity $\|J _ {\mathcal I}(x^k)^\top \lambda^k _ {\mathcal I}\| \le \|J _ {\mathcal I}(x^k)\| \cdot \|\lambda^k _ {\mathcal I}\|$.
• $(\star 3)$: the approximate-stationarity hypothesis $\|\nabla f _ {\mu^k}(x^k)\| \le \epsilon^k$.
• $(\star 4)$: $\epsilon^k \to 0$; $\|J _ {\mathcal I}(x^k)\|$ is bounded ($J _ {\mathcal I}$ continuous along $x^k \to x^\ast$); and $\|\lambda^k _ {\mathcal I}\| \to 0$ by (1).

LICQ makes $J _ {\mathcal A}(x^k)$ full row rank, so $J _ {\mathcal A}(x^k)^+ J _ {\mathcal A}(x^k)^\top = (J _ {\mathcal A} J _ {\mathcal A}^\top)^{-1}(J _ {\mathcal A} J _ {\mathcal A}^\top) = I$; left-multiplying $r^k = \nabla f(x^k) - J _ {\mathcal A}(x^k)^\top \lambda^k _ {\mathcal A}$ by $J _ {\mathcal A}(x^k)^+$ and rearranging gives $\lambda^k _ {\mathcal A} = J _ {\mathcal A}(x^k)^+ \nabla f(x^k) - J _ {\mathcal A}(x^k)^+ r^k$. Hence

where:

• $(\star 5)$: substituting $\lambda^k _ {\mathcal A} = J _ {\mathcal A}(x^k)^+ \nabla f(x^k) - J _ {\mathcal A}(x^k)^+ r^k$ and grouping the two limit-bearing pieces.
• $(\star 6)$: triangle inequality and submultiplicativity $\|J _ {\mathcal A}(x^k)^+ r^k\| \le \|J _ {\mathcal A}(x^k)^+\| \cdot \|r^k\|$.
• $(\star 7)$: the first term $\to 0$ because $J _ {\mathcal A}(x^k)^+ \nabla f(x^k) \to J _ {\mathcal A}(x^\ast)^+ \nabla f(x^\ast) = \lambda^\ast _ {\mathcal A}$ (continuity of $J _ {\mathcal A}^+$ and $\nabla f$, $x^k \to x^\ast$, and the definition of $\lambda^\ast _ {\mathcal A}$); the second term $\to 0$ because $\|J _ {\mathcal A}(x^k)^+\|$ is bounded and $\|r^k\| \to 0$ by $(\star 4)$.

Therefore $\lambda^k _ {\mathcal A} \to \lambda^\ast _ {\mathcal A}$, and since $\lambda^k > 0$ for all $k$, $\lambda^\ast _ {\mathcal A} \ge 0$ (dual feasibility).

For (3), pass to the limit in $r^k = \nabla f(x^k) - J _ {\mathcal A}(x^k)^\top \lambda^k _ {\mathcal A}$. By continuity and (2), it converges to $\nabla f(x^\ast) - J _ {\mathcal A}(x^\ast)^\top \lambda^\ast _ {\mathcal A}$. But then $(\star 4)$ says it also converges to $0$. Hence

All four ∆kkt-conditions hold at $(x^\ast, \lambda^\ast)$, so $x^\ast$ is KKT, and $\lambda^k \to \lambda^\ast$ overall.

Gradient:

So minimising $f _ \mu$ in $x$ is equivalent to seeking Lagrangian stationarity $(x, \lambda(x))$.

Hessian: (writing $C(x) := \text{diag}(c _ 1(x), \ldots, c _ p(x))$):

The “extra” term $\mu J(x)^\top C(x)^{-2} J(x)$ comes from differentiating the multiplier estimate $\lambda(x) = \mu c(x)^{-1}$: each entry $\lambda _ i(x) = \mu/c _ i(x)$ contributes $-(\mu/c _ i(x)^2)\,\nabla c _ i(x) \nabla c _ i(x)^\top$ when re-differentiated. For $i \in \mathcal A$ (active at $x^\ast$), $c _ i(x(\mu)) = O(\mu)$ along the central path, so $\mu/c _ i(x)^2 = O(1/\mu) \to \infty$ as $\mu \to 0$ — this is the mechanism behind ∆barrier-hessian-ill-conditioning.

Along the central path with $x(\mu) \to x^\ast$ and $\lambda _ i(\mu) = \mu/c _ i(x(\mu))$ staying bounded (i.e. $c _ i(x(\mu)) \to c _ i(x^\ast)$ at rate $O(\mu)$ for active constraints):

• $i \in \mathcal A$ (active, $c _ i(x^\ast) = 0$): $\mu / c _ i(x(\mu))^2 \to \infty$
• $i \in \mathcal I$ (inactive, $c _ i(x^\ast) > 0$): $\mu / c _ i(x(\mu))^2 \to 0$.

So as $\mu \to 0$ (along with a Fact from Gould), $\nabla^2 f _ {\mu}(x(\mu))$ splits into:

• $p _ a := \vert \mathcal A \vert $ eigenvalues of $\nabla^2 f _ {\mu^k}(x^k)$ tend to infinity as $k \to \infty$.
• $n - p _ a$ eigenvalues that are $O(1)$ (from $\nabla^2 _ {xx} \mathcal L$) (inactive contributions to the extra term vanish)
• The condition number of $\nabla^2 f _ {\mu^k}(x^k)$ is $O(1/\mu^k)$.

This blows up as $k \to \infty$, and we may not able to compute $x^k$ accurately.

Even though $x^k$ is an approximate minimiser of $f _ {\mu^k}$ (and hence interior), once we shift to the new barrier $f _ {\mu^{k+1}}$, the full Newton step $s^k$ from $x^k _ 0 = x^k$ is asymptotically infeasible.:

Any non-trivial decrease in $\mu$ (anything more aggressive than halving) drives the Newton iterate out of $\Omega^o$. So either:

• We must take a damped step (small $\alpha^k$), which loses the local quadratic convergence of Newton, or
• Use only small decreases $\mu^{k+1} \ge 0.5 \mu^k$, which means we have very slow outer convergence.

• The Hessian of $f _ {\mu^k}$ becomes ill-conditioned as $\mu^k \to 0$ (condition number $O(1/\mu^k)$), so e.g. Newton’s method will get worse and worse (∆barrier-hessian-ill-conditioning).
• It can be shown that the natural starting point $x^k _ 0 := x^k$ is a poor choice, since in order for iterates to remain feasible we must either dampen the Newton step (losing quadratic convergence), or only take small decreases in $\mu^k$, which means the outer loop converges very slowly (∆barrier-poor-starting-points).

The technique to address these issues is ∆primal-dual-newton-system, which treat the multiplier estimates as explicit variables.

First-order necessary conditions for $\text{iCP} _ \mu$ (which is essentially unconstrained, since $f _ \mu \to \infty$ at the boundary) give $\nabla f _ {\mu} (x(\mu)) = 0$:

With $\lambda(\mu) := \mu c(x(\mu))^{-1}$, this is Lagrangian stationarity:

Multiplying $\lambda _ i(\mu) = \mu / c _ i(x(\mu))$ by $c _ i(x(\mu))$ gives perturbed complementarity:

Strict feasibility $x(\mu) \in \Omega^o$ gives $\boxed{c > 0, \lambda > 0}$.

So $(x(\mu), \lambda(\mu))$ satisfies $\text{OPT} _ \mu$:

Contrast this with KKT for $\text{iCP}$:

There are two differences: complementarity is $c _ i \lambda _ i = \mu$ vs $0$, and the inequalities are strict vs non-strict. As $\mu \to 0$, $\text{OPT} _ \mu$ recovers KKT.

Why this matters: KKT complementarity $c _ i \lambda _ i$ is a combinatorial statement: for each $i$, either $c _ i = 0$ (active) or $\lambda _ i = 0$ (inactive), and the algorithm has to discover which. $\text{OPT} _ \mu$ loosens this to a smooth equation $c _ i \lambda _ i = \mu$, so that every constraint may be “slightly active”, and Newton’s method now applies to the joint $(x, \lambda)$ unknowns (in contrast to ∆barrier-multiplier-estimates, where $\lambda$ was derived from $x$). As $\mu \to 0$, the active set emerges automatically: active constraints have $c _ i \to 0$ with $\lambda _ i$ bounded, inactive ones $\lambda _ i \to 0$ with $c _ i$ bounded. In this view, $\text{OPT} _ \mu$ is a smooth perturbation of KKT, and primal-dual methods follow the path of solutions to $\mu = 0$ (see ∆primal-dual-newton-system).

Write the system as $g(x, \lambda) = 0$ with $g = (g _ 1; g _ 2)$, where

with $C(x) := \text{diag}(c _ 1(x), \ldots, c _ p(x))$ and $e = (1, \ldots, 1)^\top$. Linearising, we obtain

• $\partial g _ 1 / \partial x = \nabla^2 f(x) - \sum _ i \lambda _ i \nabla^2 c _ i(x) = \nabla^2 _ {xx} \mathcal L(x, \lambda)$ (cf. ∆barrier-lagrangian-connection)
• $\partial g _ 1 / \partial \lambda = -J(x)^\top$

and

• $\partial g _ 2 / \partial x = \Lambda J(x)$ where $\Lambda := \text{diag}(\lambda)$ (the $i$th component $c _ i(x) \lambda _ i$ differentiates to row $\lambda _ i \nabla c _ i(x)^\top)$
• $\partial g _ 2 / \partial \Lambda = C(x)$

Then Newton’s method solves for

Reduce by eliminating $d\lambda$: from the bottom block,

Substituting into the top block, the $J^\top \lambda$ terms cancel and we get

Why this matters: The RHS matches the primal Newton step on $f _ \mu$ (∆barrier-lagrangian-connection). The LHS differs by using explicit $\lambda$ where primal methods are forced to use the implicit $\lambda(x) := \mu c(x)^{-1}$. So the multipliers get their own Newton update $\text d \lambda$ rather than being recomputed by division by vanishing $c _ i$ each iteration. This decouples the multiplier dynamics from the singularity $\lambda(x) = \mu / c(x)$ that creates the ∆barrier-hessian-ill-conditioning and ∆barrier-poor-starting-points issues. The multipliers are now full state variables that Newton can update directly, rather than quantities derived from $x$.

Inputs: initial $\mu^0 > 0$ and $(x^0, \lambda^0)$ with $c(x^0) > 0, \lambda^0 > 0$ (e.g. $\lambda _ i^0 := \mu^0 / c _ i(x^0))$, $\mu$-shrinkage rule (often $\mu^{k+1} := (\mu^k)^2)$.

Algorithm: Until $\mu^k$ is small and KKT residual below tolerance:

• Inner Newton on $\text{OPT} _ {\mu^k}$, warm-started at $(x^k, \lambda^k)$, repeat a small number of times:
  • Compute $(\text dx, \text d\lambda)$ from ∆primal-dual-newton-system at the current $(x, \lambda)$ with parameter $\mu^k$.
  • Damped step $(x, \lambda) \leftarrow (x + \alpha \text dx , \lambda + \alpha \text d\lambda)$ with $\alpha \in (0, 1]$ chosen by linesearch / trust-region to:
    • Maintain $c(x) > 0$ and $\lambda > 0$ (say strictly inside the cone in both primal and dual);
    • Descent the $\text{OPT} _ {\mu^k}$ residual norm.
  • Repeat until approximate stationarity at $\text{OPT} _ {\mu^k}$.
• Set $(x^{k+1}, \lambda^{k+1})$ to the converged inner iterate.
• Shrink $\mu^{k+1} := \tau(\mu^k)$, commonly $\mu^{k+1} = (\mu^k)^2$.
• $k := k+1$.

Simplex method: walk along the vertices of the feasible polytope, at each pivot stepping to an adjacent vertex with better objective.

• Worst-case: exponential in problem dimension
• Average-cost: linear polynomial in practice
• Per-iteration cost: cheap pivot operation on the simplex tableau

Primal-dual IPM: cuts through the interior along the central path.

• Worst-case: polynomial in dimension and input length
• Average-cost: $O(\log (\text{LP dimension}))$ outer iterations until convergence
• Per-iteration cost: $O(N^3)$ for the primal-dual Newton solve

In practice:

• For small / medium LPs, simplex is often faster
• For very large linear programs, IPMs are often faster

• Setup: LICQ ⟹ $J _ {\mathcal A}(x^\ast)$ has full row rank, so the pseudo-inverse $J _ {\mathcal A}(x^\ast)^+ = (J _ {\mathcal A} J _ {\mathcal A}^\top)^{-1} J _ {\mathcal A}$ exists; by continuity it exists, is bounded, and is continuous at $x^k$ for large $k$. Split indices into active $\mathcal A := \{i : c _ i(x^\ast) = 0\}$ and inactive $\mathcal I$.
• Define the candidate $\lambda^\ast = (\lambda^\ast _ {\mathcal A}, 0)$ via $\lambda^\ast _ {\mathcal A} := J _ {\mathcal A}(x^\ast)^+ \nabla f(x^\ast)$. Complementary slackness $\lambda^\ast _ i c _ i(x^\ast) = 0$ is then immediate.
• Step 1 (inactive multipliers vanish): $\lambda^k _ i = \mu^k / c _ i(x^k) \to 0 / c _ i(x^\ast) = 0$ for $i \in \mathcal I$, using $\mu^k \to 0$ and $c _ i(x^\ast) > 0$.
• Step 2 (active multipliers converge): starting from $\nabla f _ {\mu^k}(x^k) = \nabla f(x^k) - J(x^k)^\top \lambda^k$, isolate the active part and bound $\|\nabla f(x^k) - J _ {\mathcal A}(x^k)^\top \lambda^k _ {\mathcal A}\| \le \epsilon^k + \|J _ {\mathcal I}(x^k)\| \|\lambda^k _ {\mathcal I}\| \to 0$. Apply $J _ {\mathcal A}(x^k)^+$ to deduce $\lambda^k _ {\mathcal A} \to J _ {\mathcal A}(x^\ast)^+ \nabla f(x^\ast) = \lambda^\ast _ {\mathcal A}$. Since $\lambda^k > 0$, $\lambda^\ast _ {\mathcal A} \ge 0$.
• Step 3 (stationarity): passing to the limit in $\nabla f(x^k) - J _ {\mathcal A}(x^k)^\top \lambda^k _ {\mathcal A}$ gives both $\nabla f(x^\ast) - J _ {\mathcal A}(x^\ast)^\top \lambda^\ast _ {\mathcal A}$ (by continuity) and $0$ (by Step 2’s bound), hence $\nabla f(x^\ast) = J _ {\mathcal A}(x^\ast)^\top \lambda^\ast _ {\mathcal A}$ — Lagrangian stationarity.

All four KKT pieces (primal feasibility, stationarity, dual feasibility, complementary slackness) hold at $(x^\ast, \lambda^\ast)$.

• LICQ ⟺ $J _ {\mathcal A}(x^\ast)$ has full row rank ⟺ $J _ {\mathcal A} J _ {\mathcal A}^\top$ is invertible, so the pseudo-inverse $J _ {\mathcal A}(x^\ast)^+ := (J _ {\mathcal A} J _ {\mathcal A}^\top)^{-1} J _ {\mathcal A}$ exists, is bounded, and (by continuity) is continuous at $x^k$ for large $k$. Without this, the pseudo-inverse blows up or is undefined.
• The pseudo-inverse is used twice: to define $\lambda^\ast _ {\mathcal A} := J _ {\mathcal A}(x^\ast)^+ \nabla f(x^\ast)$, and in Step 2 to convert the active-block residual bound into a bound on $\lambda^k _ {\mathcal A}$ itself (“apply $J _ {\mathcal A}(x^k)^+$ to both sides”).
• Without LICQ: the multiplier candidate $\lambda^\ast$ may not exist, or be non-unique; the limit $\lambda^k$ can blow up or oscillate; convergence to a KKT point cannot be proved this way.

The same role is played by LICQ in the analogous penalty-method proof (∆quadratic-penalty-method-convergence-proof) and aug-Lagrangian proof (∆aug-lagrangian-global-convergence-proof).