Operations Research

Motivation

Two-stage adaptive optimization models decisions taken in two rounds: a first-stage decision $x$ is fixed before an uncertain parameter is revealed, and a second-stage decision $y(b)$ is chosen after the parameter $b$ is observed, so that the second stage may depend on $b$ arbitrarily. The objective is the worst case over an uncertainty set $\mathcal U$ of possible parameters. Such models arise in robust network design, capacity planning and two-stage covering problems, and they generalize the two-stage robust combinatorial problems (set cover, facility location) studied by Dhamdhere, Goyal, Ravi and Singh.

Optimizing over all functions $y(\cdot)$ is intractable in general: Bertsimas and Goyal note that the optimal second stage is piecewise linear in $b$ with possibly exponentially many pieces (Bemporad, Borrelli and Morari, 2003). A standard remedy, introduced for robust linear programs by Ben-Tal, Goryashko, Guslitzer and Nemirovski (2004), restricts the second stage to affine policies (linear decision rules) $y(b)=Pb+q$; the best affine policy is computable by a single convex program and performs well empirically. The question is when this restriction loses nothing.

Timeline of the relevant results:

• 2004. Ben-Tal, Goryashko, Guslitzer and Nemirovski introduce affinely adjustable robust counterparts and show that the best affine policy is tractable for many uncertainty sets (doi:10.1007/s10107-003-0454-y).
• 2010. Bertsimas, Iancu and Parrilo prove that affine policies are optimal for a class of multistage robust problems with one-dimensional uncertainty per stage and box uncertainty sets (doi:10.1287/moor.1100.0444).
• 2012. Bertsimas and Goyal, the source of this mission, prove that an affine policy is optimal for model (1) whenever $\mathcal U$ is a simplex (Theorem 1), and show that this exactness breaks down for slightly larger sets (doi:10.1007/s10107-011-0444-4).

Setting

Let $A\in\mathbb R^{m\times n_1}$, $B\in\mathbb R^{m\times n_2}$, $c\in\mathbb R^{n_1}_+$ and $d\in\mathbb R^{n_2}_+$. The problem $\Pi_{Adapt}(\mathcal U)$ of model (1) is

$$z_{Adapt}(\mathcal U)=\min\ c^Tx+\max_{b\in\mathcal U}d^Ty(b)\quad\text{s.t.}\quad Ax+By(b)\ge b,\ \ x\ge 0,\ \ y(b)\ge 0\quad\forall b\in\mathcal U,$$

where inequalities between vectors are componentwise. A pair $(x,y)$ satisfying the constraints is feasible; its worst-case cost is $c^Tx+\max_{b\in\mathcal U}d^Ty(b)$. A feasible pair is optimal when its worst-case cost equals $z_{Adapt}(\mathcal U)$, and an affine policy is a second stage of the form $y(b)=Pb+q$ with $P\in\mathbb R^{n_2\times m}$, $q\in\mathbb R^{n_2}$, still required to be nonnegative on $\mathcal U$. The value $z_{Aff}(\mathcal U)$ is the same minimum restricted to affine policies.

A simplex in $\mathbb R^m$ is the convex hull

$$\mathcal U=\operatorname{conv}(b^1,\dots,b^{m+1})$$

of $m+1$ affinely independent points, that is, points for which $b^1-b^{m+1},\dots,b^m-b^{m+1}$ are linearly independent. The proof works with the $m\times m$ matrix $Q=[(b^1-b^{m+1})\cdots(b^m-b^{m+1})]$, the matrix $Y=[(y^*(b^1)-y^*(b^{m+1}))\cdots(y^*(b^m)-y^*(b^{m+1}))]$ of display (2), and the affine rule $\tilde y(b)=YQ^{-1}(b-b^{m+1})+y^*(b^{m+1})$. In Lean these are Qmat v, Ymat v g and interpolant v g, with vertices v : Fin (m+1) → Fin m → ℝ.

Formalization targets

Goal: Theorem 1

If $\mathcal U=\operatorname{conv}(b^1,\dots,b^{m+1})$ with affinely independent $b^j\in\mathbb R^m_+$ and $\Pi_{Adapt}(\mathcal U)$ is feasible, then there exist $\hat x$, $P\in\mathbb R^{n_2\times m}$ and $q\in\mathbb R^{n_2}$ such that

$$(\hat x,\ \hat y),\qquad \hat y(b)=Pb+q\ \ (b\in\mathcal U),$$

is an optimal solution of $\Pi_{Adapt}(\mathcal U)$, optimal among all (not only affine) two-stage solutions. In particular $z_{Aff}(\mathcal U)=z_{Adapt}(\mathcal U)$.

Milestones

The proof of Theorem 1 has no numbered lemma; the milestones are its displayed steps, in attack order:

1. $Q$ is invertible (PDF p. 6).
2. For $b=\sum_j\alpha_jb^j$ with $\sum_j\alpha_j=1$: $Q^{-1}(b-b^{m+1})=(\alpha_1,\dots,\alpha_m)^T$ (PDF p. 6).
3. $\tilde y\big(\sum_j\alpha_jb^j\big)=\sum_j\alpha_j\,y^*(b^j)$ (PDF pp. 6–7).
4. Displays (3)–(5): for any feasible $(x^*,y^*)$, the pair $(x^*,\tilde y)$ is feasible and every bound on the worst-case cost of $(x^*,y^*)$ also bounds that of $(x^*,\tilde y)$ (PDF p. 7).

Significance

The result. Theorem 1 identifies a class of uncertainty sets on which the tractable affine restriction is exact, for every constraint matrix $A$ and $B$ and every nonnegative cost. It is the positive anchor of the paper: Sections 3 and 4 show that with $m+3$ extreme points the best affine policy can already be worse by a factor $2-\delta$, and that on sets with exponentially many extreme points the gap can be $\Omega(m^{1/2-\delta})$; Section 6 uses a dominating simplex, on which affine policies are exact, to build an $O(\sqrt m)$-approximation for general $\mathcal U$. The theorem also says that on a simplex the whole adaptive problem reduces to $m+1$ scenario copies of a linear program.

Formalizing it. The result is proved on paper; no machine-checked version is known on Prove2Me. This mission produces the model (1) in Lean, the barycentric-coordinate identity for a simplex in matrix form, and a statement of optimality that asserts attainment of the minimum in (1), which the paper's proof takes for granted.

Difficulty

Two steps are not routine to formalize. First, the paper starts from "an optimal solution $x^*,y^*(b)$", that is, it assumes the minimum in (1) is attained. Over arbitrary functions $y(\cdot)$ this is not automatic; on a simplex it follows because the problem reduces to a finite linear program on the vertices, whose optimum is attained, but Mathlib has no theory of linear-programming attainment, so this reduction has to be built. Second, the affine-independence step needs the passage from affine independence of $m+1$ points to invertibility of the $m\times m$ matrix $Q$, and the identity $Q^{-1}(b-b^{m+1})=\alpha$ requires the barycentric coordinates and the inverse matrix to be matched index by index. The naive idea of comparing $z_{Aff}$ and $z_{Adapt}$ as real infima does not prove the goal: equality of the two infima says nothing about the existence of an optimal solution.

Formalization scope

Vectors in $\mathbb R^m$ are Fin m → ℝ, with the componentwise order; matrices are Matrix (Fin m) (Fin n) ℝ. The paper's indices start at $1$, Lean's at $0$: $b^j$ is v (j-1) and $b^{m+1}$ is v (Fin.last m). The simplex is convexHull ℝ (Set.range v); it is compact, convex and, by affine independence, full-dimensional, so these standing assumptions of (1) are not stated separately. Nonnegativity of $\mathcal U$ is the hypothesis that all $m+1$ vertices are nonnegative (the page writes $j=1,\dots,m$, a slip for $m+1$). Feasibility of (1) is a hypothesis, as the paper assumes. Optimality (IsOptimalAdapt) means: feasible, and every worst-case cost bound achieved by any feasible two-stage solution is achieved by this one. The values $z_{Adapt}$ and $z_{Aff}$ are infima of the sets of achievable bounds; they are provided for reference and the goal does not depend on them.

The goal must not be replaced by $z_{Aff}(\mathcal U)\le z_{Adapt}(\mathcal U)$, by optimality among affine policies only, or by a version that assumes an optimal solution exists: each of these drops the content "there is an optimal solution and it is affine". The goal does not mention $Q$, $Y$ or the interpolant.

A complete development needs: linear-programming attainment for a finite system of linear inequalities with a cost bounded below (reusable well beyond this mission), the linear-algebra lemmas relating affine independence to an invertible edge matrix (reusable for barycentric coordinates in general), and the convex-hull representation of points of a simplex. Contributions of any of these as separate lemmas are welcome.

Selected references

• D. Bertsimas, V. Goyal, On the power and limitations of affine policies in two-stage adaptive optimization, Math. Program. Ser. A, 2012. doi:10.1007/s10107-011-0444-4
• A. Ben-Tal, A. Goryashko, E. Guslitzer, A. Nemirovski, Adjustable robust solutions of uncertain linear programs, Math. Program. 99(2), 351–376, 2004. doi:10.1007/s10107-003-0454-y
• D. Bertsimas, D. A. Iancu, P. A. Parrilo, Optimality of affine policies in multistage robust optimization, Math. Oper. Res. 35(2), 363–394, 2010. doi:10.1287/moor.1100.0444
• A. Bemporad, F. Borrelli, M. Morari, Min–max control of constrained uncertain discrete-time linear systems, IEEE Trans. Autom. Control 48(9), 1600–1606, 2003. doi:10.1109/TAC.2003.816984

Motivation

A linear program is solved with data that, in practice, is rarely known exactly: coefficients come from measurements, estimates or forecasts. Robust optimization asks for a solution that remains feasible for every realization of the data in a prescribed uncertainty set, and among those the one with the best guaranteed objective value. Ben-Tal and Nemirovski introduced this framework for linear programming in Robust solutions of uncertain linear programs (Oper. Res. Lett. 25, 1999), following their treatment of robust convex optimization (Math. Oper. Res. 23, 1998) and Soyster's earlier work on inexact linear programming (Oper. Res. 21, 1973). The robust counterpart has since become the starting point of a large literature on uncertainty sets, budgets of uncertainty and adjustable policies.

A natural first objection is that the robust counterpart might be needlessly conservative: by demanding feasibility for all realizations simultaneously, it could be infeasible, or have a worse value, even when every individual realization is perfectly well behaved. This mission formalizes the paper's answer (§2.2): under two structural hypotheses, the robust counterpart is no worse than the worst realization.

Setting

Fix $c, f \in \mathbb R^n$ and write a linear program in the homogeneous form (6)

$$(P)\qquad \min\{c^{T}x \mid Ax \ge 0,\ f^{T}x = 1\},$$

where $A$ is a real $m\times n$ matrix and $Ax\ge0$ is componentwise. Every linear program can be put in this form. The matrix $A$ is uncertain: it is only known to lie in an uncertainty set $\mathcal U$ of $m\times n$ matrices. Each $A\in\mathcal U$ gives an instance $(P)$ with feasible set $\{x\mid Ax\ge0,\ f^{T}x = 1\}$ and optimal value $c^*(P)$; the family of instances is $\mathcal P$. The robust counterpart (7) is

$$(P_{\mathcal U})\qquad \min\{c^{T}x \mid x \in G_{\mathcal U}\},\qquad G_{\mathcal U} = \{x\mid Ax\ge0\ \ \forall A\in\mathcal U;\ f^{T}x = 1\},$$

and its optimal value is $c^*$. Since $G_{\mathcal U}$ does not change when $\mathcal U$ is replaced by its closed convex hull, the paper assumes throughout that $\mathcal U$ is convex and closed.

Let $\mathcal U_i\subseteq\mathbb R^n$ be the set of all realizations of the $i$-th row, the projection of $\mathcal U$ onto the data of the $i$-th constraint. The uncertainty is constraint-wise if $\mathcal U = \mathcal U_1\times\dots\times\mathcal U_m$: the rows vary independently. The Boundedness Assumption asks for a convex compact set $Q\subseteq\mathbb R^n$ that contains the feasible set of every instance.

Formalization targets

Goal: Proposition 2.1 (p. 5)

If the uncertainty is constraint-wise and the Boundedness Assumption holds, then

1. $(P_{\mathcal U})$ is infeasible if and only if some instance is infeasible:

$$G_{\mathcal U} = \emptyset \iff \exists A\in\mathcal U:\ \{x\mid Ax\ge0,\ f^{T}x=1\}=\emptyset;$$

2. if $(P_{\mathcal U})$ is feasible with optimal value $c^*$, then

$$c^* = \sup\{c^*(P)\mid (P)\in\mathcal P\}. \tag{9}$$

Milestones

The milestones follow the paper's proof: the row-wise description (8) of robust feasibility; the inclusion of $G_{\mathcal U}$ in every instance's feasible set; the reduction of the semi-infinite system (8) on $Q$ to a finite subsystem; the statement that the finite system (10) $A_1x\ge0,\dots,A_Nx\ge0,\ f^{T}x=1$ then has no solution at all; the Farkas certificate (11); the construction of one infeasible instance from it; and part (i) alone, which part (ii) uses for an augmented program.

Companions

The §2.2 example (every instance has optimal value 1, the robust counterpart is infeasible), and the two invariance remarks: $G_{\mathcal U}$ is unchanged under passing to the closed convex hull of $\mathcal U$ (§2.1) or to the product $\mathcal U_1\times\dots\times\mathcal U_m$ of its projections (§2.2).

Significance

Proposition 2.1 says that, for constraint-wise uncertainty, robustness costs nothing beyond what the worst realization already costs: the robust counterpart is feasible exactly when every instance is, and its optimal value equals the worst instance value. The §2.2 example shows the hypothesis cannot be dropped: there, correlated uncertainty in two rows makes every instance solvable with value 1 while the robust counterpart is infeasible. Together with the invariance of $G_{\mathcal U}$ under passing to the product of projections, this explains why row-wise (constraint-wise) uncertainty sets are the standard modelling choice in robust linear optimization.

The result is proved in the paper; no machine-checked version is known to exist. Formalizing it produces a reusable development of semi-infinite linear systems: the compactness reduction to finite subsystems, a homogeneous Farkas alternative, and the row-averaging argument that uses convexity and the product structure of $\mathcal U$.

Difficulty

The robust counterpart has a continuum of constraints, one for each $A\in\mathcal U$, so Farkas' Lemma cannot be applied to it directly. The step that requires care is passing from infeasibility of this semi-infinite system to infeasibility of a single instance. Compactness yields only finitely many instances whose joint system has no solution in $Q$; those instances are in general all feasible individually, and the infeasible instance has to be manufactured from their rows. Without constraint-wise uncertainty the manufactured matrix need not lie in $\mathcal U$, which is exactly what the §2.2 example exploits. Part (ii) needs the optimal values of the instances to be attained on compact feasible sets, which is where the Boundedness Assumption enters again.

Formalization scope

Vectors are Fin n → ℝ, matrices Matrix (Fin m) (Fin n) ℝ, and $Ax\ge0$ is 0 ≤ A *ᵥ x in the componentwise order. The $i$-th row of $A$ is A i and $a^{T}x$ is a ⬝ᵥ x. The projections $\mathcal U_i$ are the images of $\mathcal U$ under $A\mapsto A_i$, not free sets, and constraint-wise uncertainty is the inclusion $\mathcal U_1\times\dots\times\mathcal U_m\subseteq\mathcal U$ (the reverse inclusion always holds). The Boundedness Assumption keeps both convexity and compactness of $Q$, as on the page.

Optimal values are infima: $c^*$ is the greatest lower bound (IsGLB) of $c^{T}x$ over $G_{\mathcal U}$, and (9) states that $c^*$ is the least upper bound (IsLUB) of the set of real optimal values of the instances. No real sInf/sSup is used, so no junk value can make the statement true.

The goal carries the paper's standing assumption that $\mathcal U$ is convex and closed, and one disclosed addition: $\mathcal U$ is nonempty. The paper takes this for granted; without it part (i) fails for $f = 0$ and the supremum in (9) ranges over the empty set. The goal does not assume that the robust counterpart or any instance attains its optimum, and it does not mention finite subsystems, multipliers or the averaged matrix; those appear only in the milestones. A formalization in which the uncertainty sets $\mathcal U_i$ are arbitrary sets with $\mathcal U = \prod_i\mathcal U_i$, or in which optimal values are taken as sInf without boundedness, would not be faithful and is ruled out.

A complete development needs: compactness arguments for families of closed half-spaces, a Farkas alternative for homogeneous systems with one normalizing equation, and elementary convexity of linear images. These pieces are general and reusable beyond robust optimization. Proofs of the milestones, alternative arguments (for instance via LP duality for part (ii)) and proofs of the companion statements are welcome.

Selected references

• A. Ben-Tal, A. Nemirovski, Robust solutions of uncertain linear programs, Operations Research Letters 25(1):1–13, 1999. https://doi.org/10.1016/s0167-6377(99)00016-4 (cited here by the pages of the authors' manuscript).
• A. Ben-Tal, A. Nemirovski, Robust convex optimization, Mathematics of Operations Research 23(4):769–805, 1998. https://doi.org/10.1287/moor.23.4.769
• A. L. Soyster, Convex programming with set-inclusive constraints and applications to inexact linear programming, Operations Research 21(5):1154–1157, 1973. https://doi.org/10.1287/opre.21.5.1154

Motivation

The conditional logit model is the workhorse of discrete choice analysis in transportation, marketing, labour and industrial organization. McFadden's 1974 chapter derived it from a theory of population choice behaviour and showed how to estimate it by maximum likelihood; this line of work was recognized by his 2000 Nobel Prize in Economic Sciences, awarded for theory and methods of discrete choice analysis. Every applied logit estimation rests on a basic question: does the maximum likelihood estimate exist for the sample at hand? In small samples it may not. When one alternative is always chosen whenever it is available, the likelihood keeps increasing as a parameter tends to infinity, and numerical optimizers report diverging coefficients. This failure is known in the binary case as complete or quasi-complete separation. McFadden's Lemma 3 gives the exact condition, for the multinomial conditional logit model with general alternative sets, under which a maximizer exists.

Timeline. Berkson (1951, 1955) popularized binomial logit; multinomial versions were developed by Gurland (1960), Bloch (1967), Rassam (1971), McFadden (1968) and Theil (1969, 1970). McFadden (1974) stated the existence criterion for the conditional logit likelihood (Lemma 3) together with a quadratic-programming test for it (Lemma 4). Albert and Anderson (1984) later classified separation patterns for binary and multinomial logistic regression, and Haberman (1974) treated existence for log-linear models.

Setting

A choice experiment has $N \ge 1$ trials. Trial $n$ offers an alternative set of $J_n$ alternatives, indexed $i = 1,\dots,J_n$, each described by an attribute vector $z_{in} \in \mathbb{R}^K$ (the values of $K$ specified functions of the individual's and the alternative's characteristics). Trial $n$ is repeated $R_n \ge 1$ times, and alternative $i$ is chosen $S_{in}$ times, so $R_n = \sum_{j} S_{jn}$.

For a parameter $\theta \in \mathbb{R}^K$, with $z_{in}\theta$ the inner product, the selection probabilities are

$$P_{in}(\theta) = \frac{e^{z_{in}\theta}}{\sum_{j=1}^{J_n} e^{z_{jn}\theta}} \qquad (16)$$

and the log-likelihood of the sample is

$$L(\theta) = C - \sum_{n=1}^N \sum_{i=1}^{J_n} S_{in} \log \sum_{j=1}^{J_n} e^{(z_{jn} - z_{in})\theta}, \qquad C = \sum_{n=1}^N \Big[\log R_n! - \sum_{j=1}^{J_n}\log S_{jn}!\Big]. \qquad (18)$$

Write $\bar z_n(\theta) = \sum_i z_{in}P_{in}(\theta)$ for the probability-weighted mean attribute vector of trial $n$.

Axiom 5 (Full Rank). The $\big(\sum_n J_n\big)\times K$ matrix with rows $z_{in} - \bar z_n$ has rank $K$.

Axiom 6. There is no nonzero $\gamma \in \mathbb{R}^K$ with $S_{in}(z_{jn} - z_{in})\gamma \le 0$ for all $i, j = 1,\dots,J_n$ and $n = 1,\dots,N$. Equivalently, no nonzero direction makes every observed choice weakly best in its alternative set.

Formalization targets

Goal: Lemma 3

Under Axiom 5,

$$\big(\exists\, \hat\theta \in \mathbb{R}^K,\ \forall \theta,\ L(\theta) \le L(\hat\theta)\big) \iff \text{Axiom 6}.$$

Milestones

1. Equation (19): the gradient $\partial L/\partial\theta = \sum_n \sum_j (S_{jn} - R_nP_{jn}) z_{jn}$.
2. Equation (20): the Hessian $\partial^2L/\partial\theta\,\partial\theta' = -\sum_n R_n \sum_j (z_{jn} - \bar z_n)'P_{jn}(z_{jn} - \bar z_n)$.
3. $L$ is concave, and every critical point is a global maximizer.
4. A Hessian that is nonsingular everywhere makes $L$ strictly concave with at most one maximizer.
5. Axiom 5 holds at $\theta$ if and only if the Hessian at $\theta$ is negative definite.
6. Necessity: under Axiom 5, a maximizer forces Axiom 6.
7. Equation (21): under Axiom 6, $b(\gamma) = \max_n \max_{i,j} S_{in}(z_{jn}-z_{in})\gamma$ has a positive lower bound $b^*$ on the unit sphere.
8. The bound $L(\theta) - C \le -b^*|\theta|$ for all $\theta$.
9. Sufficiency: Axiom 6 gives a maximizer.

Significance

Lemma 3 tells the practitioner when the conditional logit maximum likelihood estimate exists, before any numerical optimization is attempted. It is a linear-inequality condition on the data alone, so it can be checked by linear or quadratic programming (Lemma 4 of the same paper). The existence of the estimator is also the first step of McFadden's asymptotic theory: Lemma 5 shows that Axiom 6 holds with probability tending to one, and Lemma 6, consistency and asymptotic normality, concerns the estimator whose existence Lemma 3 characterizes. The concavity and Hessian formulas (19)–(20) are the basis of the Newton–Raphson computation of the estimator and of its asymptotic covariance matrix.

The result has been proved since 1974 and is classical. To our knowledge it has no machine-checked proof; Mathlib has no statement about the existence of logit or softmax-regression maximum likelihood estimates. Formalizing it produces a verified existence criterion for the multinomial logit likelihood, verified gradient and Hessian formulas for log-sum-exp likelihoods with repeated observations, and a verified link between full column rank and strict concavity.

Difficulty

The likelihood is concave, and concave functions on $\mathbb{R}^K$ need not attain their supremum. Concavity alone therefore gives nothing, and existence must come from a growth condition. The obvious approach, "the likelihood is bounded above by $C$, hence attains its maximum", fails: $L$ is bounded but can approach its supremum only at infinity, which is exactly the separation case. Sufficiency needs a quantitative rate at which $L$ decreases, uniform over all directions; a direction-by-direction argument does not suffice. Necessity requires strict concavity, which is where Axiom 5 and the requirement that every trial be observed enter. A trial with $R_n = 0$ can supply the rank of Axiom 5 while contributing nothing to $L$, so with such a trial necessity fails. The calculus part, (19)–(20), involves differentiating sums of log-sum-exp terms over dependent index types and identifying the result with a weighted covariance operator.

Formalization scope

• Representation. $\mathbb{R}^K$ is EuclideanSpace ℝ (Fin K), so $|\theta| = (\theta'\theta)^{1/2}$ is the Euclidean norm and $z\theta$ is the inner product ⟪z, θ⟫. Trials are Fin N, alternatives of trial $n$ are Fin (J n), and the counts $S_{in}$ are natural numbers.
• Data structure. The structure Data K bundles $N$, $J$, $z$, $S$ and the standing assumptions $N \ge 1$ and $R_n = \sum_i S_{in} \ge 1$ for every trial; these make the trial and alternative index sets nonempty.
• Axioms 1–4 are built in. The model is the logit form (16) with $v$ linear in $\theta$ (Axiom 4), so "Suppose Axioms 1–5 hold" becomes "Data plus Axiom 5".
• Axiom 5 is read at every $\theta$. The row space of the matrix does not depend on $\theta$.
• Hessian. The Hessian is the Fréchet derivative of the gradient vector field (19), as a continuous linear map.
• The maximizer is global over all of $\mathbb{R}^K$. Neither a local maximizer nor "$L(\hat\theta) \ge L(0)$" is acceptable as the goal; that would make it trivial.
• Infrastructure. Gradients and Hessians of log-sum-exp with dependent finite index types; positive definiteness from full column rank; attainment of the maximum of a coercive continuous function on a finite-dimensional space. The calculus lemmas are reusable for any multinomial logit or softmax likelihood. Missions 4 and 5 of this series reuse the same model. Contributions of general log-sum-exp lemmas, independent of this mission's definitions, are welcome.

Selected references

• D. McFadden, Conditional logit analysis of qualitative choice behavior, in P. Zarembka (ed.), Frontiers in Econometrics, Academic Press, New York, 1974, pp. 105–142. https://eml.berkeley.edu/reprints/mcfadden/zarembka.pdf
• A. Albert and J. A. Anderson, On the existence of maximum likelihood estimates in logistic regression models, Biometrika 71(1), 1984, pp. 1–10. https://doi.org/10.1093/biomet/71.1.1
• S. J. Haberman, The Analysis of Frequency Data, University of Chicago Press, 1974.
• J. Berkson, Maximum likelihood and minimum χ² estimates of the logistic function, Journal of the American Statistical Association 50, 1955, pp. 130–162. https://doi.org/10.1080/01621459.1955.10501255

Motivation

The conditional logit model assigns to an alternative $i$ in a finite choice set the probability $e^{V_i}/\sum_j e^{V_j}$, where $V_j$ is a "representative utility" built from observed attributes of the alternative and the decision maker. It is the workhorse of discrete choice econometrics, transportation demand forecasting, marketing and revenue management, where it underlies multinomial logit assortment and pricing models. Its appeal for applied work is computational; its appeal for economics is that it can be read as the aggregate behaviour of a population of utility maximizers. This mission formalizes the result that makes that reading exact: Lemmas 1 and 2 of D. McFadden, Conditional logit analysis of qualitative choice behavior (1974), which show that, under a mild regularity condition, logit choice probabilities arise from random utility maximization exactly when the idiosyncratic taste shocks follow the extreme value (Gumbel) distribution.

Timeline:

• 1959. J. Marschak gives a nonconstructive proof that i.i.d. extreme value shocks yield logit probabilities; R. D. Luce's choice axiom appears the same year.
• 1965. Luce and Suppes publish the constructive argument, attributed to E. Holman and A. Marley, that is reproduced as the proof of Lemma 1.
• 1974. McFadden proves the converse (Lemma 2): if i.i.d. shocks with a translation complete distribution produce logit probabilities, the distribution is extreme value.
• Later. The random utility characterization was extended to correlated shocks (generalized extreme value models, McFadden 1978), which are not part of this mission.

Setting

An individual faces $J \ge 1$ alternatives with representative utilities $V_1, \dots, V_J \in \mathbb{R}$. The utility of alternative $j$ is $U_j = V_j + \varepsilon_j$, where the taste shocks $\varepsilon_1, \dots, \varepsilon_J$ are independent and identically distributed with a common law $\mu$ on $\mathbb{R}$ and distribution function $G(t) = \mu((-\infty, t])$. The individual chooses the alternative of highest utility, so the selection probability of $i$ is (Equation (2) of the paper)

$$P_i(V) = \Pr\big[\varepsilon_j - \varepsilon_i < V_i - V_j \ \text{ for all } j \ne i\big],$$

computed under the product law of the shocks. The logit formula (Equation (12)) is $L_i(V) = e^{V_i}/\sum_{j=1}^J e^{V_j}$. The extreme value law (Equation (13)) is $G(\varepsilon) = e^{-e^{-\varepsilon}}$.

A law $\mu$ is translation complete if for every function $h$ of bounded total variation on $\mathbb{R}$ with $h(\pm\infty) = 0$, the condition $\int h(e + a)\, d\mu(e) = 0$ for every real $a$ forces $h = 0$ outside a Lebesgue-null set. Laws whose characteristic function never vanishes, the extreme value law among them, are translation complete (footnote 5 of the paper).

In Lean the law is μ : Measure ℝ with [IsProbabilityMeasure μ], $G$ is ProbabilityTheory.cdf μ, the selection probability is selProb μ V i for V : Fin J → ℝ, and the logit formula is logitProb V i.

Formalization targets

Goal: the characterization

Fix a universe $X$ of alternatives with a representative utility map $u:X\to\mathbb R$ onto the real line. For a translation complete law $\mu$ normalized by $G(0) = e^{-1}$,

$$\Big(\text{for every finite }B\subseteq X,\ i\in B:\ P_i(B) = \frac{e^{u(i)}}{\sum_{j\in B} e^{u(j)}}\Big) \iff \Big(\forall \varepsilon \in \mathbb{R}:\ G(\varepsilon) = e^{-e^{-\varepsilon}}\Big).$$

The normalization only fixes the location of the shocks: without it the conclusion is the one-parameter family of Lemma 2 below.

Milestones

1. Equation (3) for i.i.d. shocks without atoms: $P_i(V) = \int \prod_{j \ne i} G(\varepsilon + V_i - V_j)\, dG(\varepsilon)$.
2. The integrand of Lemma 1's proof: under (13), the density times the other distribution functions equals $e^{-\varepsilon} \exp\big(-e^{-\varepsilon} \sum_j e^{V_j - V_i}\big)$.
3. Lemma 1: extreme value shocks give $P_i(V) = L_i(V)$ for every $J$ and $V$.
4. The functional equation of Lemma 2's proof: $G(v - \log K) = G(v)^K$ for every positive integer $K$ and real $v$.
5. Values at logarithms of rationals: with $\alpha = -\log G(0)$, $\alpha > 0$ and $G(\log(K/L)) = e^{-\alpha L / K}$ for positive integers $K, L$.
6. Lemma 2: $G(\varepsilon) = e^{-\alpha e^{-\varepsilon}}$ for some $\alpha > 0$, and $G(0) = e^{-1}$ gives (13).

Significance

The result separates two readings of the logit formula. Lemma 1 shows that it is consistent with utility maximization; Lemma 2 shows that, within the class of i.i.d. additive random utility models with translation complete shocks, the extreme value law is the only one consistent with it. Consequences drawn from the random utility reading, such as the log-sum formula for expected maximum utility used in welfare analysis, therefore apply to logit models without further distributional assumptions inside that class. The same reading supports the interpretation of multinomial logit demand in assortment optimization and revenue management.

Both lemmas are proved in the paper and in later textbooks; neither is open. As far as a search of the Prove2Me catalog shows, neither has a machine-checked proof. The mission provides Lean statements of the random utility model with i.i.d. shocks, of translation completeness and of the extreme value law that later discrete choice formalizations can reuse.

Difficulty

Lemma 1 is a computation with the extreme value density; its formal cost lies in passing from the product-measure probability (2) to the iterated integral (3) and evaluating an improper integral. Lemma 2 is harder. The natural first idea is to differentiate the logit identity in the utilities and solve a differential equation for $G$; this requires a density, which Lemma 2 does not assume. Without a density, the only handle on $G$ is the logit identity itself, an equality of integrals against $dG$ that holds for every utility vector; turning such integral identities into pointwise information about $G$ is where the hypothesis of translation completeness enters, and it yields statements only outside a Lebesgue-null set, so one-sided continuity of distribution functions is needed to recover identities at every point. A second subtlety is that the paper's (14) is written with $G$ while the event (2) is strict, so with a general law the integrals involve left limits of $G$.

Formalization scope

Alternatives are indexed by Fin J; the model is indexed by the utility vector, so the individual attributes $s$ and alternative attributes $x_j$ enter only through $V$. The shocks have joint law Measure.pi (fun _ => μ), which is what "independently identically distributed" means; a general joint law is not allowed. The event in (2) uses strict inequalities, and the selection probability is defined for every law, with no density. Translation completeness quantifies over BoundedVariationOn h Set.univ with limits $0$ at both ends; such $h$ are bounded and measurable, so the integrals are genuine, and "measure zero" is Lebesgue measure.

Lemma 2's hypothesis is stated on every finite subset of the paper's alternative universe, with a surjective utility map. Distinct alternatives may have the same utility. The printed proof uses $K$ equal-utility alternatives, which surjectivity alone need not supply; proving the stated theorem requires an additional continuity argument. A trivializing formalization is ruled out: the selection probability is a genuine product-measure probability, and the hypotheses of the goal are met by the extreme value law, which is translation complete with $G(0) = e^{-1}$.

A complete development needs: Fubini for Measure.pi over Fin J split at one coordinate; the Gumbel density and the improper integral $\int e^{-\varepsilon} e^{-c e^{-\varepsilon}} d\varepsilon = 1/c$; the facts that bounded-variation functions are bounded and measurable and that distribution functions are right-continuous with left limits. The integral representation (milestone 1) and the Gumbel computations are reusable in any random utility formalization. Proofs of the milestones, of the footnote-5 fact that the extreme value law is translation complete, and alternative proofs of Lemma 2 are welcome.

Selected references

• D. McFadden, Conditional logit analysis of qualitative choice behavior, in P. Zarembka (ed.), Frontiers in Econometrics, Academic Press, New York, 1974, pp. 105–142.
• J. Marschak, Binary choice constraints and random utility indicators, in K. Arrow, S. Karlin, P. Suppes (eds.), Mathematical Methods in the Social Sciences, Stanford University Press, 1960 (Stanford Symposium, 1959).
• R. D. Luce and P. Suppes, Preference, utility, and subjective probability, in R. D. Luce, R. Bush, E. Galanter (eds.), Handbook of Mathematical Psychology, Vol. III, Wiley, 1965.
• R. D. Luce, Individual Choice Behavior: A Theoretical Analysis, Wiley, 1959.
• W. Feller, An Introduction to Probability Theory and Its Applications, Vol. II, Wiley, 1966, p. 479.
• D. McFadden, Modelling the choice of residential location, in A. Karlqvist et al. (eds.), Spatial Interaction Theory and Planning Models, North-Holland, 1978, pp. 75–96.

Motivation

Linear programming, $\min\{c^\top x : Ax=b,\ x\ge0\}$ with $A\in\mathbb R^{d\times n}$, is the basic model of operations research, and the complexity of solving it is a central question of algorithm theory. Interior-point methods follow the central path: primal–dual pairs $(x,s)$ with $x,s>0$ and $x_is_i=t$ for every $i$, as the path parameter $t$ decreases to $0$. A classical short-step method needs $O(\sqrt n\log(n/\delta))$ iterations, each solving a linear system with the matrix $A\frac XSA^\top$, for a total of roughly $n^{2.5}$ operations or more.

Cohen, Lee and Song (J. ACM 68(1), 2021; arXiv:1810.07896) showed that linear programs can be solved in time $n^{\omega+o(1)}\log(n/\delta)$ (for the current values of the matrix multiplication exponent $\omega$ and its dual $\alpha$), matching the cost of multiplying two $n\times n$ matrices. The analysis has two halves: a data structure that maintains the projection matrix lazily, and the stochastic central path method, which replaces each Newton step by a sparse random step and proves that the iterates still stay close to the central path. This mission formalizes the second half.

Timeline: Karmarkar's projective method (1984) gave the first polynomial interior-point method; Renegar (1988) gave the $O(\sqrt n\log(1/\delta))$ path-following bound; Vaidya (1989) reduced the per-iteration cost with low-rank updates; Lee and Sidford (2014–2015) reduced the iteration count to $\widetilde O(\sqrt d)$; Cohen, Lee and Song (STOC 2019, J. ACM 2021) reached $n^\omega$; van den Brand (2020) derandomized the result.

Setting

Vectors are in $\mathbb R^n$ and products, quotients and roots of vectors are coordinatewise. For $\epsilon$ and vectors $a,b$, $a\approx_\epsilon b$ means $(1-\epsilon)b_i\le a_i\le(1+\epsilon)b_i$ for all $i$; $a\approx_\epsilon t$ for a scalar $t$ is defined likewise. The number of variables is $n\ge10$ and $A$ has full row rank $d\le n$.

The potential is $\Phi_\lambda(r)=\sum_{i=1}^n\cosh(\lambda r_i)$, evaluated at $r=\mu/t-1$ with $\mu=xs$; it is small exactly when every $x_is_i$ is close to $t$.

StochasticStep (Algorithm 1) takes positive $x,s$, a direction $\delta_\mu$, a sampling parameter $k$ and the output $\widetilde v$ of a data structure with $x/s\approx_{\epsilon_{\mathrm{mp}}}\widetilde v$. It rescales to $\overline x=x\sqrt{\widetilde v/w}$, $\overline s=s\sqrt{w/\widetilde v}$ ($w=x/s$), draws a sparse vector $\widetilde\delta_\mu$ with independent coordinates, $\widetilde\delta_{\mu,i}=\delta_{\mu,i}/p_i$ with probability $p_i=\min(1,k(\delta_{\mu,i}^2/\|\delta_\mu\|_2^2+1/n))$ and $0$ otherwise, and computes the step $(\widetilde\delta_x,\widetilde\delta_s)$ through the projection $\overline P=\sqrt{\overline X/\overline S}A^\top(A\frac{\overline X}{\overline S}A^\top)^{-1}A\sqrt{\overline X/\overline S}$. The draw is repeated until $\|\overline s^{-1}\widetilde\delta_s\|_\infty$ and $\|\overline x^{-1}\widetilde\delta_x\|_\infty$ are at most $1/(100\log n)$; the output is $(x+\widetilde\delta_x,s+\widetilde\delta_s)$.

Main (Algorithm 2) sets $\epsilon=\frac1{40000\log n}$, $\epsilon_{\mathrm{mp}}=\frac1{40000}$, $k=1000\epsilon\sqrt n\log^2n/\epsilon_{\mathrm{mp}}$, $\lambda=40\log n$, starts at $t=1$, and in each iteration sets $t^{\mathrm{new}}=(1-\frac{\epsilon}{3\sqrt n})t$, takes the direction

$$\delta_\mu=\Big(\frac{t^{\mathrm{new}}}{t}-1\Big)xs-\frac\epsilon2t^{\mathrm{new}}\frac{\nabla\Phi_\lambda(\mu/t-1)}{\|\nabla\Phi_\lambda(\mu/t-1)\|_2},$$

runs StochasticStep, and falls back to a deterministic ClassicalStep whenever $\Phi_\lambda(\mu^{\mathrm{new}}/t^{\mathrm{new}}-1)>n^3$.

Formalization targets

Goal: Lemma 4.14

For every iteration $j$, almost surely Assumption 4.1 holds for the input of iteration $j$ (in particular $x^js^j\approx_{0.1}t_j$ and $\|\delta_\mu\|_2\le\epsilon t_j$), almost surely the resampling loop of iteration $j$ succeeds with positive probability, and

$$\mathbb P(\text{ClassicalStep is used in iteration }j)\le\frac{10}{n^2}.$$

The paper writes $O(1/n^2)$; its proof gives the constant $10$.

Milestones

Lemma A.1 (variance of a product), Lemma 4.12 (properties of $\Phi_\lambda$), Lemma 4.2 (explicit step), Lemma 4.3 and Claim 4.7 (moments and success probability of the sampled step), Lemma 4.8 (moments of $\mu^{\mathrm{new}}$), and Lemma 4.13:

$$\mathbf E\Big[\Phi_\lambda\Big(\frac{\mu^{\mathrm{new}}}{t^{\mathrm{new}}}-1\Big)\Big]\le\Phi_\lambda\Big(\frac\mu t-1\Big)-\frac{\lambda\epsilon}{15\sqrt n}\Big(\Phi_\lambda\Big(\frac\mu t-1\Big)-10n\Big).$$

Significance

Lemma 4.14 is what makes the randomized method usable: the iterates stay in the $0.1$-neighbourhood of the central path along the whole run, and the expensive fallback is rare enough that its expected cost, $\widetilde O(n^{2.5})\cdot 10/n^2$, is negligible. The paper's cost bound (Lemma 4.16) and its main theorem rest on it. The same potential-based "stochastic central path" analysis was reused in later solvers, for instance for empirical risk minimization (Lee, Song and Zhang, COLT 2019).

The result is proved in the paper; no machine-checked version exists. A formalization pins down the probabilistic model that the paper leaves implicit (independence of the sampled coordinates, the law of the resampling loop, a data structure and fallback that see only the past) and checks the constants, several of which are tight against printed slack (Remark 4.4).

The running-time claims of the paper (Theorem 2.1's expected time $n^{\omega+o(1)}$, Lemma 4.16, Section 5) are not part of this mission: they live in an arithmetic cost model that Lean does not have. The accuracy guarantee of Theorem 2.1 (Lemma A.6, ClassicalStep from [57]) is also outside the mission.

Difficulty

The obvious argument would bound each quantity under the product law of the sparse direction. But StochasticStep resamples, so the step actually taken is distributed according to that law conditioned on a success event, and expectations and variances shift. A second difficulty is that $\Phi_\lambda$ is controlled only in expectation, while Assumption 4.1 must hold surely at every iteration; this is reconciled by the deterministic ClassicalStep fallback, which caps $\Phi_\lambda$ at $n^3$, and by an induction over iterations of $\mathbf E[\Phi]\le10n$ under the trajectory law. Claim 4.7 needs a Bernstein inequality, which Mathlib does not yet provide.

Formalization scope

Coordinates are Fin n, vectors Fin n → ℝ, $A$ a Matrix (Fin d) (Fin n) ℝ with A.rank = d, and $\log$ the natural logarithm. $\|\cdot\|_2$ is written out as $\sqrt{\sum_iv_i^2}$; $\|\cdot\|_\infty\le c$ is stated coordinatewise. The sampled direction has law Measure.pi of two-point laws; the step taken by StochasticStep has that law conditioned (ProbabilityTheory.cond) on the success event, and every $\mathbf E$, $\mathbf{Var}$ of Lemmas 4.3, 4.8 and 4.13 is under this conditioned law. mp.Query is replaced by its value $\overline P(\overline X\overline S)^{-1/2}\widetilde\delta_\mu$; the data structure and ClassicalStep are arbitrary measurable functions $U_j$, $C_j$ of the history with the only properties the paper uses. The trajectory is Mathlib's Ionescu-Tulcea measure, with kernels equal to the step law of Main. $n$ is the number of variables of the program the loop runs on.

Deviations from the page, all recorded in the items: Assumption 4.1 is used with $\epsilon\le1/(40000\log n)$ instead of the printed $<$, because Main sets $\epsilon$ to exactly that value; $O(1/n^2)$ is instantiated as $10/n^2$, the constant of the paper's proof; the conclusions of Lemma 4.14 are stated for every iteration index rather than while $t>\delta^2/(32n^3)$; at $\nabla\Phi_\lambda=0$ the second term of $\delta_\mu$ is $0$. No hypothesis $k\le n$ is imposed.

A trivializing formalization is ruled out: every statement that integrates against the conditioned law also concludes that this law is a probability measure (so it cannot be the zero measure), the goal concludes that each resampling loop succeeds with positive probability, the oracles $U_j$, $C_j$ cannot see the coins of the current iteration, and the goal is about the whole iterated process from the initial point, not one step from an arbitrary law.

Contributions welcome: a Bernstein inequality for bounded independent sums, conditional-law lemmas for cond of Measure.pi, and Markov-kernel measurability for the step law; these are reusable beyond this mission.

Selected references

• M. B. Cohen, Y. T. Lee, Z. Song, Solving Linear Programs in the Current Matrix Multiplication Time, J. ACM 68(1), Article 3, 2021. https://doi.org/10.1145/3424305 (arXiv:1810.07896, https://arxiv.org/abs/1810.07896)
• N. Karmarkar, A new polynomial-time algorithm for linear programming, Combinatorica 4, 1984. https://doi.org/10.1007/BF02579150
• J. Renegar, A polynomial-time algorithm, based on Newton's method, for linear programming, Math. Programming 40, 1988. https://doi.org/10.1007/BF01580724
• P. M. Vaidya, Speeding-up linear programming using fast matrix multiplication, Proc. 30th FOCS, 1989.
• Y. T. Lee, A. Sidford, Path finding methods for linear programming, FOCS 2014. https://doi.org/10.1109/FOCS.2014.52
• Y. T. Lee, Z. Song, Q. Zhang, Solving Empirical Risk Minimization in the Current Matrix Multiplication Time, COLT 2019. https://arxiv.org/abs/1905.04447
• J. van den Brand, A deterministic linear program solver in current matrix multiplication time, SODA 2020. https://doi.org/10.1137/1.9781611975994.16