Operations Research

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

Why partial balance

Equilibrium distributions of Markov models of networks are rarely computed by solving the full equilibrium equations directly. In the classical product-form results (migration processes, Jackson and Kelly networks, loss networks, clustering processes) the equilibrium distribution satisfies a stronger, local family of equations, and that is what makes it computable. The strongest such family is detailed balance, which characterizes reversibility. Many models that are not reversible still satisfy an intermediate family, partial balance: the probability flux balances not pair by pair but across a chosen set of transitions. F. P. Kelly's Reversibility and Stochastic Networks (Wiley, 1979) uses partial balance throughout: quasi-reversibility in Chapter 3 is a form of it, and the models that display it tend to be insensitive, meaning their equilibrium distribution does not change when exponential holding times are replaced by general ones with the same mean.

Section 9.4 of the book collects what partial balance means in a single statement. Theorem 9.5 summarizes Exercises 1.6.2–1.6.4 and 1.7.7–1.7.8 and gives five operational characterizations of partial balance. Corollaries 9.6–9.8 translate them into the language of spatial processes, and Theorem 9.9 turns partial balance of a coarse description into a product-form equilibrium for a finer one. This last step is the mechanism behind insensitivity.

Setting

A Markov process on a state space $\mathcal S$ has transition rates $q(j,k)\ge0$ for $j\ne k$, with $q(j,j)=0$. It is irreducible: every state can be reached from every other through transitions of positive rate. An equilibrium distribution is a collection of positive numbers $\pi(j)$ summing to one that satisfies the equilibrium equations

$$\pi(j)\sum_{k\in\mathcal S}q(j,k)=\sum_{k\in\mathcal S}\pi(k)q(k,j),\qquad j\in\mathcal S.$$

For an irreducible process on a finite state space it exists and is unique.

Given a set $\mathcal A\subseteq\mathcal S$, $\pi$ satisfies partial balance with respect to $\mathcal A$ if

$$\pi(j)\sum_{k\in\mathcal A}q(j,k)=\sum_{k\in\mathcal A}\pi(k)q(k,j),\qquad j\in\mathcal A.$$

Truncating the process to $\mathcal A$ deletes every transition out of $\mathcal A$, and the result is required to be irreducible within $\mathcal A$. The time-reversed process of a process with equilibrium distribution $\pi$ has rates $\pi(k)q(k,j)/\pi(j)$.

A spatial process has $J$ sites, the vertices of a graph $G$. Site $j$ carries an attribute $n_j$ from a finite set $\mathcal N_j$, and the state space is $\mathcal S=\mathcal N_1\times\cdots\times\mathcal N_J$. Write $T_j^m\mathbf n$ for the state $\mathbf n$ with the attribute of site $j$ changed to $m$. The process must satisfy three conditions: only one site changes at a time; the rate $q(\mathbf n,T_j^m\mathbf n)$ depends on the other sites only through the neighbours of $j$; and any $T_j^m\mathbf n$ can be reached from $\mathbf n$ without changing the other sites. The conditional distribution of site $j$ given the rest is $P(n_j\mid\mathbf n_{G-j})=\pi(\mathbf n)/\sum_m\pi(T_j^m\mathbf n)$. A Markov field is a positive distribution whose conditional distributions depend only on the neighbours of $j$.

Formalization targets

Goal: Theorem 9.5

For an irreducible process on a finite $\mathcal S$ with equilibrium distribution $\pi$, a nonempty $\mathcal A$ within which the truncated process is irreducible, and a constant $c>0$, $c\ne1$, the following are equivalent:

1. partial balance with respect to $\mathcal A$;
2. the equilibrium distribution of the truncated process is $\pi(j)/\sum_{k\in\mathcal A}\pi(k)$;
3. multiplying the rates $q(j,k)$, $j,k\in\mathcal A$, by $c$ leaves the equilibrium distribution unchanged;
4. multiplying the rates $q(j,k)$, $j\in\mathcal A$, $k\notin\mathcal A$, by $c$ changes the equilibrium distribution to

$$B\pi(j)\ (j\in\mathcal A),\qquad Bc\pi(j)\ (j\notin\mathcal A),\qquad B^{-1}=\sum_{j\in\mathcal A}\pi(j)+c\sum_{j\notin\mathcal A}\pi(j);$$

5. time reversal and truncation to $\mathcal A$ commute.

When $\mathcal S-\mathcal A$ is nonempty, these are also equivalent to:

6. the chain observed just before each exit from $\mathcal A$ and the chain observed just after each entry into $\mathcal A$ have the same equilibrium distribution.

Milestones

• Corollary 9.6: the equivalences for the sets on which all sites but $j$ are frozen, i.e. partial balance (9.26) at a site.
• Corollary 9.7: on a state space with at least two states, partial balance at every site makes $\pi$ a Markov field, with $0<\pi(\mathbf n)<1$ as in the book's definition.
• Corollary 9.8: the equivalences for the set on which site $j$ is frozen at one attribute (9.27).
• Theorem 9.9: if a reduced description $\mathbf n=f(\mathbf x)$ has a distribution $\pi(\mathbf n)$ in partial balance for the rates (9.31), then the finer process has equilibrium distribution $\pi(\mathbf x)=\pi(\mathbf n)\prod_jP_j(x_j\mid n_j)$.

What the results give

Theorem 9.5 makes partial balance testable by operations on the process itself: truncation, speeding up or slowing down transitions, and time reversal. The book points to close relationships between statement (iv) and the product form of Section 2.3, and between statement (ii) and part (iii) of Theorem 3.12. Corollary 9.6 (iv) explains why the reversed migration process has such a simple form. Corollary 9.7 strengthens Theorem 9.3 by replacing reversibility with partial balance at each site. Theorem 9.9 is the step from partial balance to insensitivity: it is what the book uses to show that a spatial process keeps its equilibrium distribution when the lifetimes of attributes are mixtures of gamma distributions.

All of these results were proved in 1979. None has a machine-checked proof. The platform has the reversible special case of statement (ii) (KellyStochasticNetworks.truncated_reversible) and the rate-level objects for time reversal and truncation, which this mission reuses. Formalizing Theorem 9.5 also produces a reusable account of embedded exit and entry chains of a finite Markov process.

Difficulty

Most of the equivalences (i)–(v) are short manipulations of the equilibrium equations, but each direction from a property of an altered process back to partial balance needs uniqueness of equilibrium distributions for irreducible finite processes, and (v) ⇒ (i) also needs their existence. Mathlib has neither in the form needed here. Statement (vi) is a different kind of claim. The equilibrium distribution of the exit chain is proportional to the exit flux $\pi(j)\sum_{k\notin\mathcal A}q(j,k)$, and that of the entry chain to the entry flux. Proving this requires the hitting distributions and the Green's function of the jump chain killed on leaving a set, and the convergence of the series that define them. Corollary 9.8 (v) inherits that work. Theorem 9.9 requires uniqueness for the reduced frozen processes and careful bookkeeping of the fibres $\{x_j:f_j(x_j)=n_j\}$.

Formalization scope

State spaces are finite types, and every sum is an unconditional sum over a finite type. Because the state space is finite, the book's extra condition for (vi), a finite flux out of $\mathcal A$, holds automatically. Rates are real functions with $q(j,j)=0$ built into the hypotheses; the reduced rates of (9.31) also have zero self-rates. "The equilibrium distribution of a process is $X$" means: $X$ is positive, sums to one, satisfies the equilibrium equations, and every distribution with these properties equals $X$. The book's "$c\ne0$ or $1$" is read as $c>0$, $c\ne1$, so that altered rates remain rates. Truncation to $\mathcal A$ is the published truncatedRates, a process on the subtype $\mathcal A$, and the reversed rates are the published reversedRates. The exit and entry chains are defined from the jump chain through series of restricted matrix powers. Spatial processes live on $\prod_j\mathcal N_j$ with $T_j^m$ given by Function.update. In Theorem 9.9 the graph is complete, as the book assumes from p. 202 on.

The equilibrium distribution of the truncated process must be the unique positive normalized solution of the truncated equilibrium equations. It must not be defined as the conditional distribution, which would make (ii) a tautology. For the same reason (iii) and (iv) are stated through the equilibrium equations of the altered rates, not by assumption.

Theorem 9.10 (p. 207) is not a target. Its hypothesis, that a nominal lifetime "can have any distribution with unit mean", is not defined on the page, and the point-process and lifetime description it needs lies outside this rate-level development.

Contributions are welcome at every level: uniqueness and existence of equilibrium distributions for irreducible finite rate matrices (reusable across the series), convergence of the killed Green's function, and proofs of the corollaries from the goal.

Selected references

• F. P. Kelly, Reversibility and Stochastic Networks, Wiley, Chichester, 1979; reissued Cambridge University Press, 2011. https://doi.org/10.1017/CBO9781139171724 (§9.4, pp. 200–208; §1.6, pp. 25–27)
• F. P. Kelly and E. Yudovina, Stochastic Networks, Cambridge University Press, 2014. https://doi.org/10.1017/CBO9781139565363

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