Probability

Why a fluid limit for matching markets

Two-sided platforms (ride hailing, delivery, freelancing, logistics) collect demand and supply requests over a matching cycle and then clear them with an optimization problem. When such a platform runs an experiment, treated and control units compete for the same supply, so the outcome of one unit depends on the treatment of others. Interference bias of this kind is the subject of Bright, Delarue and Lobel, Reducing Marketplace Interference Bias Via Shadow Prices (arXiv:2205.02274v4; Management Science, 2025). Their analysis of estimators, including the shadow-price estimator built from linear-programming duals, is carried out in a deterministic fluid limit. Theorem 1 of the paper justifies that step: as the market becomes dense, the expected matching value per unit of density converges to the optimal value of a deterministic linear program with the arrival rates as inputs.

The same kind of statement, the expected value of a random linear program converging to the value of the program at the mean data, is also used in network revenue management and in stochastic matching, where fluid linear programs serve as benchmarks.

Setting

Network. There are $n_d$ demand types and $n_s$ supply types. A directed network has one demand node $u^d_i$ per demand type (a sink), one supply node $u^s_j$ per supply type (a source) and a finite set $\mathcal U_*$ of intermediate nodes. Every arc $(u,w)$ carries a value $v_{u,w}\in\mathbb R$ and a capacity $k_{u,w}\ge 0$.

Matching program. For a density $\tau$, demand vector $d$ and supply vector $s$, $\Phi_\tau(d,s)$ is the maximum of $\sum_{(u,w)} v_{u,w}x_{u,w}$ over flows $x$ with inflow at most $d_i$ at $u^d_i$, outflow at most $s_j$ at $u^s_j$, flow balance at every intermediate node and $0\le x_{u,w}\le\tau k_{u,w}$ (display (2), p. 7). Write $\Phi:=\Phi_1$. Without intermediate nodes this is a capacitated bipartite transportation problem.

Random counts. Demand units of type $i$ arrive at rate $\lambda_i\ge 0$ and supply units of type $j$ at rate $\pi_j\ge 0$. At density $\tau$ the numbers of units collected in one cycle are independent Poisson counts

$$D^{\tau,\lambda}_i\sim\mathrm{Poisson}(\lambda_i\tau),\qquad S^\tau_j\sim\mathrm{Poisson}(\pi_j\tau).$$

Treatment effect. If a treatment changes the demand rates from $\lambda$ to $\lambda+\beta$ (with $\lambda+\beta\ge 0$), the global treatment effect at density $\tau$ is $\Delta^\tau=\frac1\tau\mathbb E[\Phi_\tau(D^{\tau,\lambda+\beta},S^\tau)-\Phi_\tau(D^{\tau,\lambda},S^\tau)]$ (Definition 1, p. 8).

Formalization targets

Goal: Theorem 1, display (7), p. 12

$$\lim_{\tau\to\infty}\mathbb E\Big[\frac1\tau\,\Phi_\tau\big(D^{\tau,\lambda},S^\tau\big)\Big]=\Phi(\lambda,\pi),$$

for every network, every $\lambda\ge 0$ and every $\pi\ge 0$, along integer densities $\tau$. No rate of convergence and no constant is fixed.

Milestones

1. Lemma 1 (p. 37), scaling. For integer $d,s\ge 0$ and $\tau>0$, $\frac1\tau\Phi_\tau(d,s)=\Phi(\frac1\tau d,\frac1\tau s)$, and $x$ is optimal on the right if and only if $\tau x$ is optimal on the left.
2. Poisson tail bound (proof of Lemma 2, p. 38). For $V\sim\mathrm{Poisson}(\lambda_0)$, $\mathbb P(V\ge t)\le e^{-t}$ for every $t\ge 7\lambda_0$ (Janson, Łuczak and Ruciński 2000, Corollary 2.4).
3. Lemma 2 (p. 38). If $X_n\sim\mathrm{Poisson}(\lambda n)$ for $n\ge 1$, then $Z=\sup_{n\ge1}X_n/n$ has finite mean, whatever the joint law of the $X_n$.

Companion: Corollary 1, display (8), p. 13

$$\lim_{\tau\to\infty}\Delta^\tau=\Phi(\lambda+\beta,\pi)-\Phi(\lambda,\pi).$$

Significance

Theorem 1 is the bridge between the paper's stochastic model and its deterministic analysis. Corollary 1 follows from it and identifies the estimand $\Delta=\Phi(\lambda+\beta,\pi)-\Phi(\lambda,\pi)$. Every fluid-limit result of the paper is stated about this deterministic quantity: that the RCT estimator overestimates the treatment effect (Theorem 2), and that in a symmetric experiment the shadow-price estimator is no more biased than the RCT estimator (Theorem 4). Without Theorem 1 those results describe a linear program, not the experiment.

The result is proved in the paper (Appendix C.1, pp. 37–39). As far as a search of the Prove2Me catalogue shows, it has not been formalized: the catalogue has fluid limits of queueing networks and several Poisson concentration bounds, but no law of large numbers for the value of a random network-flow program. A formal development would provide a reusable definition of the capacitated network-flow program with supply and demand bounds, its scaling identity, and an integrability tool for suprema of scaled Poisson counts. The two other missions of this series, on the shadow-price estimator in the fluid limit and in finite samples, use the same network program.

Difficulty

The obvious argument is: $\frac1\tau D^\tau\to\lambda$ and $\frac1\tau S^\tau\to\pi$, $\Phi$ is continuous, so the expectations converge. Each step needs care.

• The counts at different densities are defined only through their laws. A pathwise limit needs a coupling of the counts across $\tau$, which is not part of the statement.
• Continuity of $\Phi$ in its right-hand side, including at data where some $\lambda_i$ or $\pi_j$ is zero, is a statement about the value function of a parametric linear program. Mathlib has no theory of LP value functions.
• Passing the limit inside the expectation needs a dominating function. The paper builds one from the supremum of scaled Poisson counts (Lemma 2), whose integrability rests on the Poisson tail bound. That bound is not in Mathlib either.

Formalization scope

Representation.

• Demand types are Fin nd and supply types Fin ns; Lean index $i$ is the paper's type $i+1$.
• A network is a structure with arc type E and intermediate-node type I. Each arc's tail is a supply or intermediate node and its head is a demand or intermediate node, which builds in the source/sink structure. Values are real and capacities nonnegative.
• $\Phi_\tau(d,s)$ is the supremum of the objective over the feasible set of (2). For nonnegative data this set contains $0$ and lies in the box $[0,\tau k]$, so the supremum is a maximum.

Probability.

• The count law at density $\tau$ is the product of independent Poisson laws on $\mathbb N^{n_d}\times\mathbb N^{n_s}$. Independence across types and between demand and supply is the paper's reading of a Poisson process per type; it is used by the proofs but not written out.
• Expectations are Bochner integrals. The integrand of the goal is bounded by $\sum|v_{u,w}|k_{u,w}$, and every function on the countable count space is measurable, so the expectation is genuine and not the junk value $0$.
• The limit is taken along integer $\tau\to\infty$, as on p. 12.
• In Lemma 2, $\max_{n\ge1}$ is read as a supremum in $[0,\infty]$, and "finite mean" as a finite lower Lebesgue integral.

Ruled out. A trivializing formalization is excluded: the expectation is taken against the Poisson laws with means $\lambda_i\tau$ and $\pi_j\tau$, with the factor $1/\tau$ inside the expectation, not against point masses at the rates or in the rescaled form $\mathbb E[\Phi(\frac1\tau D^\tau,\frac1\tau S^\tau)]$. No hypothesis beyond $\lambda,\pi\ge 0$ is imposed.

Not stated.

• The second sentence of Theorem 1, the almost-sure convergence of optimal primal and dual solutions. It assumes unique optima only inside its proof and needs a coupling across $\tau$ that the paper does not construct.
• Lemma 1's "$z^*=\tau x^*$" presupposes unique optima. It is stated as the equivalence "$x$ is optimal for $\Phi(\frac1\tau d,\frac1\tau s)$ if and only if $\tau x$ is optimal for $\Phi_\tau(d,s)$".

The statements follow arXiv:2205.02274v4. Its printed slips are corrected as disclosed in each item's Formalization Note.

Infrastructure needed.

• Monotonicity and continuity of the value of a capacitated linear program in its right-hand side.
• The Poisson Chernoff bound.
• Convergence in probability of scaled Poisson counts.

The network-program definition and the Poisson tail bound are reusable beyond this mission. Contributions toward any of these are welcome.

Selected references

• I. Bright, A. Delarue, I. Lobel, Reducing Marketplace Interference Bias Via Shadow Prices, arXiv preprint, v4, 2024; Management Science, 2025. https://arxiv.org/abs/2205.02274v4
• S. Janson, T. Łuczak, A. Ruciński, Random Graphs, Wiley, 2000, Corollary 2.4. https://doi.org/10.1002/9781118032718

Motivation

Global sensitivity analysis asks how much of the uncertainty in the output $Y$ of a model, a simulator or a learned predictor is attributable to an input $X$. A standard way to answer is to compare the law of $Y$ with its law after $X$ is revealed, and to average that discrepancy over $X$. Variance-based indices (Sobol' 1993; Wagner 1995) compare conditional and unconditional means. Density-based and distribution-based indices compare whole distributions. Rényi's postulates for measures of dependence ask, among other things, that such an index vanish exactly under independence (zero-independence) and reach its maximum exactly when $Y$ is a function of $X$ (max-functionality, Postulate E). Variance-based indices fail the first postulate, because the conditional mean can be constant while the conditional law moves.

Borgonovo, Figalli, Plischke and Savaré (Management Science, 2025) propose indices in which the discrepancy is an optimal transport cost. The squared 2-Wasserstein case is the numerator of the Wasserstein correlation coefficient of Wiesel (2022), whose Theorem 2.2 gives max-functionality for that coefficient. The paper extends max-functionality to general continuous costs, through a strict-convexity property of the transport cost on Dirac masses (Lemma 4). This mission formalizes that result, Theorem 5, together with the zero-independence property, Proposition 3.

Setting

Let $(\Omega,\mathcal B,P)$ be a probability space and $\mathsf X$, $\mathsf Y$ Polish spaces with their Borel $\sigma$-algebras. The input $X:\Omega\to\mathsf X$ and the output $Y:\Omega\to\mathsf Y$ are measurable, with laws $\mu$ and $\nu$.

A cost is a continuous function $k:\mathsf Y\times\mathsf Y\to[0,+\infty]$ with $k(y,y')=0$ if and only if $y=y'$. It need not be symmetric, a metric, or finite. For probability measures $\nu,\nu'$ on $\mathsf Y$, a coupling is a probability measure $\pi$ on $\mathsf Y\times\mathsf Y$ with marginals $\nu$ and $\nu'$, and the optimal transport cost is

$$K(\nu,\nu')=\inf_{\pi\in\Pi(\nu,\nu')}\iint k(y,y')\,d\pi(y,y').$$

For a sub-$\sigma$-algebra $\mathcal F$ of the Borel sets of $\mathsf X$, a version of the conditional distribution of $Y$ given $(X,\mathcal F)$ is a Markov kernel $x\mapsto\nu^{\mathcal F}_x$ that is $\mathcal F$-measurable in $x$ and satisfies $P(X\in A,\,Y\in B)=\int_A\nu^{\mathcal F}_x(B)\,d\mu(x)$ for $A\in\mathcal F$ and Borel $B$. The OT-based sensitivity measure of Definition 1 is

$$\xi^K(Y,X\mid\mathcal F)=\int_{\mathsf X}K(\nu,\nu^{\mathcal F}_x)\,d\mu(x),$$

written $\xi^K(Y,X)$ when $\mathcal F$ is the full Borel $\sigma$-algebra. The candidate maximum is

$$\mathbb M^K[Y]=\mathbb E[k(Y,Y')]=\iint k(y,y')\,d\nu(y)\,d\nu(y'),$$

with $Y'$ an independent copy of $Y$; for $k(y,y')=\|y-y'\|_2^2$ it equals twice the total variance of $Y$.

Formalization targets

Goal: Theorem 5 (max-functionality)

If $\mathbb M^K[Y]<+\infty$, then for every sub-$\sigma$-algebra $\mathcal F$ and every version of the conditional distribution,

$$\xi^K(Y,X\mid\mathcal F)\le\mathbb M^K[Y],$$

and

$$\xi^K(Y,X\mid\mathcal F)=\mathbb M^K[Y]\iff Y=f(X)\ P\text{-a.s. for some }\mathcal F\text{-measurable }f:\mathsf X\to\mathsf Y.$$

Milestones

The steps of the paper's proof, in order:

1. (57): $K(\nu,\nu')\le\iint k\,d\nu\,d\nu'$.
2. $K(\nu,\delta_{y'})=\int k(y,y')\,d\nu(y)$.
3. Lemma 4: for $y_1\ne y_2$ in $\operatorname{supp}\nu$ and $t\in(0,1)$,

$$K\big(\nu,(1-t)\delta_{y_1}+t\delta_{y_2}\big)<(1-t)K(\nu,\delta_{y_1})+tK(\nu,\delta_{y_2}).$$

4. (62): the disintegration $\int K(\nu,\delta_y)\,d\nu(y)=\int_{\mathsf X}\int K(\nu,\delta_y)\,d\nu^{\mathcal F}_x(y)\,d\mu(x)$.
5. The bound $K(\nu,\nu')\le\int K(\nu,\delta_y)\,d\nu'(y)$.
6. (63): at the maximum, that bound is an equality for $\mu$-a.e. $x$.
7. (65): Jensen's inequality for $K(\nu,\cdot)$ over the representation (64) of a mixture.

Companion: Proposition 3 (zero-independence)

$$\xi^K(Y,X)=0\iff X\text{ and }Y\text{ are independent}.$$

Significance

The result. Together, Proposition 3 and Theorem 5 place $\xi^K$ on a scale from $0$ (no information about $Y$ in $X$) to $\mathbb M^K[Y]$ (complete information), so the ratio $\xi^K(Y,X\mid\mathcal F)/\mathbb M^K[Y]$ is a normalized measure of association. This normalization is what lets practitioners rank inputs of a simulator or features of a model on a common scale. The result covers all continuous costs vanishing exactly on the diagonal, including the powers of a metric that give the Wasserstein indices, and arbitrary sub-$\sigma$-algebras $\mathcal F$, which model partial information about $X$ (for example the information carried by a function $g(X)$).

Formalizing it. The theorem is proved in the paper; the work here is to formalize the known proof. As far as is known, no machine-checked version of this result or of Wiesel's Wasserstein-correlation theorem exists. Along the way the development needs facts about optimal transport with extended-real costs that are not in Mathlib: the product-coupling bound, the transport cost to a Dirac mass, convexity and lower semicontinuity of $K(\nu,\cdot)$, and cyclical monotonicity of optimal plans.

Difficulty

The inequality and the "functional dependence ⇒ equality" direction are direct computations with the product coupling and with Dirac conditional laws. The converse is the hard direction. Equality of the integrals gives only the pointwise equality (63) almost everywhere, and the step from there to "$\nu^{\mathcal F}_x$ is a Dirac mass" needs strict convexity of $K(\nu,\cdot)$. Strict convexity fails for general measures: $K(\nu,\cdot)$ is affine along many segments. It holds only between Dirac masses at distinct points of $\operatorname{supp}\nu$ (Lemma 4), and even there it fails when the costs involved are infinite. A non-Dirac conditional law is not a two-point mixture, so the strict inequality must be carried over to general measures with an integrand that takes values in $[0,+\infty]$. The last step, building an $\mathcal F$-measurable $f$ from the conditional laws, is measure-theoretic bookkeeping that is easy to get wrong.

Formalization scope

• Spaces and values. $\mathsf X$ and $\mathsf Y$ are Polish with their Borel $\sigma$-algebras; the paper's supports in $\mathbb R^{n_X}$, $\mathbb R^{n_Y}$ are such spaces. All of $k$, $K$, $\mathbb M^K$ and $\xi^K$ take values in $[0,+\infty]$ (ℝ≥0∞), so $+\infty$ is allowed as on the page, and integrals are lower Lebesgue integrals.
• Transport cost. $K$ is the published definition RWPI.SqrtLasso.transportCost. Its first argument is the marginal $\nu$, and the cost is evaluated as $k(y,y')$ with $y\sim\nu$.
• Conditional laws. The conditional distribution is not a chosen object. Any Markov kernel satisfying the defining property IsCondDistrib is admissible, and every theorem quantifies over all versions and over all sub-$\sigma$-algebras $\mathcal F$.
• Ruled-out encodings. The goal must not fix $\mathcal F$ to the Borel $\sigma$-algebra or fix a version of the kernel. Its conclusion must not be replaced by "$\nu^{\mathcal F}_x$ is a Dirac mass for $\mu$-a.e. $x$", which is an intermediate step of the proof, not the theorem.
• Disclosed correction. Lemma 4 carries the extra hypotheses $K(\nu,\delta_{y_1}),K(\nu,\delta_{y_2})<+\infty$. As printed, both sides of (19) can be $+\infty$ under the lemma's assumptions. Theorem 5 is stated as printed.
• Not posed. The cyclical monotonicity of optimal plans, cited from Villani (2008, Ch. 5), and the existence of optimal plans used in Proposition 3 are standard optimal-transport facts. They are not milestones but are welcome as supporting lemmas.

Selected references

• E. Borgonovo, A. Figalli, E. Plischke, G. Savaré, Global Sensitivity Analysis via Optimal Transport, Management Science, 2025 (authors' manuscript of 7 March 2024 used here). https://doi.org/10.1287/mnsc.2023.01796
• J. Wiesel, Measuring Association with Wasserstein Distances, Bernoulli 28(4):2816–2832, 2022. https://doi.org/10.3150/21-BEJ1438
• C. Villani, Optimal Transport: Old and New, Springer, 2008. https://doi.org/10.1007/978-3-540-71050-9
• A. Rényi, On Measures of Dependence, Acta Mathematica Academiae Scientiarum Hungaricae 10:441–451, 1959. https://doi.org/10.1007/BF02024507
• I. M. Sobol', Sensitivity Estimates for Nonlinear Mathematical Models, Mathematical Modelling and Computational Experiments 1:407–414, 1993.
• H. M. Wagner, Global Sensitivity Analysis, Operations Research 43(6):948–969, 1995. https://doi.org/10.1287/opre.43.6.948

Motivation

Thompson sampling (TS) is the standard randomized heuristic for the Bayesian multi-armed bandit: at each decision it samples a parameter vector from the current posterior and pulls the arm that would be best if the sample were true. It ignores the remaining horizon. Min, Maglaras and Moallemi (arXiv:1902.04251, Management Science, 2025) obtain horizon-aware generalizations of TS from the information relaxation duality of Brown, Smith and Sun (Oper. Res. 2010). Each penalty function yields an upper bound on the Bayes-optimal value and an IRS policy. Such a policy samples a whole future outcome and plays the first action of the penalized clairvoyant optimum. With the zero-information penalty it is exactly TS.

Two of these policies, IRS.FH and IRS.V-Zero, are designed to interpolate between TS and the myopic rule. When one period remains they pull the arm with the highest predictive mean (Remark 5 of the paper). This mission formalizes the other end, Proposition 1: when many periods remain, they randomize exactly as Thompson sampling does. The paper presents this as the sense in which IRS policies "behave like TS during the initial decision epochs" and shift smoothly from exploration to exploitation without a tuning parameter.

Setting

There are $K \ge 1$ arms $\mathcal A = \{1,\dots,K\}$. Arm $a$ has an unknown parameter $\theta_a$. Given $\theta_a$, its rewards $R_{a,1}, R_{a,2}, \dots$ are i.i.d. with law $\mathcal R_a(\theta_a)$. The prior of $\theta_a$ is $\mathcal P_a(y_a)$, a member of a conjugate family indexed by a belief $y_a$. A Bayesian update $\mathcal U_a$ maps $(y_a, r)$ to the posterior belief after a reward $r$. The arms are independent. Throughout, rewards are integrable: $\mathbb E_{y_a}|R_{a,1}| < \infty$ for every belief (display (3)).

• The conditional mean reward is $\mu_a(\theta_a) = \mathbb E[R_{a,n} \mid \theta_a]$, and the predictive mean is $\bar\mu_a(y_a) = \mathbb E_{y_a}[\mu_a(\theta_a)]$.
• The outcome $\omega = (\theta, (R_{a,n})_{a,n})$ has law $\mathcal I(y)$.
• The posterior predictive mean after the first $n$ rewards of arm $a$ is $\hat\mu_{a,n}(\omega; y_a) = \mathbb E_{y_a}[\mu_a(\theta_a) \mid R_{a,1},\dots,R_{a,n}]$.

At remaining horizon $T$ and belief $y$, the three policies act as follows:

1. TS samples $\tilde\theta \sim \mathcal P(y)$ and plays $A^{TS}(y) = \operatorname{argmax}_a \mu_a(\tilde\theta_a)$.
2. IRS.FH samples $\tilde\omega \sim \mathcal I(y)$ and plays $A^{IRS.FH}(T,y) = \operatorname{argmax}_a \hat\mu_{a,T-1}(\tilde\omega)$, the arm with the best posterior mean after $T-1$ sampled rewards of every arm.
3. IRS.V-Zero samples $\tilde\omega \sim \mathcal I(y)$ and solves the allocation problem

$$\max_{n \in \mathcal N_T} \sum_{a=1}^K \sum_{i=1}^{n_a} \hat\mu_{a,i-1}(\tilde\omega), \qquad \mathcal N_T = \Big\{n \in \mathbb N_0^K : \sum_a n_a = T\Big\}.$$

It then plays the arm $A^{IRS.V\text{-}Zero}(T,y) = \operatorname{argmax}_a \tilde n^*_a$ to which the optimal allocation $\tilde n^*$ assigns the most pulls.

Each action is a random variable driven by the policy's own sample.

Formalization targets

Goal: Proposition 1 (p. 25)

Assume $\mu_i(\theta_i) \ne \mu_j(\theta_j)$ almost surely under $\mathcal P(y)$, for all $i \ne j$. Then for every arm $a$,

$$\lim_{T\to\infty} \mathbb P[A^{IRS.FH}(T,y) = a] = \mathbb P[A^{TS}(y) = a], \qquad \lim_{T\to\infty} \mathbb P[A^{IRS.V\text{-}Zero}(T,y) = a] = \mathbb P[A^{TS}(y) = a].$$

This holds for every choice of the tie-breaking rules in the argmaxes and in the allocation problem.

Milestones

1. Proposition 4 (p. 53): $\mu_a(\theta_a)$ and $\max_a \mu_a(\theta_a)$ are integrable under the prior.
2. Proposition 5 (p. 54): $\hat\mu_{a,n}$ and the sample mean $\frac1n\sum_{i\le n} R_{a,i}$ converge almost surely to $\mu_a(\theta_a)$.
3. Display (109) (p. 54): if a sequence of vectors converges to a vector with distinct entries, every argmax rule eventually selects the argmax of the limit.
4. The allocation argument (113)–(123) (pp. 54–56): if the payoff rows converge to distinct limits, then for large $T$ every optimal allocation gives its largest share only to the arm with the largest limit.
5. Display (55) alone: the IRS.FH half of the goal.

Significance

Proposition 1 identifies the long-horizon behaviour of the two IRS policies with that of Thompson sampling. Together with Remark 5 it pins down both ends of the horizon. It explains why their simulated regret curves track TS early and the myopic policy late (§5 of the paper). It also shows that the design choice "play the arm with the largest allocation" in IRS.V-Zero is the one that recovers TS. Other selection rules, such as any arm with a positive allocation, need not.

The paper gives a complete proof in Appendix D.1–D.2. A machine-checked version adds three things. It fixes the probabilistic model of a conjugate Bayesian bandit with infinite reward sequences in Lean, which the companion missions on the same paper also need. It makes the role of tie-breaking explicit. And it repairs a step of the printed proof: the bound $n^\circ(A^\circ_T) \ge \lfloor T/2 \rfloor$ in §D.2 holds only for two arms, and the general case needs $\lceil T/K \rceil$ and a threshold growing with $K$. No formalization of this result is known on Prove2Me or in Mathlib.

Difficulty

The IRS.FH half reduces to almost-sure convergence of posterior means, which is a martingale argument (Lévy's upward theorem). That argument needs $\mu_a(\theta_a)$ to be measurable with respect to the limit $\sigma$-algebra of the reward sequence, a fact supplied by the strong law of large numbers. Building the outcome law, a mixture over $\theta$ of infinite product measures, and identifying the conjugate update with conditional expectation are where the measure theory lies.

The IRS.V-Zero half does not follow from convergence of the posterior means alone. The selected arm is a function of an optimal allocation of a combinatorial problem whose size grows with $T$. The naive expectation that the arm with the best limit receives "most" pulls has to be proved by an exchange argument that handles the finitely many early terms, where the posterior means can still be badly ordered. The printed version of this argument is correct only for $K = 2$.

Formalization scope

• Arms are Fin K with NeZero K; the reward sequence of an arm is ℕ → ℝ, 0-based.
• One parameter type and one belief type serve all arms; reward laws and priors are Markov kernels.
• Conjugacy is stated as Bayes' rule on rectangles for every belief. The update is assumed jointly measurable, and the standing integrability (3) is assumed for every belief.
• The outcome law $\mathcal I(y)$ is a product over arms of mixtures $\int \delta_\theta \otimes \mathcal R_a(\theta)^{\otimes \mathbb N}\, \mathcal P_a(y_a)(d\theta)$.
• $\hat\mu_{a,n}$ is the pathwise conjugate version $\bar\mu_a(y_{a,n})$ of the conditional expectation (10), defined at every outcome so that it can sit inside an argmax.
• $\mathbb P[A = a]$ is the mass of $\{a\}$ under the pushforward of the policy's sampling law. All tie-breaking rules are required to be measurable; otherwise the pushforward would be the zero measure.
• The assumption of the proposition is required almost surely under the prior, not as injectivity of $\mu$.
• A formalization in which the selectors were non-measurable, the outcome law degenerate, or the distinctness hypothesis unsatisfiable would make the goal vacuous. The definitions rule out the first two, and a deterministic two-arm instance shows the hypotheses are satisfiable.

The development needs Lévy's upward theorem (MeasureTheory.tendsto_ae_condExp in Mathlib) and the strong law of large numbers (ProbabilityTheory.strong_law_ae). It also needs a lemma identifying iterated conjugate updates with the posterior given finitely many rewards, and the measurability of $\theta \mapsto \mathcal R_a(\theta)^{\otimes\mathbb N}$. These last two are reusable for any Bayesian bandit or sequential-estimation development. Contributions of the two deterministic milestones are self-contained and welcome first.

Selected references

• S. Min, C. Maglaras, C. C. Moallemi, Thompson Sampling with Information Relaxation Penalties, Management Science, 2025; preprint arXiv:1902.04251v2, 2021. https://arxiv.org/abs/1902.04251
• D. B. Brown, J. E. Smith, P. Sun, Information Relaxations and Duality in Stochastic Dynamic Programs, Operations Research 58(4), 2010. https://doi.org/10.1287/opre.1090.0796
• W. R. Thompson, On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples, Biometrika 25(3/4), 1933. https://doi.org/10.2307/2332286
• D. Russo, B. Van Roy, Learning to Optimize via Posterior Sampling, Mathematics of Operations Research 39(4), 2014. https://doi.org/10.1287/moor.2014.0650

Motivation

Multi-server queues model call centres, hospital wards, cloud server pools and any system in which $n$ identical servers work through one first-come-first-served (FCFS) line. Capacity planning for such systems needs bounds on the steady-state queue length and on the probability of delay (the chance that an arriving job finds every server busy) that are explicit, hold for every $n$, and stay informative as the load approaches capacity.

For a single server, Kingman's bound (Kingman 1962, ref. [89] of the paper) gives exactly this: the mean wait is at most a multiple of $1/(1-\rho)$ that depends only on the first two moments of the inter-arrival and service times. For $n$ servers with general distributions (the GI/GI/n queue), no comparable explicit bound was known.

Timeline.

• 1962: Kingman's single-server bound.
• 1981: Halfin and Whitt identify the heavy-traffic regime $\sqrt n(1-\rho)\to\beta$ in which the probability of delay is nondegenerate (Halfin and Whitt 1981).
• 2013: Gamarnik and Goldberg prove (for the GI/G/n queue) that the steady-state number waiting is stochastically dominated by the supremum of a process built from independent renewal processes (Gamarnik and Goldberg 2013). Their bound is not explicit and needs more than two service moments.
• 2017–2025: Li and Goldberg derive explicit bounds with $1/(1-\rho)$ scaling under $\mathbb E[A^2]<\infty$ and $\mathbb E[S^{2+\epsilon}]<\infty$ (arXiv:1706.04628, Mathematics of Operations Research 2025).

Setting

A GI/GI/n queue has $n\ge1$ identical servers and one FCFS line. Jobs arrive with i.i.d. inter-arrival times distributed as a nonnegative random variable $A$. Each job needs an i.i.d. service time distributed as a nonnegative random variable $S$, independent of the arrivals. An arriving job starts service at the first server to become free, in order of arrival.

Write $\mu_A=1/\mathbb E[A]$ and $\mu_S=1/\mathbb E[S]$ for the arrival and service rates, and $\rho=\mu_A/(n\mu_S)$ for the traffic intensity. The moment index is

$$r^*=\sup\{r:\mathbb E[S^r]<\infty\},$$

which may be $+\infty$. Normalised moments are $\mathbb E[(A\mu_A)^2]$ and $\mathbb E[(S\mu_S)^r]$.

$Q(t)$ is the number of jobs in the system (in service and waiting) at time $t$. "$Q(\infty)$ exists" means that, whatever the initial condition, $Q(t)$ converges in distribution as $t\to\infty$ to a random variable $Q(\infty)$. The steady-state number waiting is $L(\infty)=\max(0,Q(\infty)-n)$, and the steady-state probability of delay is $\mathbb P(Q(\infty)\ge n)$.

The milestones also use equilibrium renewal processes. For a nonnegative $X$ with finite mean, its equilibrium law $R(X)$ satisfies $\mathbb P(R>y)=\frac1{\mathbb E[X]}\int_y^\infty\mathbb P(X>z)\,dz$. An equilibrium renewal process with renewals distributed as $X$ has its first renewal interval distributed as $R(X)$ and the later ones distributed as $X$.

$A_e(t)$ counts the renewals in $[0,t]$ of such a process for $A$. $N_{e,1}(t),\dots,N_{e,n}(t)$ count those of $n$ independent such processes for $S$, also independent of $A_e$.

Formalization targets

Goal: Theorem 2 (pp. 9–10)

Assume

• $\mathbb E[A^2]<\infty$;
• $r^*>2.5$;
• $\mu_A<n\mu_S$;
• $Q(\infty)$ exists.

Then for all $x>0$ and all $r\in[2.5,r^*)$,

$$\mathbb P\Big(L(\infty)\ge\frac x{1-\rho}\Big)\le 2\cdot10^4\,\mathbb E[(S\mu_S)^2]\,(10^6)^r\Big(\mathbb E[(S\mu_S)^2]^{r-1}r^{2.5r}+r^{1.5r}\mathbb E[(S\mu_S)^r]\Big)x^{-\frac r2}+1.1\,e^{-.0225\,\mathbb E[(A\mu_A)^2]^{-1}x},$$ $$\mathbb P\big(Q(\infty)\ge n\big)\le 2\cdot10^4\,\mathbb E[(S\mu_S)^2]\,(10^7)^r\Big(\mathbb E[(S\mu_S)^2]^{r-1}r^{2.5r}+r^{1.5r}\mathbb E[(S\mu_S)^r]\Big)\big(n(1-\rho)^2\big)^{-\frac r2}+1.1\,e^{-.0028\,\mathbb E[(A\mu_A)^2]^{-1}n(1-\rho)^2}.$$

The paper states the infimum over $r$. With $r$ universally quantified the statement is equivalent, because the exponential term does not depend on $r$.

Milestones

They follow the paper's proof, in attack order:

• Lemma 1, the Gamarnik–Goldberg comparison $\mathbb P(L(\infty)\ge x)\le\mathbb P(\sup_{t\ge0}(A_e(t)-\sum_{i=1}^nN_{e,i}(t))\ge x)$;
• Lemma 2, which splits that supremum into an arrival part and a service part;
• Lemmas 8, 9 and 10, which reduce the arrival part to Kingman's bound with an explicit exponent;
• Lemma 3, the exponential bound for the arrival part;
• Lemma 4 and Theorem 10, conditional bounds from central moments of the pooled process $\sum_iN_{e,i}(t)$;
• Lemmas 5 and 6, those central-moment bounds;
• the explicit bound for $r\in(2.5,r^*)$ displayed on p. 28;
• Lemma 14, a transient comparison for a queue started with all servers busy;
• Lemma 7 and Corollary 5, comparisons with fewer servers that make the bound informative at the level $Q(\infty)\ge n$.

Significance

The result. The bounds hold for every $n$ and every load with $\rho<1$, and they depend on the distributions only through a few normalised moments. In the Halfin–Whitt regime, $n(1-\rho)^2$ stays bounded away from zero, so the delay bound stays informative. As $\rho\to1$ with $n$ fixed, the queue tail decays at the Kingman scale $1/(1-\rho)$. These are the first such bounds for general GI/GI/n queues; earlier explicit results needed Poisson arrivals or exponential service. The paper derives from Theorem 2 bounds at a fixed $r$ (Corollaries 1–2), a refinement beyond $1/(1-\rho)$ scaling when $n(1-\rho)^2$ is large (Corollary 3), and optimised bounds under growth conditions on $\mathbb E[(S\mu_S)^r]$ (Theorem 3).

Formalization. The theorem is proved in the paper, partly in a supplemental appendix, and none of it is formalized; there are no machine-checked results for multi-server queues with general distributions anywhere. A complete development would add to the platform the first FCFS multi-server queue with general distributions, its steady state, equilibrium renewal processes, and the Gamarnik–Goldberg comparison chain. Each is reusable for later work on many-server queues.

Difficulty

The obvious route is to compare the $n$-server queue with a single-server one, or to apply a heavy-traffic limit. The first loses a factor of $n$; the second gives limits but no bounds at a fixed $n$.

The paper's route needs two hard ingredients:

• Comparison with renewal processes. The stochastic comparison of Lemmas 1 and 14 is a theorem about the sample-path dynamics of FCFS service, not a calculation.
• Uniform moment bounds. Explicit control of $\mathbb E|\sum_{i=1}^nN_{e,i}(t)-nt|^r$ at the scale $(nt)^{r/2}$, uniformly in $n$ and $t$, is required. Classical renewal moment bounds are either asymptotic in $t$ or lose the dependence on $n$.

For the probability of delay, Lemma 1 at level $0$ gives only the trivial bound $1$. A second comparison across server counts (Lemma 7) is needed.

Formalization scope

Laws. $A$ and $S$ are represented by their laws $\mu_A,\mu_S$: probability measures on $\mathbb R$ carried by $[0,\infty)$, integrable, with positive mean. Positivity of the means is implicit in the paper, which needs $\mu_X=1/\mathbb E[X]$ finite; it is stated explicitly. $\mathbb E[(X\mu_X)^r]$ is a real integral and is used only where it is finite: $r<r^*$ for $S$, and $r=2$ with $\mathbb E[A^2]<\infty$ for $A$. $r^*$ is an extended real.

The queue.

• The queue is the Kiefer–Wolfowitz workload recursion, driven by an initial condition and i.i.d. pairs of inter-arrival and service times.
• An initial condition gives the first-arrival time and the residual service time at each server. No job waits at time $0$.
• "$Q(\infty)$ exists with law $\pi$" quantifies over every initial law of that kind, independent of the inputs. This class is narrower than the paper's, so the hypothesis is weaker and the stated theorem is stronger. The paper's proof uses only such initial conditions (Lemma 7 through Lemma 14).
• $L(\infty)$ is the truncated difference $\max(0,Q(\infty)-n)$.

Renewal processes.

• Counts are natural numbers, computed pathwise from the renewal intervals.
• Bounding families are arbitrary probability spaces carrying mutually independent intervals with the prescribed laws.
• Suprema of processes are taken in the extended reals.

Ruled out. None of the following trivializing formalizations is used:

• a steady-state hypothesis that no queue satisfies, or one restricted to an initial-condition class that drops what Lemma 7 needs;
• a real-valued supremum or infimum, which collapses to $0$ on unbounded families;
• a Bochner moment of a non-integrable function, which is $0$ in Lean.

The infimum over $r$ is stated as a universal quantifier.

Needed infrastructure. Equilibrium distributions and renewal counting processes, moment and maximal inequalities for sums of independent renewal processes, Kingman's exponential martingale bound, and the sample-path comparison of FCFS multi-server queues. Contributions on any of these are welcome: the renewal and random-walk results (Lemmas 5, 6, 8–10) are independent of queueing.

Selected references

• Y. Li and D. A. Goldberg, Simple and explicit bounds for multi-server queues with 1/(1−ρ) scaling, Mathematics of Operations Research, 2025; arXiv:1706.04628v3. https://arxiv.org/abs/1706.04628
• D. Gamarnik and D. A. Goldberg, Steady-state GI/G/n queue in the Halfin–Whitt regime, Annals of Applied Probability 23(6), 2382–2419, 2013. https://doi.org/10.1214/12-AAP905
• S. Halfin and W. Whitt, Heavy-traffic limits for queues with many exponential servers, Operations Research 29(3), 1981. https://doi.org/10.1287/opre.29.3.567
• J. F. C. Kingman, On queues in heavy traffic, Journal of the Royal Statistical Society B 24, 383–392, 1962 (ref. [89] of Li–Goldberg).
• J. F. C. Kingman, A martingale inequality in the theory of queues, Mathematical Proceedings of the Cambridge Philosophical Society 60(2), 359–361, 1964 (ref. [90] of Li–Goldberg; the source of Lemma 9).

Motivation

Many demand, traffic and click-through series are seasonal: their mean repeats after a fixed number of periods (a week of days, a year of months), while the length of that cycle is not known in advance. A retailer choosing among $K$ promotions, or a platform choosing among $K$ recommendations, then faces a multi-armed bandit whose arms have time-varying but periodic mean rewards. Standard stochastic bandit algorithms assume stationary means and fail here, and fully adversarial or arbitrarily drifting models give up the structure that periodicity offers.

N. Chen, C. Wang and L. Wang (arXiv:2005.08088, Operations Research, 2025) propose a two-stage policy for this problem: a spectral first stage that estimates every arm's period from a short exploration phase, and a confidence-bound second stage that learns and exploits the resulting "effective arms". Their main upper bound, Theorem 1, shows that the policy's regret grows like $\sqrt{T\sum_k T_k}$ up to logarithmic factors. This mission formalizes that theorem and the chain of lemmas behind it. The matching lower bound (Theorem 3 of the paper) is the subject of a separate mission of this series.

Setting

At each epoch $t=1,\dots,T$ a decision maker pulls one of $K$ arms. Pulling arm $k$ at epoch $t$ yields the reward $Y_{k,t}=\mu_{k,t}+\epsilon_t$. The mean $\mu_{k,t}$ has a minimum period $T_k\ge1$: $\mu_{k,t+T_k}=\mu_{k,t}$ for all $t\ge1$, and no smaller positive shift has this property. The standing assumptions are:

• Assumption 1: $\mu_{k,t}\in[0,1]$.
• Assumption 2: the noises $\epsilon_t$ are independent and sub-Gaussian with parameter $\sigma$, i.e. $\mathbb E[e^{\lambda\epsilon_t}]\le e^{\sigma^2\lambda^2/2}$ for all real $\lambda$ and $\mathbb P(|\epsilon_t|>x)\le2e^{-x^2/(2\sigma^2)}$.

The pseudo-regret of an arm sequence $\pi$ is $R^\pi_T=\sum_{t=1}^T(\max_k\mu_{k,t}-\mu_{\pi_t,t})$, and the number of effective arms is $d=\sum_{k=1}^KT_k$: arm $k$ at phase $t \bmod T_k$ is one effective arm.

Stage one (Algorithm 1). Each arm is pulled $n$ times in a row. From its $n$ rewards the policy computes the periodogram $|\tilde y_k(v)|$, the modulus of the discrete Fourier transform $\tilde y_k(v)=\frac1n\sum_sY_{k,s}e^{-2\pi ivs}$, on $v\in[0,1/2]$. A data-driven threshold $\tau_k$ is computed from the constants $U_1,U_2$ (sums of sinc-type maxima), $\bar\epsilon_v=\frac{2\sigma H}{1-\pi/24}\sqrt{\log(n)/n}$ and $\sup_v|\tilde y_k(v)|$. The algorithm repeatedly takes the highest peak above $\tau_k$ in $[g/n,1/2]$, snaps it to the nearest fraction $j_1/j_2$ with $j_2<n/(2g)$, and removes a $2g/n$-neighbourhood. The period estimate $\hat T_k$ is the least common multiple of the selected reciprocals.

Stage two (Algorithm 2). At every later epoch an $S$-round screening tournament, $S=\lfloor\log_2T\rfloor$, compares the arms at their current phases modulo $\hat T_k$. It uses estimated means $m^{(s)}_{k,t}$ and widths $w^{(s)}_{k,t}$ pooled from the stage-one samples and from round-specific index sets $\Psi^{(s)}(t)$. The weak oracle runs the same policy with the true periods in place of $\hat T_k$.

The parameters are fixed by (7): $n=\lfloor\sqrt{T/K}\rfloor$, $g=\lceil\sqrt n\rceil$, $H=\sqrt{1+\log n}$. Assumption 3 asks $T>4K$, $g\ge\max\{2,\sqrt n\}$ and $T_k<n/(2g)$. Assumption 4 asks the smallest present Fourier coefficient $b_k$ of each arm to dominate a combination of $\bar\epsilon_v$ and the largest coefficient $B_k$.

Formalization targets

Goal: Theorem 1

Under Assumptions 1–4 and (7), the two-stage policy satisfies

$$\mathbb E\big[R^\pi_T\big]\le C\sqrt{T\,d\,\log^2(T)\,\log(T/d)} ,$$

with $C$ depending only on $\sigma$ and $K$, and with Algorithm 2 run at $\delta=8/T$.

Stage one: correct period identification

• Lemma 1: uniform noise bound for the periodogram, $48n\exp(-(n\delta^2/4\sigma^2)(1-\pi/24)^2)$.
• Lemmas 2–3: leakage bounds $\pi B_kU_1$, $\pi B_kU_2$, and $B_k\le(\sup|\tilde\epsilon_k|+\sup|\tilde y_k|)/(1-\pi U_2)$.
• Lemmas 4–5: the threshold separates the neighbourhoods of the present frequencies from the rest, with probability $1-48/n^{H^2-1}$.
• Lemmas 6–7: the sharper noise bound near possible present frequencies, and the location of the peaks within $1/n$ of them.
• Theorem 2: $\mathbb P(\Lambda)\ge1-48K/n^{H^2-1}-200K/n^{0.867H^2-1}-200K/n^{0.6934H^2-1}$.

Stage two: the oracle's regret

• Lemma 15: an anytime sub-Gaussian deviation bound.
• Lemma 9: the confidence event $\mathcal E$ has probability at least $1-\delta S$.
• Lemma 10: on $\mathcal E$ the optimal arm survives screening, and each epoch's gap is at most $1$, $8\sigma/2^s$ or $2\sigma/\sqrt T$.
• Lemma 16: the summed widths of a round are at most $4\sigma\sqrt{|\Psi^{(s)}(T)|\,d\log(8T/\delta)\log(T/d)}$.
• Lemma 11: the oracle's expected regret is at most $C\sqrt{Td\log^2(T)\log(T/d)}$ with $C$ depending only on $\sigma$.

Significance

The theorem shows that unknown seasonality costs only the factor $\sqrt{d}=\sqrt{\sum_kT_k}$ that would be paid if every phase of every arm were a separate stationary arm. Up to logarithmic factors this matches the $\sqrt{TK}$ rate of the classical stationary bandit with $d$ arms in place of $K$. Together with the paper's lower bound $\Omega(\sqrt{T\max_kT_k})$ it locates the minimax regret of periodic bandits within a factor $\sqrt K$. Stage one is a non-asymptotic, finite-sample period-identification guarantee for a thresholded periodogram, which is of independent use in spectral estimation.

The paper proves these results; none of them is machine-checked. A complete development yields a verified finite-sample periodogram theory (noise bounds for random trigonometric polynomials, Dirichlet-kernel leakage estimates), a verified anytime sub-Gaussian inequality, and a verified analysis of a nested elimination algorithm with data reuse. The formalization also pins down several points the printed text leaves open: the constant $0.233$ of Lemma 6 and the exponent $0.694$ of Lemma 7 and Theorem 2 are slightly stronger than what the proofs give, and the printed claim that Theorem 1's constant is independent of $K$ is not established by the printed argument.

Difficulty

Two steps resist the obvious approach. In stage one, the periodogram has to be controlled at uncountably many frequencies at once. A union bound over $v$ is unavailable, and the leakage of each present frequency spreads over the whole spectrum, so both the noise and the signal must be bounded uniformly in $v$ with explicit constants. In stage two, the samples used to estimate an effective arm are selected with the estimated periods, which were computed from the same stage-one data, so the usual independence behind confidence bounds fails. Lemmas 9–11 are therefore stated for the weak oracle, whose sample selection does not depend on the data; on the event $\Lambda$ the definitions make the oracle's run and the policy's run the same function of the noise path.

Formalization scope

Epochs are natural numbers $t\ge1$; arms are Fin K (the paper's arm $k$ is index $k-1$), and $K\ge1$. The noise is a sequence $\epsilon:\mathbb N\to\Omega\to\mathbb R$ on a probability space, independent, with Mathlib's HasSubgaussianMGF at variance proxy $\sigma^2$, plus the printed tail clause. Both algorithms are defined in full, as deterministic functions of the noise path, with canonical tie-breaking: smallest maximizer, smallest closest candidate, smallest arm index. Algorithm 1's loop runs at most $n$ passes, and $\mathrm{LCM}(\hat v_i^{-1})$ is the least common multiple of the reduced denominators of the selected frequencies. Algorithm 2 takes its period vector as an argument; the two-stage policy and the oracle are its two instances. Expected regret is an outer expectation: an integrable majorant of the regret whose integral obeys the bound. Probability bounds bound the outer measure of the bad event. Every DFT uses the absolute stage-one epochs, as in Step 8 of Algorithm 1.

Disclosed deviations: Theorem 1's constant depends on $\sigma$ and $K$ (the printed proof does not give uniformity in $K$), and $\delta=8/T$ is fixed as in the proof of Lemma 11. Lemma 6 uses $0.2328$ for the printed $0.233$; Lemma 7 and Theorem 2 use $0.6934$ for the printed $0.694$. Lemma 7 covers the nonzero present frequencies. Lemma 11 states the standing Assumption 3 and the parameters (7) as hypotheses, as its proof uses them.

The algorithm is the object of the theorem. Replacing Algorithm 1 by "any estimator that is correct with high probability", or Algorithm 2 by "any policy with the oracle's regret", would make the goal a different and much weaker statement. So would placing the constant after the horizon or the instance, since the regret is at most $T$, or using a bare Lebesgue integral of a possibly non-measurable regret.

A complete development needs discrete Fourier analysis of finite sequences, Bernstein's inequality for complex polynomials on the unit disc (Lemma 14 of the paper, from Nagy and Totik), sinc-function estimates, anytime concentration for sub-Gaussian sums, and a bookkeeping layer for the nested index sets of Algorithm 2. The concentration and polynomial-inequality parts are reusable beyond this mission. Proofs of individual lemmas, helper lemmas such as the closed form (6) of the frequency components, and alternative arguments are welcome.

Selected references

• N. Chen, C. Wang, L. Wang, Learning and Optimization with Seasonal Patterns, arXiv:2005.08088v4, 2021; Operations Research, 2025. https://arxiv.org/abs/2005.08088
• P. Auer, Using Confidence Bounds for Exploitation-Exploration Trade-offs, Journal of Machine Learning Research 3, 397–422, 2002. https://www.jmlr.org/papers/v3/auer02a.html
• L. Li, Y. Lu, D. Zhou, Provably Optimal Algorithms for Generalized Linear Contextual Bandits, ICML 2017. https://arxiv.org/abs/1703.00048
• B. Nagy, V. Totik, Bernstein's Inequality for Algebraic Polynomials on Circular Arcs, Constructive Approximation 37, 223–232, 2013.

Motivation

Probabilistic (chance) constraints require that a decision keeps a random constraint satisfied with probability at least $p$. They are the standard way to model reliability requirements in engineering design, energy networks and the control of systems driven by uncertain data, and they are increasingly imposed on infinite-dimensional decisions such as the control of a partial differential equation. Geiersbach and Henrion (arXiv:2306.03965v2) derive first-order optimality conditions for such problems. This mission formalizes the general part of their paper (§2.1): a necessary and sufficient multiplier rule for a convex objective under a Gaussian probabilistic constraint on a Banach space.

The difficulty with chance constraints is that the probability function $\varphi(u)=\mathbb P(g(u,\xi)\le0)$ is given only implicitly and is in general neither smooth nor convex. Prékopa's theory of log-concave measures (Prékopa 1995) shows that $\log\varphi$ is concave under joint convexity, and Farshbaf-Shaker, Henrion and Hömberg (2018) carried this to infinite dimensions. For Gaussian data, van Ackooij and Henrion (2014) gave gradient formulae via the spherical-radial decomposition, and Hantoute, Henrion and Pérez-Aros (2019) proved local Lipschitz continuity of $\varphi$ and an estimate of its Clarke subdifferential under a moderate growth condition. The paper combines the two lines into an optimality system.

Setting

Let $U$ be a reflexive separable Banach space with dual $U^*$, let $F:U\to\mathbb R$ be convex and Fréchet differentiable with derivative $DF(u)\in U^*$, and let $\xi$ be an $m$-dimensional centered Gaussian random vector with nondegenerate covariance matrix $\Sigma$. A constraint function $g:U\times\mathbb R^m\to\mathbb R$ is locally Lipschitzian with $g(u,\cdot)$ convex for every $u$ (Assumption 2.1). For $p\in(0,1]$ the problem is

$$\min_{u\in U}F(u)\quad\text{s.t.}\quad\varphi(u):=\mathbb P(\omega\mid g(u,\xi(\omega))\le0)\ge p. \tag{9}$$

The Clarke directional derivative of a locally Lipschitzian $f:U\to\mathbb R$ is $f^\circ(u;h)=\limsup_{u'\to u,\,t\downarrow0}(f(u'+th)-f(u'))/t$, and its Clarke subdifferential is $\partial^Cf(u)=\{x^*\in U^*\mid\langle x^*,h\rangle\le f^\circ(u;h)\ \forall h\}$. At a point $u^*$ the paper assumes $g(u^*,0)<0$ (10) and either the moderate growth condition (12) or boundedness of $M=\{z\mid g(u^*,z)\le0\}$ (11).

The radial probability function is $e(u,v)=\mu_\eta\{r\ge0\mid g(u,r\Sigma^{1/2}v)\le0\}$ for $v$ on the unit sphere $\mathbb S^{m-1}$, where $\mu_\eta$ is the chi distribution with $m$ degrees of freedom; with $\mu_\zeta$ the uniform distribution on $\mathbb S^{m-1}$, $\varphi(u)=\int_{\mathbb S^{m-1}}e(u,v)\,d\mu_\zeta(v)$ (14). Integrals of set-valued maps are taken in the sense of Aumann.

Formalization targets

Goal: Proposition 2.6, statement (21)

If $g$ is jointly convex, there is a generalized Slater point $\hat u$ with $\varphi(\hat u)>p$, and $u^*$ is feasible with (10) and either (12) or $M$ bounded, then $u^*$ is a global minimizer of (9) if and only if, for some $\lambda\ge0$,

$$DF(u^*)\in\lambda\,\partial^C\varphi(u^*),\qquad\varphi(u^*)\ge p,\qquad\lambda(\varphi(u^*)-p)=0. \tag{21}$$

Milestones

1. Theorem 2.2, Lipschitz part: $e(\cdot,v)$ is uniformly locally Lipschitzian and $\varphi$ is locally Lipschitzian around $\bar u$.
2. Prékopa's statement (p. 7): joint quasi-convexity of $g$ and a log-concave density give a log-concave $\varphi$.
3. Lemma 2.5: $\tilde\varphi=-\log\varphi+\log p$ is convex and lower semicontinuous, and (9) is equivalent to $\min\{F\mid\tilde\varphi\le0\}$ (20).
4. The convex multiplier rule for (20): $-DF(u^*)\in\lambda\partial\tilde\varphi(u^*)$, $\lambda\tilde\varphi(u^*)=0$.
5. Convex and Clarke subdifferentials agree for locally Lipschitz convex functions.
6. The chain rule $\partial^C\tilde\varphi(u^*)=-\varphi(u^*)^{-1}\partial^C\varphi(u^*)$.

Companion targets

The spherical-radial representation (14); the inclusion (15) $\partial^C\varphi(\bar u)\subseteq\int_{\mathbb S^{m-1}}\partial^C_ue(\bar u,v)\,d\mu_\zeta(v)$ with equality under (16); the radius representation (17)–(19); Lemma 2.3 (Clarke regularity of $-e(\cdot,v)$); Theorem 2.4 (equality in (15) when $M$ is bounded); and Proposition 2.6's second statement (22), where $\partial^C\varphi(u^*)$ is replaced by the spherical integral.

Significance

Statement (21) is a complete first-order characterization of optimality for a nonconvex, nonsmooth constraint in infinite dimensions: the constraint function is neither differentiable nor convex in general, yet the multiplier rule is both necessary and sufficient. The companion formula (22), with the exact subdifferential of Theorem 2.4, makes the system computable from the explicit radial function $e$; the paper's later sections evaluate it for an elliptic PDE with a random source and derive explicit optimality conditions for a probabilistic state constraint.

The results of §2.1 are proved on paper, partly by citation. None of them has a machine-checked proof in Mathlib or on Prove2Me. Formalizing them requires a Clarke calculus on Banach spaces, Prékopa-type log-concavity for Gaussian measures, the spherical-radial decomposition of a Gaussian vector, and a convex multiplier rule with an extended-valued constraint. Each of these is reusable well beyond chance constraints.

Difficulty

The direct route fails because $\varphi$ is not convex: probability functions bounded between 0 and 1 are essentially never concave, and the Clarke subdifferential of a nonconvex function supports only necessary conditions. The paper's argument needs the convex reformulation through $-\log\varphi$, which is extended-valued where $\varphi=0$. It then needs three facts that do not come for free in a Banach space: local Lipschitz continuity of $\varphi$ (which fails for merely locally Lipschitz $g$ without (12) or a bounded $M$), the identification of the convex and Clarke subdifferentials, and a chain rule through the logarithm. The companion formula (22) needs Clarke's theorem on subdifferentials of integral functionals, whose equality case requires a regularity that the radial function only has after a sign change (Lemma 2.3).

Formalization scope

$\mathbb R^m$ is EuclideanSpace ℝ (Fin m) with $m\ge1$. The law of $\xi$ is Mathlib's multivariateGaussian 0 Sigma with Sigma.PosDef, and $\Sigma^{1/2}$ is the positive semidefinite square root used there. $\varphi(u)$ is the real number $\mathbb P(\{g(u,\xi)\le0\})$. Reflexivity is surjectivity of the canonical map into the bidual. The Clarke directional derivative is the limsup over $u'\to u$, $t\downarrow0$ with denominator $t$ (the page prints $h$, a misprint). Clarke regularity includes local Lipschitz continuity. $\tilde\varphi$ is valued in EReal with $-\log0=+\infty$, and its convexity is convexity of its epigraph. A solution is a global minimizer over the feasible set. "$DF(u^*)\in\lambda S$" is encoded as "$DF(u^*)=\lambda x^*$ for some $x^*\in S$", which is unambiguous at $\lambda=0$. The ball in (12) is open. The Aumann integral uses selections whose scalar components $v\mapsto\langle\hat x^*_v,h\rangle$ are integrable.

The convex multiplier rule (milestone 4) carries one disclosed addition: the extended-valued constraint is finite and continuous near $u^*$. The paper cites [39, Theorem 3.66] without spelling out a qualification condition, and in the application this property follows from milestone 1. Prékopa's statement is formalized with $\log0=-\infty$ (the page prints $+\infty$), as the multiplicative inequality $\varphi(u)^t\varphi(w)^{1-t}\le\varphi(tu+(1-t)w)$.

A formalization that makes any of these statements trivial is ruled out: the Clarke limsup is only used with a local Lipschitz hypothesis, so the junk value of an unbounded limsup never enters; the Aumann integral demands integrable selections; the multiplier carries its sign and complementarity condition; and optimality is global, not stationarity.

Contributions are welcome at every level: a Clarke calculus on normed spaces (subdifferential nonemptiness, sum and chain rules), log-concavity of Gaussian measures and Prékopa's theorem, polar coordinates for Gaussian vectors, and a Slater-type multiplier rule for extended-valued convex constraints.

Selected references

• C. Geiersbach and R. Henrion, Optimality conditions in control problems with random state constraints in probabilistic or almost-sure form, arXiv:2306.03965v2, 2024 (Math. Oper. Res., 2025). https://arxiv.org/abs/2306.03965
• A. Hantoute, R. Henrion and P. Pérez-Aros, Subdifferential characterization of probability functions under Gaussian distribution, Math. Program. 174, 167–194, 2019. https://doi.org/10.1007/s10107-018-1237-9
• W. van Ackooij and R. Henrion, Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussian-like distributions, SIAM J. Optim. 24(4), 1864–1889, 2014. https://doi.org/10.1137/130922689
• M. H. Farshbaf-Shaker, R. Henrion and D. Hömberg, Properties of chance constraints in infinite dimensions with an application to PDE constrained optimization, Set-Valued Var. Anal. 26, 821–841, 2018. https://doi.org/10.1007/s11228-017-0452-5
• A. Prékopa, Stochastic Programming, Kluwer, 1995. https://doi.org/10.1007/978-94-017-3087-7
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, 1983. https://doi.org/10.1137/1.9781611971309
• J. Peypouquet, Convex Optimization in Normed Spaces: Theory, Methods and Examples, SpringerBriefs in Optimization, Springer, 2015.

Motivation

Optimal control of partial differential equations under uncertainty arises whenever a physical system is steered by a control while its sources, coefficients or loads are only known in distribution. A typical requirement is a state constraint: the state (a temperature, a concentration, a displacement) must stay below a threshold everywhere in the domain. When the source term is random, demanding this for every realization is often impossible or too conservative, and the natural relaxation is a probabilistic (chance) constraint: the state must stay below the threshold uniformly on the domain with probability at least $p$.

Probabilistic constraints are classical in operations research (Prékopa, Stochastic Programming, 1995), but in a PDE setting the constraint is an infinite family of random inequalities, one per point of the domain, and the probability function is not differentiable in general. C. Geiersbach and R. Henrion (arXiv:2306.03965v2, Math. Oper. Res. 2025) derive optimality conditions for such problems in terms of the Clarke subdifferential, building on the spherical-radial decomposition of Gaussian vectors developed by Hantoute, Henrion and Pérez-Aros (Math. Program. 2019) and van Ackooij and Henrion. This mission formalizes the PDE half of §2 of that paper: explicit formulas for the subdifferential of the probability function and the resulting optimality conditions.

Setting

Let $D\subseteq\mathbb R^d$, $d\in\{1,2,3\}$, be a nonempty, open, bounded set of class S: there are $\gamma\in(0,1)$ and $r_0>0$ with $\operatorname{meas}(B_r(x)\setminus D)\ge\gamma\operatorname{meas}(B_r(x))$ for all $x\in\partial D$ and $0<r<r_0$. For $F\in L^2(D)$ the Dirichlet problem $-\Delta y=F$ in $D$, $y=0$ on $\partial D$ has a unique weak solution in $H^1_0(D)\cap C(\bar D)$; the solution operator $A^{-1}:L^2(D)\to C(\bar D)$ is linear and continuous.

Fix $f_0,\varphi_1,\dots,\varphi_m\in L^2(D)$ and a threshold $\alpha\in\mathbb R$. The control-to-state operator is

$$S(u,z)=A^{-1}\Big(u+f_0+\sum_{i=1}^m z_i\varphi_i\Big)=P(u,z)+y_0,\qquad P(u,z)=A^{-1}u+\sum_i z_iA^{-1}\varphi_i,\ y_0=A^{-1}f_0.$$

The random vector $\xi$ on a probability space $(\Omega,\mathcal F,\mathbb P)$ is Gaussian $\mathcal N(0,\Sigma)$ with $\Sigma$ positive definite, and $y(\cdot,\omega)=S(u,\xi(\omega))$ is the random state. With $g(u,z)=\max_{x\in\bar D}[S(u,z)](x)-\alpha$ and the probability function $\varphi(u)=\mathbb P(\omega\mid g(u,\xi(\omega))\le0)$, the problem is

$$\min_{u\in L^2(D)}F(u)\quad\text{s.t.}\quad \varphi(u)=\mathbb P\big(y(x,\omega)\le\alpha\ \ \forall x\in D\big)\ge p,\tag{24}$$

with $F$ convex and Fréchet differentiable and $p\in(0,1]$.

The analysis uses the spherical-radial decomposition. Let $\mu_\zeta$ be the uniform distribution on the unit sphere $\mathbb S^{m-1}$, $\mu_\eta$ the chi distribution with $m$ degrees of freedom and density $\chi(t)=Kt^{m-1}e^{-t^2/2}$, $\Sigma^{1/2}$ the square root of $\Sigma$. The radial probability function is $e(u,v)=\mu_\eta\{r\ge0\mid g(u,r\Sigma^{1/2}v)\le0\}$ and the radius function $\rho(u,v)\in[0,\infty]$ is the supremum of that set. For $x\in\bar D$, $u_x\in L^2(D)$ is defined by $(u_x,h)_{L^2(D)}=(A^{-1}h)(x)$, and

$$M^*(\bar u,v)=\{x\in\bar D\mid[P(0,\Sigma^{1/2}v)](x)>0,\ \bar y(x)+\rho(\bar u,v)[P(0,\Sigma^{1/2}v)](x)=\alpha\},\qquad \bar y=S(\bar u,0).$$

$\partial^C$ is the Clarke subdifferential, with $L^2(D)^*$ identified with $L^2(D)$, and integrals of set-valued maps are Aumann integrals.

Formalization targets

Goal: Theorem 2.16

Assume a control $\hat u$ with $\varphi(\hat u)>p$ (41), and let $u^*$ be feasible with expected state $\tilde y=S(u^*,0)$ satisfying $\tilde y(x)<\alpha$ on $\bar D$ (42). If $u^*$ solves (24), there is $\lambda\ge0$ with

$$DF(u^*)\in-\lambda\int_{\{v\mid\rho(u^*,v)<\infty\}}\chi(\rho(u^*,v))\cdot\operatorname{clco}\Big\{\frac{u_x}{[P(0,\Sigma^{1/2}v)](x)}\ \Big|\ x\in M^*(u^*,v)\Big\}\,d\mu_\zeta(v),$$ $$\lambda\big(\mathbb P(y^*(x,\omega)\le\alpha\ \forall x\in D)-p\big)=0,$$

and these conditions are sufficient when $M^*(u^*,v)$ is a singleton for $\mu_\zeta$-a.e. $v$ with $\rho(u^*,v)<\infty$, or the set $\{z\mid g(u^*,z)\le0\}$ is bounded. The goal has no hard-coded constants; it asserts the form of the optimality system.

Milestones

In attack order: well-posedness of the state equation with the bound $\|\hat y\|_{C(\bar D)}\le C(\|z\|_2+\|u+f_0\|_{L^2(D)})$ (Lemma 1.2); measurability of the random state (Lemma 1.4); convexity and global Lipschitz continuity of $g$ (Lemma 2.9); moderate growth (Lemma 2.10); the Lipschitz and subdifferential conclusions of the general theory under (26) (Corollary 2.11); the optimality conditions (21), (22) for (24) (Corollary 2.12); $\partial^C_ue(\bar u,v)=\{0\}$ when $\rho(\bar u,v)=\infty$ (Lemma 2.13); the explicit formula $\partial^C_ue(\bar u,v)=-\chi(\rho(\bar u,v))\cdot\operatorname{clco}\{u_x/[P(0,\Sigma^{1/2}v)](x)\mid x\in M^*(\bar u,v)\}$ when $\rho(\bar u,v)<\infty$ (Lemma 2.14); and the upper estimate (40) of $\partial^C\varphi(\bar u)$ with its equality cases (Theorem 2.15).

Significance

The result turns an abstract statement — $DF(u^*)$ lies in a multiple of the Clarke subdifferential of a probability function that has no closed form — into a condition written in the problem data: Green-type functions $u_x$, the expected state, and the radius function along each direction of the sphere. This is what makes the conditions usable for checking candidate controls and for designing algorithms, and Example 2.18 of the paper computes a solution this way. The sufficiency half identifies when the nonconvex chance-constrained problem behaves like a convex one.

The paper's proofs are on paper only; to our knowledge none of these results, nor the spherical-radial decomposition, Clarke calculus in Hilbert spaces, or Aumann integrals in $L^2$, has a machine-checked proof. Formalizing Theorem 2.16 requires connecting elliptic regularity, Gaussian measure theory and nonsmooth analysis, and the milestones are reusable on their own: Lemma 1.2 is a regularity statement for the Poisson equation, and Lemmas 2.13 and 2.14 are max-function subdifferential formulas in an infinite-dimensional space.

Difficulty

The obvious route — differentiate $\varphi$ and use a Lagrange multiplier — fails because $\varphi$ is generally not differentiable even for smooth data (Example 2.17 of the paper), and because $g$ is only Lipschitz, being a maximum over $\bar D$. The general subdifferential estimate (Theorem 2.2, from Hantoute–Henrion–Pérez-Aros) needs a growth condition and only gives an inclusion, since the radial probability function is not Clarke regular; the radius function $\rho(\cdot,v)$ is an infimum over the open set where $P(0,\Sigma^{1/2}v)>0$, which is neither compact nor independent of $u$, so standard max-function formulas do not apply directly. Directions with $\rho=\infty$ require a separate limiting argument through the Mordukhovich subdifferential (Lemma 2.13). The regularity input, continuity of solutions on $\bar D$ for class-S domains (Kinderlehrer–Stampacchia 1980), is not in Mathlib.

Formalization scope

All declarations are in the namespace ChanceStateCtrl.PDE. $\mathbb R^d$ and $\mathbb R^m$ are EuclideanSpace ℝ (Fin d) and EuclideanSpace ℝ (Fin m), indices are Fin m, $L^2(D)$ is Lp ℝ 2 (volume.restrict D). Weak derivatives are the published HunterPDE.Shared.HasWeakDeriv; a weak solution is a function continuous on $\bar D$, zero on $\partial D$, in $L^2(D)$ with weak first derivatives in $L^2(D)$, satisfying the weak form against $C_c^\infty(D)$.

Standing hypotheses on every theorem: $D$ open, bounded and nonempty, $d\in\{1,2,3\}$, class S, and CompactSpace (closure D) (a consequence of boundedness). Nonemptiness of $D$ is the one added hypothesis: Lean's supremum over an empty domain is $0$, not $-\infty$, and Lemma 2.9's identification of events fails for $D=\emptyset$. $A^{-1}$ is a continuous linear map to $C(\bar D)$ together with the hypothesis that it returns the weak solution; it is never an arbitrary operator. $m\ge1$, as the uniform distribution on $\mathbb S^{m-1}$ requires. $\Sigma^{1/2}$ is the positive semidefinite root; $\|\Sigma^{1/2}\|$ is its operator norm. The Clarke directional derivative is $\limsup_{u'\to u,t\downarrow0}(f(u'+th)-f(u'))/t$ (the page's denominator $h$ is a misprint). $\rho$ takes values in $[0,\infty]$. clco is the norm closure of the convex hull. The integral in (40)/(43) is over the restriction of $\mu_\zeta$ to $\{\rho<\infty\}$, not a renormalization. A solution is a global minimizer over $\{u\mid\varphi(u)\ge p\}$; the page's "solution to (9)" refers to this instance of the abstract problem (9).

Trivializing formalizations are ruled out: every Clarke subdifferential is taken at a point where the function is locally Lipschitz (proved in Corollary 2.11 and Theorem 2.15), the state operator is tied to the PDE, the multiplier carries its sign and complementarity, Aumann selections must be scalarly integrable, and optimality is global.

Needed infrastructure: elliptic regularity in $C(\bar D)$ for class-S domains; the spherical-radial decomposition of $\mathcal N(0,\Sigma)$; Clarke calculus (chain rule, max functions) on Hilbert spaces; Clarke's theorem on subdifferentials of integral functionals. Contributions of any of these as stand-alone theorems are welcome.

Selected references

• C. Geiersbach, R. Henrion, Optimality conditions in control problems with random state constraints in probabilistic or almost-sure form, arXiv:2306.03965v2, 2024; Math. Oper. Res., 2025. https://arxiv.org/abs/2306.03965
• A. Hantoute, R. Henrion, P. Pérez-Aros, Subdifferential characterization of probability functions under Gaussian distribution, Math. Program. 174 (2019), 167–194. https://doi.org/10.1007/s10107-018-1237-9
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, 1983; SIAM reprint 1990. https://doi.org/10.1137/1.9781611971309
• D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, Academic Press, 1980; SIAM reprint 2000. https://doi.org/10.1137/1.9780898719451
• A. Prékopa, Stochastic Programming, Kluwer, 1995. https://doi.org/10.1007/978-94-017-3087-7