Optimization

Motivation

In the contextual bandit problem a learner repeatedly observes a context, picks one of $K$ actions, and sees the reward of that action only. Against a finite class $\Pi$ of policies, statistically optimal regret of order $\sqrt{KT\ln|\Pi|}$ has been known since EXP4 (Auer et al., 2002), but EXP4 maintains a weight per policy and costs $\Omega(|\Pi|)$ time per round. For the large policy classes used in practice (linear classifiers, trees), that is prohibitive.

The oracle-efficient line of work accesses $\Pi$ only through a cost-sensitive classification oracle (an arg max oracle, AMO). The RandomizedUCB algorithm of Dudík et al. (2011) obtains optimal regret with polynomially many oracle calls by solving a convex program in each round, but the number of calls is large. Agarwal, Hsu, Kale, Langford, Li and Schapire (2014) replace that solver by a coordinate descent method whose number of iterations, and hence of oracle calls, is bounded independently of $|\Pi|$. Their algorithm, ILOVETOCONBANDITS, and its practical variant are now standard references for oracle-based exploration.

This mission formalizes the optimization half of that paper: Algorithm 2 solves the per-epoch problem (OP) after at most $4\ln(1/(K\mu))/\mu$ coordinate steps.

Setting

Let $A=\{0,\dots,K-1\}$ be the actions, $X$ any set of contexts, and $\Pi\subseteq A^X$ a finite nonempty set of policies. A history $H_t$ is a sequence of $t\ge1$ records $(x_i,a_i,r_i(a_i),p_i(a_i))$ with $r_i(a_i)\in[0,1]$ the observed reward and $p_i(a_i)\in(0,1]$ the probability with which $a_i$ was chosen. Write $\widehat{\mathbb E}_{x\sim H_t}[f(x)]=\frac1t\sum_i f(x_i)$.

The inverse propensity scoring estimate (Eq. (1)) is

$$\widehat{\mathcal R}_t(\pi)=\frac1t\sum_{i=1}^t\frac{r_i(a_i)\,\mathbb 1\{\pi(x_i)=a_i\}}{p_i(a_i)},$$

the estimated regret is $\widehat{\mathrm{Reg}}_t(\pi)=\max_{\pi'\in\Pi}\widehat{\mathcal R}_t(\pi')-\widehat{\mathcal R}_t(\pi)$, and for a minimum probability $\mu$ one sets $b_\pi=\widehat{\mathrm{Reg}}_t(\pi)/(\psi\mu)$ with $\psi=100$.

Weights are vectors $Q\in\mathbb R^\Pi$; $\Delta^\Pi$ is the set of nonnegative $Q$ with $\sum_\pi Q(\pi)\le1$. The smoothed projection of $Q$ is

$$Q^\mu(a\mid x)=(1-K\mu)\sum_{\pi:\ \pi(x)=a}Q(\pi)+\mu .$$

The optimization problem (OP) asks for $Q\in\Delta^\Pi$ with

$$\sum_{\pi\in\Pi}Q(\pi)b_\pi\le2K\quad(2),\qquad \widehat{\mathbb E}_{x\sim H_t}\Bigl[\frac1{Q^\mu(\pi(x)\mid x)}\Bigr]\le2K+b_\pi\ \ \forall\pi\in\Pi\quad(3).$$

Algorithm 2 starts from $Q_{\mathrm{init}}$ and loops. With $V_\pi(Q)=\widehat{\mathbb E}[1/Q^\mu(\pi(x)\mid x)]$, $S_\pi(Q)=\widehat{\mathbb E}[1/Q^\mu(\pi(x)\mid x)^2]$ and $D_\pi(Q)=V_\pi(Q)-(2K+b_\pi)$: if $\sum_\pi Q(\pi)(2K+b_\pi)>2K$ it rescales $Q$ by $c=2K/\sum_\pi Q(\pi)(2K+b_\pi)$ (Eq. (4)); then, if some $\pi$ has $D_\pi(Q)>0$, it adds

$$\alpha_\pi(Q)=\frac{V_\pi(Q)+D_\pi(Q)}{2(1-K\mu)S_\pi(Q)}$$

to $Q(\pi)$ (Step 8) and repeats; otherwise it halts and outputs $Q$.

The analysis uses the potential (Eq. (6)), with $\tau=t$ and $\mathcal U_A$ uniform on $A$,

$$\Phi_m(Q)=\tau\mu\left(\frac{\widehat{\mathbb E}_x[\mathrm{RE}(\mathcal U_A\,\|\,Q^\mu(\cdot\mid x))]}{1-K\mu}+\frac{\sum_\pi Q(\pi)b_\pi}{2K}\right),\qquad \mathrm{RE}(p\|q)=\sum_a\bigl(p_a\ln\tfrac{p_a}{q_a}+q_a-p_a\bigr).$$

Formalization targets

Goal: Theorem 3 (p. 6)

For $0<\mu\le1/(2K)$, Algorithm 2 with $Q_{\mathrm{init}}=\mathbf 0$ satisfies: every run executes Step 8 at most

$$\frac{4\ln(1/(K\mu))}{\mu}$$

times, whatever policy each Step 8 chooses among those with $D_\pi>0$; and when it halts, its output solves (OP). The bound depends on $K\mu$ only, not on $|\Pi|$ or $t$.

Milestones

1. Lemma 5 (p. 10). If Algorithm 2 halts and outputs $Q$, then $Q$ satisfies (2), (3) and $\sum_\pi Q(\pi)\le1$.
2. Lemma 6 (p. 10). If $\sum_\pi Q(\pi)(2K+b_\pi)>2K$ and $c$ is as in Eq. (4), then $\Phi_m(cQ)\le\Phi_m(Q)$.
3. Lemma 7 (p. 10). If $D_\pi(Q)>0$ and $Q'$ adds $\alpha_\pi(Q)$ to $Q(\pi)$, then

$$\Phi_m(Q)-\Phi_m(Q')\ge\frac{\tau\mu^2}{4(1-K\mu)}.$$

Significance

The result. Theorem 3 is what makes ILOVETOCONBANDITS computationally efficient: each call of Algorithm 2 is implemented with one AMO call per iteration (Lemma 1 of the paper), so the oracle complexity of an epoch is $O(\ln(1/(K\mu))/\mu)$. Combined with the epoch schedule and warm start, this gives the paper's total of $\tilde O(\sqrt{KT/\ln(|\Pi|/\delta)})$ oracle calls over $T$ rounds. Theorem 3 also gives a constructive proof that (OP) is feasible for every history, which the regret analysis (a separate mission in this series) assumes.

Formalizing it. The result is proved in the paper, with complete proofs of Lemmas 5–7 in Appendix D. No machine-checked version is known. The formalization would give a checked termination bound for a coordinate descent method on a non-smooth feasibility problem, with a fully explicit constant, and a verified definition of the unnormalized relative entropy potential that is reusable for other smoothed-projection analyses (e.g. RandomizedUCB-type convex programs).

Difficulty

Termination cannot be read off the constraints. Step 8 raises one weight and can push the total weight above 1, after which Step 5 shrinks every coordinate, so no constraint and no single weight moves monotonically along a run. The number of policies that violate (3) can also go up after a step. Bounding the number of iterations therefore needs a global quantity that tracks progress through both kinds of step. The rescaling step is the harder of the two: it lowers every $Q^\mu(a\mid x)$ at once, which pushes the relative-entropy term the wrong way, and it must be offset by the drop in the regret term. Knowing that (OP) is feasible, or that some convex function has a minimizer, bounds nothing about how many steps a particular method takes; that is the obvious approach, and it gives no count.

Formalization scope

• Actions are Fin K with $K\ge1$; contexts form an arbitrary type (no measure is needed: Theorem 3 is deterministic). $\Pi$ is a nonempty Finset (X → Fin K); weights are real functions on its subtype.
• Histories are indexed by Fin t with $t\ge1$ (0-based indices). The paper allows $p_i(a_i)\in[0,1]$; the formalization requires $p_i(a_i)\in(0,1]$, since Eq. (1) divides by it.
• $\widehat{\mathrm{Reg}}_t(\pi)$ is written as $\max_{\pi'}\widehat{\mathcal R}_t(\pi')-\widehat{\mathcal R}_t(\pi)$, which equals $\widehat{\mathcal R}_t(\pi_t)-\widehat{\mathcal R}_t(\pi)$ for any maximizer $\pi_t$; $\psi=100$ is hard-wired in $b_\pi$.
• $Q^\mu$, $V_\pi$, $S_\pi$, (OP) and $\Phi_m$ all use the smoothed projection of the unnormalized weights; there is no default policy in this mission.
• $\mu$ ranges over $(0,1/(2K)]$, the range of $\mu_m$ in Algorithm 1 that the printed theorem refers to. $\tau$ in $\Phi_m$ is the history length $t$.
• Algorithm 2 is encoded relationally. A run of length $n$ from $Q_{\mathrm{init}}$ is a sequence $Q^{(0)}=Q_{\mathrm{init}},\dots,Q^{(n)}$ in which each $Q^{(k+1)}$ is Step 8, for some policy with $D_\pi>0$, applied to the rescaled $Q^{(k)}$. It halts at $Q^{(n)}$ when no policy has $D_\pi>0$ after rescaling, and it then outputs the rescaled $Q^{(n)}$. "Iterations" means executions of Step 8. The last pass, which halts at Step 10, is not counted: the paper's proof bounds "the number of times Step 8 is executed". The bound is compared in $\mathbb R$, without rounding.
• Lemmas 5–7 are stated for nonnegative weight vectors without a bound on their sum, because Algorithm 2 rescales vectors whose sum may exceed 1. Lemma 7's "$\alpha=\alpha_\pi(Q)>0$" is part of its conclusion.
• A trivializing formalization is ruled out: the goal quantifies over every run from $\mathbf 0$ and every choice in Step 8, not over some run, and the potential, $b_\pi$ and $D_\pi$ are computed from the history rather than taken as free parameters.
• The auxiliary facts $\Phi_m\ge0$ and $\Phi_m(\mathbf 0)\le\tau\mu\ln(1/(K\mu))/(1-K\mu)$ are inline claims in the paper and are not stated separately; contributions stating and proving them are welcome, as are general lemmas on the unnormalized relative entropy.
• Out of scope: the regret bound (Theorem 2) and the probabilistic model (mission I of this series), the AMO implementation (Lemma 1), warm start and epoch-level oracle counts (Lemmas 2, 3, 8), and the support lower bound (Theorem 4).

Selected references

• A. Agarwal, D. Hsu, S. Kale, J. Langford, L. Li, R. E. Schapire, Taming the Monster: A Fast and Simple Algorithm for Contextual Bandits, ICML 2014; arXiv:1402.0555v2. https://arxiv.org/abs/1402.0555
• M. Dudík, D. Hsu, S. Kale, N. Karampatziakis, J. Langford, L. Reyzin, T. Zhang, Efficient Optimal Learning for Contextual Bandits, UAI 2011. https://arxiv.org/abs/1106.2369
• P. Auer, N. Cesa-Bianchi, Y. Freund, R. E. Schapire, The Nonstochastic Multiarmed Bandit Problem, SIAM J. Comput. 32(1), 2002. https://doi.org/10.1137/S0097539701398375

Motivation

A polyhedral cone can be described in two ways: as the set of nonnegative combinations of finitely many vectors (a finitely generated cone), or as the intersection of finitely many closed half-spaces through the origin. That the two descriptions give the same class of sets is the Minkowski–Weyl theorem. It is the structural basis of linear programming: the simplex method, LP duality, Farkas' lemma, and the vertex/facet description of polytopes used throughout combinatorial optimization all rest on it.

Hermann Weyl's 1935 paper Elementare Theorie der konvexen Polyeder (Comment. Math. Helv. 7, 290–306) gives an elementary, self-contained proof of both directions. Its first result, which Weyl calls the Hauptsatz (main theorem, Satz 1), is the direction "finitely generated ⇒ finite intersection of half-spaces", in a sharp form: the half-spaces needed are exactly the extreme supports of the generating set, i.e. its facets. Its sharpening, Satz 2, bounds the number of generators needed to represent a point by the dimension $n$. This mission formalizes §§1–2 of the paper (pp. 290–295): the Hauptsatz, its sharpening, and the steps of Weyl's inductive proof.

Timeline:

• 1896, H. Minkowski, Geometrie der Zahlen: polytopes as bounded intersections of half-spaces and as convex hulls of finitely many points.
• 1911, C. Carathéodory: a point in the convex hull of a set in $\mathbb{R}^d$ is a convex combination of at most $d+1$ of its points (Rend. Circ. Mat. Palermo 32).
• 1935, H. Weyl: the present paper; Satz 1 and Satz 2 for cones, with the dual statements in §3 and the polytope theorem in §4.

Setting

Points of $\mathbb{R}^n$ are $n$-tuples $x = (x_1, \ldots, x_n)$, and $\langle \alpha, x \rangle = \alpha_1 x_1 + \cdots + \alpha_n x_n$. A vector $\alpha \ne 0$ determines the half-space $\{x : \langle\alpha,x\rangle \ge 0\}$; positive multiples of $\alpha$ give the same half-space.

A point system $S$ is a finite set of points of $\mathbb{R}^n$. It is non-degenerate if its points do not all satisfy one equation $\langle\alpha,x\rangle = 0$ with $\alpha \neq 0$, i.e. the only $\alpha$ orthogonal to every point of $S$ is $0$.

A half-space $\langle\alpha,x\rangle\ge 0$ ($\alpha\ne 0$) is a support of $S$ if every point of $S$ lies in it. It is an extreme support if, in addition, equality $\langle\alpha,x\rangle = 0$ holds at $n-1$ linearly independent points $x$ of $S$.

A point $x$ is representable by $S$ if it is a nonnegative combination of the points of $S$:

$$x = \sum_{s\in S} c_s\, s, \qquad c_s \ge 0 .$$

The set of points lying in all extreme supports of $S$ is Weyl's konvexe Pyramide. In the Lean development these objects are Representable, NonDegenerate, IsSupport and IsExtremeSupport in the namespace WeylPolyhedra.Pyramid, with points of type Fin n → ℝ and $\langle\alpha,x\rangle$ written α ⬝ᵥ x.

Formalization targets

Goal: Satz 2 (Verschärfung des Hauptsatzes), p. 295

For a finite non-degenerate $S \subset \mathbb{R}^n$ and a point $x$ with $\langle\alpha,x\rangle\ge 0$ for every extreme support $\alpha$ of $S$,

$$\exists\, T \subseteq S,\quad |T| \le n,\quad x = \sum_{t\in T} c_t\, t,\ \ c_t \ge 0 .$$

Satz 1 (Hauptsatz), p. 291

Under the same hypotheses, $x$ is representable by $S$. Satz 2 contains Satz 1.

Steps of the proof (§1–§2)

1. A finite non-degenerate $S$ has only finitely many extreme supports, up to positive scaling (p. 291).
2. The reduction step of case a) (p. 292): if $S$ has an extreme support $\beta$ and $p$ satisfies all extreme supports, there are $e \in S$ with $\langle\beta,e\rangle>0$ and $\lambda\ge 0$ such that $q = p-\lambda e$ still satisfies all extreme supports and lies on the plane of one of them.
3. The lifting step (p. 293): with $x_n \ge 0$ an extreme support of $S$ and $S_0$ the points on $x_n = 0$, every extreme support $\beta$ of $S_0$ in $\mathbb{R}^{n-1}$ lifts to the extreme support $\beta_1x_1+\cdots+\beta_{n-1}x_{n-1} - \mu x_n \ge 0$ of $S$ (inequality (6)).
4. Case b) (p. 291, proved pp. 293–294): if $S$ has no extreme support, every point of $\mathbb{R}^n$ is representable by $S$.

Significance

Satz 1 together with its trivial converse identifies the cone generated by $S$ with the intersection of its extreme-support half-spaces. This is one half of the Minkowski–Weyl theorem for cones, and it names the half-spaces: they are the facets of the cone. Satz 2 adds the conic form of Carathéodory's theorem: every point of a cone generated by a finite spanning set in $\mathbb{R}^n$ is a nonnegative combination of at most $n$ generators. In linear programming this is the statement that a feasible system has a basic feasible solution. The second mission in this series, on §§3–4 of the paper, uses Satz 1 to prove that a bounded region cut out by finitely many inequalities is the convex hull of finitely many points, and conversely.

On formalization status: Mathlib defines finitely generated and dually finitely generated pointed cones (PointedCone, PointedCone.DualFG) and proves Carathéodory's theorem for convex hulls (convexHull_eq_union), but, at the pinned revision, it does not prove the Minkowski–Weyl theorem or the facet description of a finitely generated cone. The results are classical and proved in the paper; this mission produces machine-checked proofs of them, in Weyl's formulation with extreme supports, together with the intermediate steps of his induction.

Difficulty

The hypothesis only controls $x$ against the extreme supports, not against every support. Showing that $x$ lies in the cone generated by $S$ whenever $\langle\alpha,x\rangle\ge 0$ holds for every support is the conic Farkas lemma, which follows from a separating hyperplane argument. Here that argument is not enough: a separating hyperplane is a support, but in general not an extreme one, and the statement is about the finitely many extreme ones. The proof has to produce, for a point outside the cone, a violated extreme support, which requires control over the facet structure of the cone.

The dimension count of Satz 2 is a second difficulty. An induction on the dimension naturally gives $n$ generators in one case and $n+1$ in another (a point of a half-space needs one generator on each side), and Weyl notes that he could not avoid a detour to recover the bound $n$. The case where $S$ has no extreme support at all must also be handled separately; it is not vacuous, since $S$ can then generate all of $\mathbb{R}^n$.

Formalization scope

Conventions committed to in Lean:

• $\mathbb{R}^n$ is Fin n → ℝ; points and normals share this type (the dual space is identified with $\mathbb{R}^n$, as in the paper). The pairing is dotProduct, written α ⬝ᵥ x.
• A point system is a Finset (Fin n → ℝ). The zero vector is not excluded.
• A support normal satisfies α ≠ 0. Extreme supports require a subset T ⊆ S with T.card = n - 1 whose elements are linearly independent in the vector space $\mathbb{R}^n$.
• "All extreme support equations are satisfied" in Satz 1 is read as the inequalities $\langle\alpha,x\rangle\ge0$ for every extreme normal $\alpha$, as the proof and Satz 2 make explicit. The hypothesis quantifies over all extreme normals, so no representatives are chosen.
• "Positive-linear" combinations have nonnegative coefficients (display (3)). In Satz 2 the subset $T$ is not required to be linearly independent.
• Finiteness of extreme supports is stated up to positive scaling.
• The lifting step is stated in the coordinates Weyl fixes on p. 293: $\mathbb{R}^n$ is Fin (m+1) → ℝ, the extreme support is $x_n \ge 0$ (Fin.last m), $S_0$ is projected by Fin.init, and $\mu$ is given together with hypotheses that it is the attained minimum. The hypothesis $n \ge 2$ is made explicit.

Replacing extreme supports by all supports in the hypothesis of Satz 1 or Satz 2 would turn the goal into a much weaker theorem (the conic Farkas lemma plus Carathéodory) and is not an admissible formalization. Dropping non-degeneracy makes Satz 1 false: for $S = \{e_1\} \subset \mathbb{R}^2$ the extreme supports are $\pm x_2 \ge 0$, and $x = (-1, 0)$ satisfies both without being a nonnegative multiple of $e_1$.

A complete development needs basic linear algebra over Fin n → ℝ (hyperplanes through $n-1$ independent points, projection to a coordinate hyperplane) and finite minimisation. The facet description of finitely generated cones, conic Carathéodory and the finiteness of facets are reusable beyond this mission, including for the second mission of the series. Contributions of lemmas on PointedCone that connect Representable with PointedCone.span are welcome.

Selected references

• H. Weyl, Elementare Theorie der konvexen Polyeder, Commentarii Mathematici Helvetici 7 (1935), 290–306. https://doi.org/10.1007/BF01292722
• C. Carathéodory, Über den Variabilitätsbereich der Fourier'schen Konstanten von positiven harmonischen Funktionen, Rendiconti del Circolo Matematico di Palermo 32 (1911), 193–217. https://doi.org/10.1007/BF03014795
• A. Schrijver, Theory of Linear and Integer Programming, Wiley, 1986, §7.2 (the Farkas–Minkowski–Weyl theorem). ISBN 978-0-471-98232-6
• G. M. Ziegler, Lectures on Polytopes, Springer GTM 152, 1995, Lecture 1. https://doi.org/10.1007/978-1-4613-8431-1

Motivation

The traveling salesman problem asks for a shortest closed route visiting every node of a weighted complete graph exactly once. It is NP-hard, so practitioners use fast heuristics, and the basic question about a heuristic is how far from optimal its tour can be. Rosenkrantz, Stearns and Lewis (SIAM J. Comput. 6(3), 1977) gave the first systematic worst-case analysis of the simple constructive heuristics under the triangle inequality: nearest neighbor, the family of insertion methods, and several variants.

Insertion methods build a tour by growing it one node at a time. They are among the most widely used construction heuristics in practice and in textbooks, and they differ only in the rule that chooses which node to insert next: the nearest one, the cheapest one, the farthest one, a random one, or any other. This mission formalizes the paper's result that holds for the whole family at once, regardless of that rule: every insertion method produces a tour at most $\lceil \lg n\rceil + 1$ times longer than an optimal one (Theorem 3, p. 571).

Timeline. 1977: Rosenkrantz, Stearns and Lewis prove $\lceil\lg n\rceil+1$ for every insertion method (Theorem 3), $\tfrac12(\lceil\lg n\rceil+1)$ for nearest neighbor (Theorem 1), both from a shared counting lemma (Lemma 1), and the constant $2$ for nearest and cheapest insertion (Theorem 4). 1994: Bafna, Kalyanasundaram and Pruhs (Theoretical Computer Science 125, 1994) give instances on which some insertion methods reach ratio $\Omega(\log n/\log\log n)$, so the logarithmic growth cannot be replaced by a constant for the family as a whole.

Setting

A traveling salesman graph with $n$ nodes consists of a finite node set $N$ with $|N|=n$ and a distance $d:N\times N\to\mathbb R$ with $d(i,j)=d(j,i)$, $d(i,j)\ge 0$ and $d(i,j)+d(j,k)\ge d(i,k)$ for all nodes (the triangle inequality). A tour visits every node once and returns to its start; its length is the sum of its edge lengths, and OPTIMAL is the least length of a tour.

A subtour is a tour on a subset of the nodes; a single node is a tour without edges. Given a subtour $T$ and a node $k\notin T$, TOUR$(T,k)$ is obtained by choosing an edge $(x,y)$ of $T$ minimizing

$$d(x,k)+d(k,y)-d(x,y)$$

and replacing it by the edges $(x,k)$ and $(k,y)$; if $T$ is a single node $i$, TOUR$(T,k)$ is the two-node tour $(i,k),(k,i)$. COST$(T,k)$ is the length of TOUR$(T,k)$ minus the length of $T$.

An insertion method constructs subtours $T_1,\dots,T_n$ with $T_1=\{a_0\}$ a single node and $T_{i+1}=\mathrm{TOUR}(T_i,a_i)$ for some node $a_i\notin T_i$, $1\le i<n$. The final tour $T_n$ is the approximation, and INSERT denotes its length. No rule for choosing the $a_i$ is fixed, and ties between minimizing edges are broken arbitrarily.

Write $\lg$ for the logarithm to base 2 and $\lceil x\rceil$ for the least integer $\ge x$.

Formalization targets

Goal: Theorem 3

For every traveling salesman graph with $n\ge 1$ nodes and every run of every insertion method,

$$\mathrm{INSERT}\ \le\ \bigl(\lceil\lg n\rceil+1\bigr)\cdot\mathrm{OPTIMAL}.$$

Milestones

1. (2.2), shortcutting: visiting a subset of the nodes in the order of a tour gives a tour of the subset that is no longer.
2. (2.1): if the numbers $l_1\ge\dots\ge l_n$ satisfy $d(p,q)\ge\min(l_p,l_q)$ for distinct $p,q$, then $\mathrm{OPTIMAL}\ge 2\sum_{i=k+1}^{\min(2k,n)} l_i$ for $1\le k\le n$.
3. Lemma 1: if $d(p,q)\ge\min(l_p,l_q)$ for distinct nodes and $l_p\le\frac12\mathrm{OPTIMAL}$ for all $p$, then

$$\sum_p l_p\le\tfrac12\bigl(\lceil\lg n\rceil+1\bigr)\mathrm{OPTIMAL}.$$

4. Lemma 2: $\mathrm{COST}(T,k)\le 2\,d(k,j)$ for every node $j$ of $T$.
5. (3.7): $\mathrm{INSERT}=\sum_{i=1}^{n-1}\mathrm{COST}(T_i,a_i)$.
6. (3.10): $\mathrm{COST}(T_i,a_i)\le 2\,d(a_i,a_j)$ whenever $j<i$.
7. (3.12): $\mathrm{COST}(T_i,a_i)\le\mathrm{OPTIMAL}$ for $1\le i<n$.

Significance

The result. Theorem 3 is a guarantee for an entire class of algorithms rather than for one. Any rule for choosing the next node, including rules designed for speed or for empirical quality, inherits a worst-case ratio of $\lceil\lg n\rceil+1$ from the insertion step alone. The rule matters only for improving on that: nearest and cheapest insertion achieve the constant $2(1-1/n)$ (Theorem 4 and its corollary, the subject of the third mission of this series), while the logarithmic bound remains the best general statement for other rules, such as farthest or arbitrary insertion. Lemma 1 is reusable on its own: it converts "every node carries a charge bounded by half the optimum and by its distance to other nodes" into a logarithmic bound, and the same lemma yields the nearest neighbor bound of Theorem 1.

Formalizing it. The theorem has been proved since 1977; the work here is a machine-checked proof of the known argument together with a reusable library for subtours, insertion and insertion costs. The companion nearest neighbor bound (Theorem 1) is already on the platform as SupplyChainTheory.nearest_neighbor_bound (proved), and nearest insertion with constant 2 as SupplyChainTheory.nearest_insertion_bound; neither covers arbitrary insertion methods or states Lemma 1 separately.

Difficulty

The per-step facts are local: each insertion is cheap relative to a node already present (Lemma 2) and relative to OPTIMAL (3.12). The obvious way to combine them, adding up $n-1$ costs each at most OPTIMAL, gives only the ratio $n-1$. The logarithm comes from a global counting argument over all nodes simultaneously (Lemma 1), in which OPTIMAL is compared with tours on nested subsets of nodes of doubling size, and the per-node charges must be matched against the edges of those tours. Formally, the delicate parts are the bookkeeping of subtours as they grow (that every earlier node lies on the current subtour, and that the insertion cost equals the length increase), the shortcutting of a tour to an arbitrary subset, and the ceiling-of-logarithm arithmetic.

Formalization scope

Nodes are Fin n; a tour of all nodes is a permutation τ : Equiv.Perm (Fin n), and OPTIMAL is the minimum of the tour length over the finite, nonempty set of permutations. Subtours are duplicate-free lists of nodes, with closed length $d(x_0,x_1)+\dots+d(x_{m-1},x_0)$. TOUR$(T,k)$ is encoded as inserting $k$ at a list position whose resulting length is minimal among all positions; inserting at a position removes exactly one edge of $T$ and raises the length by exactly $d(x,k)+d(k,y)-d(x,y)$, so this is the paper's rule, with every tie-breaking allowed. COST is the minimum length increase over positions. The paper's 1-based subtour index is kept ($T_1=[a_0]$, $T_n$ final). $\lceil\lg n\rceil$ is Nat.clog 2 n. All quantities are real.

Conventions and deviations, each disclosed in the item statements:

• The distance satisfies $d(i,i)=0$, a normalization not in the paper; a loop never enters any length.
• Ratios are multiplied out ($\mathrm{INSERT}\le c\cdot\mathrm{OPTIMAL}$), so the paper's exclusion of the identically zero distance (1.1) is not needed.
• Condition a) of Lemma 1 is required for distinct nodes only. The page says "for all nodes $p$ and $q$", which for $p=q$ would force every $l_p\le 0$ and make the lemma inapplicable in the proof of Theorem 3; the proof uses the condition only on edges of a tour.
• (2.2) is stated for every subset of the nodes and every tour, which is what the shortcut argument shows; the paper applies it to one specific subset and an optimal tour.
• (2.1) uses 0-based node labels, so its range $k+1,\dots,\min(2k,n)$ becomes $k,\dots,\min(2k,n)-1$.

The goal quantifies over every run: any choice of the inserted nodes $a_i$ and any minimizing insertion position. Adding a selection rule (nearest, cheapest) or fixing a tie-breaking would state a weaker, different theorem; restricting to instances with OPTIMAL $=0$ or to a fixed small $n$ would trivialize it.

Reusable beyond this mission: the subtour and insertion library (closed length of a list, TOUR, COST, insertion runs) and Lemma 1, which also yields Theorem 1. Contributions welcome: proofs of the milestones, general lemmas about the closed length of List.insertIdx and of filtered lists, and a proof of Theorem 1 from this mission's Lemma 1.

Selected references

• D. J. Rosenkrantz, R. E. Stearns, P. M. Lewis II, An Analysis of Several Heuristics for the Traveling Salesman Problem, SIAM Journal on Computing 6(3):563–581, 1977. https://doi.org/10.1137/0206041
• V. Bafna, B. Kalyanasundaram, K. Pruhs, Not all insertion methods yield constant approximate tours in the Euclidean plane, Theoretical Computer Science 125(2):345–353, 1994.

Motivation

In robust stochastic optimization a decision maker chooses a decision $x$ whose outcome depends on a random vector $\mathbf R$, but knows only the first two moments of $\mathbf R$: its mean vector $\mu$ and its covariance matrix $\Sigma$. The decision is evaluated by its worst-case expected utility over every distribution consistent with those moments. This model is standard in portfolio selection, where estimated means and covariances are the usual inputs, and in pricing and inventory problems with mean-variance information. It goes back to Scarf's min-max newsvendor (1958) and the Chebyshev-type moment bounds of Bertsimas and Popescu (2005).

For a linear outcome $x'\mathbf R$, such as the return of a portfolio with weights $x$, the robust objective is

$$U(x) = \min_{\mathbf R \sim (\mu,\Sigma)} E[u(x'\mathbf R)],$$

an optimization over an infinite-dimensional set of $n$-variate distributions. Popescu (2007) showed that this problem depends on $\mu$ and $\Sigma$ only through the scalar mean $\mu_x = x'\mu$ and variance $\sigma_x^2 = x'\Sigma x$. The multivariate robust problem then reduces to a univariate moment problem, and for many utilities to a parametric quadratic program. The reduction rests on one structural fact, the general projection property, which this mission formalizes.

Setting

Fix a dimension $n$. A law on $\mathbb R^n$ is a Borel probability measure on $\mathbb R^n$. For a vector $\mu \in \mathbb R^n$ and a real $n\times n$ matrix $\Sigma$, the mean-covariance class $\mathbb M^n_{(\mu,\Sigma)}$ is the set of laws $P$ under which every coordinate $R_i$ has a finite second moment and

$$\int R_i\,dP(R) = \mu_i, \qquad \int (R_i-\mu_i)(R_j-\mu_j)\,dP(R) = \Sigma_{ij} \qquad (1\le i,j\le n).$$

Writing $\mathbf R \sim (\mu,\Sigma)$ means that the law of $\mathbf R$ lies in $\mathbb M^n_{(\mu,\Sigma)}$. For $n=1$ the superscript is dropped: for real $m$ and $v$, $\mathbb M_{(m,v)}$ is the set of laws on $\mathbb R$ with finite second moment, mean $m$ and variance $v$.

For a vector $x \in \mathbb R^n$, the $x$-projection sends the law $P$ of $\mathbf R$ to the law of the scalar $\mathbf r = x'\mathbf R$, that is, to the pushforward of $P$ under $R \mapsto x'R$. Write $\mu_x = x'\mu$ and $\sigma_x^2 = x'\Sigma x$. The matrix $\Sigma$ is positive semidefinite, $\Sigma \succeq 0$, when $x'\Sigma x \ge 0$ for all $x$ (and $\Sigma$ is symmetric); $\Sigma^{1/2}$ denotes its positive semidefinite square root.

Formalization targets

Goal: Theorem 1 (General Projection Property)

For every $\mu \in \mathbb R^n$, every $\Sigma \succeq 0$ and every nonzero $x \in \mathbb R^n$, the $x$-projection maps $\mathbb M^n_{(\mu,\Sigma)}$ into and onto $\mathbb M_{(\mu_x,\sigma_x^2)}$:

$$\bigl\{\, \text{law of } x'\mathbf R \;:\; \mathbf R \sim (\mu,\Sigma) \bigr\} \;=\; \mathbb M_{(x'\mu,\; x'\Sigma x)}.$$

The "into" half says every projected law has the right mean and variance. The "onto" half says that every univariate law with mean $\mu_x$ and variance $\sigma_x^2$, however heavy-tailed or irregular, is the law of $x'\mathbf R$ for some $\mathbf R \sim (\mu,\Sigma)$. The degenerate case $x'\Sigma x = 0$ is included.

Milestones

1. The into half (§2.1, justification of (4)): $x'\mathbf R$ has mean $x'\mu$ and variance $x'\Sigma x$.
2. The degenerate case: if $x'\Sigma x = 0$ then $x'\mathbf R = x'\mu$ almost surely.
3. Standardization: if $\mathbf r \sim (m, v)$ with $v > 0$, then $v^{-1/2}(\mathbf r - m) \sim (0,1)$.
4. Normalization: for $x'\Sigma x > 0$, the vector $y = (x'\Sigma x)^{-1/2}\Sigma^{1/2}x$ satisfies $y'y = 1$.
5. Isotropic lift: if $y'y = 1$ and $\mathbf z \sim (0,1)$, there is $\mathbf Z \sim (0, I_n)$ with $y'\mathbf Z$ distributed as $\mathbf z$.
6. Affine image: if $\mathbf Z \sim (0,I_n)$ then $\mu + \Sigma^{1/2}\mathbf Z \sim (\mu,\Sigma)$, and $x'(\mu + \Sigma^{1/2}Z) = x'\mu + (x'\Sigma x)^{1/2}\,y'Z$ for every $Z$.

Significance

The result. Theorem 1 immediately yields Proposition 1 of the paper: for every objective $u$,

$$\min_{\mathbf R\sim(\mu,\Sigma)} E[u(x'\mathbf R)] = \min_{\mathbf r\sim(\mu_x,\sigma_x^2)} E[u(\mathbf r)],$$

with minima in the wide sense of infima. The robust objective is therefore a function of $(\mu_x, \sigma_x)$ alone, which makes every robust mean-covariance problem with a linear outcome a bicriteria mean-variance problem. The paper's later results use this: the two-point and one-point support properties, the parametric quadratic programming solution, and the portfolio applications (bonus schemes, value at risk). The projection property holds with no assumption on $u$, so it serves non-concave, discontinuous and quantile-based objectives alike.

Formalizing it. The theorem is proved in the paper; no machine-checked version is known. The mission produces a formal definition of mean-covariance classes that treats integrability honestly, a proof of the projection property, and through it a formally verified reduction of multivariate moment-robust problems to univariate ones. The paper's own construction of the lifted vector has a gap (see Difficulty), so a formal proof also records a corrected argument.

Difficulty

The into half is a computation with linearity of expectation. The difficulty is entirely in the onto half. Given an arbitrary univariate law with prescribed mean and variance, one must build an $n$-variate law with a prescribed full covariance matrix whose one-dimensional marginal in direction $x$ is exactly the given law. This is a coupling problem: the obvious approach, taking independent coordinates, fixes the marginal in direction $x$ as a convolution and cannot reproduce an arbitrary target. Taking $\mathbf R$ supported on the line through $\mu$ in a single direction reproduces the target law but has a rank-one covariance and fails whenever $\Sigma$ has rank above one.

The paper's appendix constructs the lift through conditional distributions of the remaining coordinates given the projected one. As printed, the conditional second-moment requirement it imposes cannot hold for unbounded targets, so that argument does not go through verbatim. The milestone for the lift states only the claim, not the printed construction.

The integrability bookkeeping is real work: every intermediate law must be shown to have finite second moments before its moments can be computed.

Formalization scope

• $\mathbb R^n$ is EuclideanSpace ℝ (Fin n) with its Borel σ-algebra; $x'R$ is the inner product $\langle x, R\rangle$; $x'\Sigma x$ is x.ofLp ⬝ᵥ S *ᵥ x.ofLp, where the matrix $\Sigma$ is named S (the symbol Σ is reserved in Lean).
• Laws are probability measures. Both classes require finite second moments (MemLp … 2), so that means and covariances are genuine integrals, not the default value $0$ that Lean assigns to non-integrable functions. The univariate class is parametrized by the variance $v = \sigma^2$, not by $\sigma$.
• The projection is the pushforward P.map (fun R => ⟪x, R⟫) under a continuous map. "Pathwise" identities in the paper become equalities of pushforward laws, or pointwise algebraic identities.
• $\Sigma^{1/2}$ is CFC.sqrt S, acting through Matrix.toEuclideanCLM, as in Mathlib's multivariateGaussian.
• The goal is stated as Set.MapsTo ∧ Set.SurjOn with both classes explicit. Its only hypotheses are $\Sigma \succeq 0$ and $x \ne 0$, as in the paper. No bound on $n$, no invertibility of $\Sigma$ and no positivity of $x'\Sigma x$ is assumed. Restricting the target to Gaussian, bounded or finitely supported laws, or dropping the finite-second-moment clause (which would admit Cauchy laws as "mean 0, variance 0"), would trivialize or change the theorem and is ruled out.
• Milestones 3, 4 and 6 assume $x'\Sigma x > 0$ (or $v > 0$), the case the proof treats after its first sentence; milestone 2 covers the complementary case.

Needed infrastructure: moments of pushforwards under linear and affine maps, a covariance calculus for coordinates of random vectors, and a coupling that realizes the isotropic lift. Mathlib's multivariateGaussian, stdGaussian and CFC.sqrt are available. A reusable lemma "the covariance of $A\mathbf Z + b$ is $A\,\mathrm{Cov}(\mathbf Z)A'$" would serve beyond this mission. Related platform work on moment-based ambiguity sets: Wasserstein Distributionally Robust Optimization II. Contributions of any milestone, and alternative proofs of the lift, are welcome.

Selected references

• I. Popescu, Robust Mean-Covariance Solutions for Stochastic Optimization, Operations Research 55(1):98–112, 2007. https://doi.org/10.1287/opre.1060.0353
• D. Bertsimas, I. Popescu, Optimal Inequalities in Probability Theory: A Convex Optimization Approach, SIAM Journal on Optimization 15(3):780–804, 2005. https://doi.org/10.1137/S1052623401399903
• H. Scarf, A Min-Max Solution of an Inventory Problem, in Studies in the Mathematical Theory of Inventory and Production, Stanford University Press, 1958.
• W. W. Rogosinski, Moments of Non-Negative Mass, Proceedings of the Royal Society A 245:1–27, 1958. https://doi.org/10.1098/rspa.1958.0062

Motivation

Interior point methods for a linear program $\min\{c^\top x : Ax\ge b\}$ with $A\in\mathbb R^{m\times n}$ follow the central path of the logarithmic barrier $-\sum_i\log s_i$, where $s=Ax-b$ is the slack vector. Renegar's path-following analysis (1988) gives $O(\sqrt m\,L)$ iterations, and for decades this was the best bound for methods whose iterations cost a linear system solve. Vaidya's volumetric barrier $-\log\det(A^\top S^{-2}A)$ and the hybrid volumetric barriers of Vaidya and of Anstreicher (references [45] and [2] of the paper) reached $O((m\,\mathrm{rank}(A))^{1/4}L)$ iterations at the price of more expensive linear algebra. Nesterov and Nemirovski showed that a universal barrier gives $O(\sqrt n\,L)$ iterations, but that barrier cannot be evaluated efficiently.

Lee and Sidford (FOCS 2014; full version arXiv:1312.6677) obtained $\tilde O(\sqrt{\mathrm{rank}(A)}\,L)$ iterations, each costing $\tilde O(1)$ linear system solves, by following a weighted central path whose weights are recomputed from the slacks. The weights come from a weight function $g$, defined as the minimizer of a regularized D-optimal-design problem. This mission is about that weight function and the theorem (Theorem 1 of the paper) certifying its properties. The companion mission, Path-Finding Methods for Linear Programming I, formalizes the path-following framework (Theorem 5 of §IV.C) that consumes these properties.

Setting

Fix $A\in\mathbb R^{m\times n}$ with full column rank, $\mathrm{rank}(A)=n$, and $1\le n<m$. For vectors $s,w\in\mathbb R^m_{>0}$ write $S=\mathrm{diag}(s)$, $W=\mathrm{diag}(w)$, $W^\alpha=\mathrm{diag}(w_i^\alpha)$, and $A_s=S^{-1}A$. For a matrix $M$ let $\|v\|_M=\sqrt{v^\top Mv}$.

Projection matrix and slack sensitivity (Definition 2, p. 428). The projection matrix is $P_{S^{-1}A}(w)=W^{1/2}S^{-1}A\,(A^\top S^{-1}WS^{-1}A)^{-1}A^\top S^{-1}W^{1/2}$, and the slack sensitivity is

$$\gamma(s,w)=\max_{i\in[m]}\big\|W^{-1/2}\mathbb 1_i\big\|_{P_{S^{-1}A}(w)} .$$

Weight function (Definition 4, p. 428). A map $g:\mathbb R^m_{>0}\to\mathbb R^m_{>0}$ is a weight function with constants $c_1,c_\gamma,c_r$ if it is differentiable and, for every $s>0$, with $G(s)=\mathrm{diag}(g(s))$, $G'(s)$ the Jacobian of $g$ at $s$, and $\|y\|_{G(s)}=\sqrt{\sum_ig_i(s)y_i^2}$:

1. Size: $\|g(s)\|_1\le c_1$;
2. Slack sensitivity: $c_\gamma\ge1$ and $\gamma(s,g(s))\le c_\gamma$;
3. Step consistency: $c_r\ge1$ and for all $r\ge c_r$, $y\in\mathbb R^m$: $\|(I+r^{-1}G^{-1}G'S)y\|_{G(s)}\le\|y\|_{G(s)}$ and $\|y+r^{-1}G^{-1}G'Sy\|_\infty\le\|y\|_\infty+c_r\|y\|_{G(s)}$;
4. Uniformity: $\|g(s)\|_\infty\le2$.

The regularized objective (6), p. 429. For $\alpha,\beta\in\mathbb R$,

$$\hat f(s,w)=\mathbb 1^\top w-\frac1\alpha\log\det\big(A_s^\top W^\alpha A_s\big)-\beta\sum_{i\in[m]}\log w_i ,\qquad g(s)=\arg\min_{w\in\mathbb R^m_{>0}}\hat f(s,w).$$

At $\alpha=1,\beta=0$ this is the D-optimal design problem, dual to computing the John ellipsoid of the polytope $\{y:|[A(y-x)]_i|\le s_i\}$ (§V.B).

Formalization targets

Goal: Theorem 1 (Properties of Weight Function), §V.A, p. 429

With

$$\alpha=1-\Big(\log_2\frac{2m}{\mathrm{rank}(A)}\Big)^{-1},\qquad \beta=\frac{\mathrm{rank}(A)}{2m},$$

the objective $\hat f(s,\cdot)$ has a unique minimizer over $\mathbb R^m_{>0}$ for every $s>0$, and the resulting $g$ is a weight function with

$$c_1(g)=2\,\mathrm{rank}(A),\qquad c_\gamma(g)=2,\qquad c_r(g)=2\log_2\frac{2m}{\mathrm{rank}(A)} .$$

Milestones: the three bullets of Theorem 1

• Size: every minimizer $w$ of $\hat f(s,\cdot)$ satisfies $\|w\|_1\le2\,\mathrm{rank}(A)$.
• Slack sensitivity: every minimizer $w$ satisfies $\gamma(s,w)\le2$.
• Step consistency: any map $g$ selecting a minimizer at every $s>0$ is differentiable on $\mathbb R^m_{>0}$ and satisfies the two step-consistency inequalities for every $r\ge2\log_2\frac{2m}{\mathrm{rank}(A)}$.

A supporting (non-milestone) item states the existence and uniqueness of the minimizer on its own.

Significance

The result. Theorem 1 is the input that turns the weighted path-following framework into an $\tilde O(\sqrt{\mathrm{rank}(A)}\,L)$-iteration method: the framework needs $O(c_\gamma^{-1}c_r^{-3}c_1^{-1/2})$-sized steps in $t$ (p. 428), and Theorem 1 makes that $\tilde\Omega(1/\sqrt{\mathrm{rank}(A)})$. The step consistency bound is what allows the weights to be recomputed after each Newton step without losing centrality. The same construction underlies later work on Lewis-weight barriers and on fast approximate John ellipsoids and maximum flow (§VIII of the paper).

Formalizing it. The theorem is proved in the full version of the paper (arXiv:1312.6677); the FOCS extended abstract contains no proofs. No part of it has a machine-checked proof. A complete formalization would give a verified account of leverage-score calculus (sums of leverage scores equal the rank; derivatives of projection matrices), of the convexity of $w\mapsto-\log\det(A^\top W^\alpha A)$ for $\alpha\in(0,1)$, and of differentiability of an argmin via the implicit function theorem, none of which is currently packaged in Mathlib in this form.

Difficulty

Size and slack sensitivity are statements about the minimizer, which is only characterized implicitly; they require precise matrix calculus for $\log\det(A_s^\top W^\alpha A_s)$ and a comparison between the matrices $A^\top WA$ (which defines $\gamma$) and $A^\top W^\alpha A$ (which defines $g$). The specific values of $\alpha$ and $\beta$ matter here: the unregularized choice $\alpha=1$, $\beta=0$ makes the problem degenerate (p. 429).

The hard part is step consistency. The Jacobian $G'$ of an argmin is available only implicitly, as the solution of a linear system obtained by differentiating the optimality condition. A bound on $\|G'\|$ that depends on $m$ is easy to get and useless: the theorem needs the operator norm of $I+r^{-1}G^{-1}G'S$ in the $G(s)$-norm to be at most $1$ as soon as $r$ exceeds $2\log_2(2m/\mathrm{rank}(A))$, and an $\ell_\infty$ bound with only an additive $c_r\|y\|_{G(s)}$ loss.

Existence and differentiability of the minimizer are conclusions, not hypotheses. The minimization is over an open orthant on which the objective is not obviously coercive or strictly convex for $\alpha<1$, and differentiability of $g$ requires the Hessian of $\hat f$ at the minimizer to be invertible.

Formalization scope

Vectors are Fin m → ℝ, matrices Matrix (Fin m) (Fin n) ℝ; inverses are Matrix.inv, $\log\det$ is Real.log (Matrix.det …), $w_i^\alpha$ is Real.rpow, $\log_2$ is Real.logb 2, the Jacobian is fderiv ℝ g s, and $\|\cdot\|_\infty$ is Mathlib's sup norm on Fin m → ℝ.

Conventions and pinned hypotheses:

• Full column rank A.rank = n is assumed in every theorem. The paper never states it, but without it $A_s^\top W^\alpha A_s$ is singular and every formula is undefined (in Lean, Matrix.inv and Real.log would return junk $0$).
• $1\le n<m$. $\beta=\mathrm{rank}(A)/(2m)$ and $\log_2(2m/\mathrm{rank}(A))$ need $\mathrm{rank}(A)\ge1$; at $m=\mathrm{rank}(A)$ the page's $\alpha$ is $0$ and $1/\alpha$ in (6) is undefined.
• Reading of $\alpha$: the exponent $-1$ is the reciprocal of $\log_2\frac{2m}{\mathrm{rank}(A)}$, giving $\alpha\in(0,1)$.
• Size is an upper bound $\|g(s)\|_1\le c_1$ (the paper's weight function has $\|g(s)\|_1=\tfrac32\mathrm{rank}(A)$, while Theorem 1 reports $c_1=2\,\mathrm{rank}(A)$).
• The first step-consistency bullet (an operator-norm bound) is stated for every vector $y$.
• $g$ is any map $\mathbb R^m\to\mathbb R^m$ whose value at each positive $s$ minimizes $\hat f(s,\cdot)$ over $\mathbb R^m_{>0}$. Only its values on the open orthant matter. The goal also asserts that such minimizers exist and are unique, so it is not vacuous.

Ruling out trivializations: the goal does not assume $g$ to be a weight function or to be differentiable, and it does not replace $g$ by an arbitrary weight function; differentiability is a conclusion (a predicate using fderiv without it would make step consistency hold vacuously wherever $g$ fails to be differentiable).

Useful infrastructure, reusable beyond this mission: leverage scores and their sum; derivatives of $w\mapsto\log\det(A^\top WA)$ and of projection matrices; convexity of $-\log\det(A^\top W^\alpha A)$ in $w$ (related to the published ConvexOptimization.log_det_concaveOn); differentiability of the argmin of a strictly convex smooth function. Contributions of these as separate theorems are welcome, as is a proof of any single bullet of Theorem 1.

Selected references

• Y. T. Lee, A. Sidford, Path Finding Methods for Linear Programming: Solving Linear Programs in Õ(√rank) Iterations and Faster Algorithms for Maximum Flow, FOCS 2014, pp. 424–433. https://doi.org/10.1109/FOCS.2014.52
• Y. T. Lee, A. Sidford, Path Finding I: Solving Linear Programs with Õ(√rank) Linear System Solves, arXiv, 2013. https://arxiv.org/abs/1312.6677
• J. Renegar, A polynomial-time algorithm, based on Newton's method, for linear programming, Mathematical Programming 40 (1988). https://doi.org/10.1007/BF01580724