Operations Research

Motivation

Airlines sell seats on a single flight in several fare classes. Discount classes book first, so the carrier must decide how many seats to protect for later, higher-paying passengers. The classical answer, from Littlewood's rule for two classes to the nested optimality conditions of Brumelle and McGill (1993), computes optimal protection levels from a forecast of the demand distribution of each class. Those forecasts must be built from censored booking data, which is where most of the practical difficulty lies.

van Ryzin and McGill (2000) observed that the optimality conditions themselves can drive an adaptive rule. After each flight one only records whether a fill event occurred (whether demand reached each protection level) and moves the levels by a stochastic-approximation step. No forecast and no optimization is needed. Their Theorem 1 asserts that this rule converges almost surely to the optimal protection levels, with an explicit mean-square rate. It is one of the standard references for model-free revenue management. The convergence proof combines Robbins–Monro theory, the Robbins–Siegmund supermartingale lemma and an induction over fare classes.

Setting

There are $k+1$ fare classes with fares $f_1 > f_2 > \cdots > f_{k+1} > 0$ and discount ratios $r_{i+1} = f_{i+1}/f_1$. On a flight the class demands $X_1, \dots, X_{k+1}$ are independent, nonnegative and continuously distributed. Successive flights $X^1, X^2, \dots$ are independent with the same law.

A protection vector is $\theta = (\theta_1, \dots, \theta_k)$. The fill events are

$$A_i(\theta, X) = \{X_1 > \theta_1,\ X_1 + X_2 > \theta_2,\ \dots,\ X_1 + \cdots + X_i > \theta_i\}, \qquad i = 1, \dots, k.$$

A vector $\theta^*$ is characterised by the optimality condition (3), $P(A_i(\theta^*, X)) = r_{i+1}$ for $i = 1, \dots, k$. The adjustments are $H_i(\theta, X) = r_{i+1} - \mathbf 1(A_i(\theta, X))$, and their means are $h_i(\theta) = r_{i+1} - P(A_i(\theta, X))$. From an arbitrary $\theta^1$ the algorithm updates

$$\theta^{n+1} = \theta^n - \gamma_n H(\theta^n, X^n), \qquad \gamma_n = \frac{A}{n + B},\ A > 0,\ B \ge 0.$$

The interim protection levels $p_i(\theta) = \max\{\theta_j : 1 \le j \le i\}$ are what the airline actually uses to control bookings. The updates themselves monitor $A_i(\theta^n, X^n)$.

Formalization targets

Goal: Theorem 1 (p. 765)

Assume:

• A1: every $X_i$ has bounded support.
• A2 (windowed): for each $R > 0$ there is $\delta_R > 0$ with $(\theta_i - \theta^*_i)\,h_i(\theta_i, \theta^*_{i-1}, \dots, \theta^*_1) \ge \delta_R|\theta_i - \theta^*_i|^2$ whenever $|\theta_i - \theta^*_i| \le R$.
• A3: the distribution functions of the partial sums $X_1 + \cdots + X_i$ are Lipschitz.

Then for $i = 1, \dots, k$,

$$\theta^n_i \to \theta^*_i \ \text{a.s.}, \qquad p_i(\theta^*) = \theta^*_i,$$

and there are $\beta > 0$ and $C$ with

$$E|\theta^n_i - \theta^*_i|^2 \le C\,\gamma_n^{\beta/2^{i-1}}.$$

The exponent is left existential, as in the paper, so any improvement of the rate constant still proves the goal.

Milestones

The milestones follow the paper's own proof, in attack order:

• the base case $i = 1$ (p. 765);
• the Lipschitz bound on $h_{i+1}$ (p. 766);
• Lemma 3, boundedness of the iterates (p. 764);
• the almost-supermartingale inequality (12) (p. 766);
• Lemma 1 (p. 764);
• the a.s. summability of $\gamma_n|\theta^n_i - \theta^*_i|$ (p. 766);
• Lemma 2, Robbins–Siegmund (p. 764);
• the induction step for (9) (p. 766);
• Lemma 4, the rate recursion (p. 764).

Significance

Theorem 1 says that the nested optimality conditions of Brumelle and McGill can be reached by observing binary fill events alone. It also says that the limit is automatically ordered, $\theta^*_1 \le \cdots \le \theta^*_k$, so the interim levels coincide with the limit. The algorithm is distribution-free and needs neither demand forecasts nor uncensoring, which is why it is used in practice and as a baseline for later data-driven and learning-based revenue management.

The theorem is proved in the paper. As far as is known it has no machine-checked proof. Formalizing it requires a vector Robbins–Monro argument with coupled coordinates, the Robbins–Siegmund lemma, and an explicit polynomial rate. Two of the milestones are general probability facts reusable far beyond this mission: Lemma 1, and Robbins–Siegmund, which no platform item states yet. The formalization also corrects the record on two printed slips, described under Formalization scope.

Difficulty

Coordinate $i+1$ is not a Robbins–Monro process by itself. Its fill event $A_{i+1}(\theta, X)$ depends on all of $\theta_1, \dots, \theta_{i+1}$, so the drift of $\theta^n_{i+1}$ is perturbed by the errors of the lower coordinates. Applying a scalar convergence theorem coordinate by coordinate fails at this point. The perturbation has to be bounded via the Lipschitz assumption A3 and shown to be summable along the path, which uses the mean-square rate already established for class $i$. The rate and the almost-sure convergence must therefore be carried through the induction together. The exponent halves at each class. A second difficulty is that A2 is only usable on a bounded set, so boundedness of the iterates (Lemma 3) has to come first and must be uniform.

Formalization scope

Classes are indexed from $1$ in Lean ($\mathbb N$-indexed vectors whose index $0$ and indices above $k+1$ are never read). A flight's demand vector has the product law $\bigotimes_i \nu_i$, and the sample path has the countable product of copies of it. Lean's iterate … n is the paper's $\theta^{n+1}$, so the rate pairs it with $\gamma_{n+1}$. Continuity of a demand law is "every singleton is null". Expectations are Bochner integrals, and integrability is part of each rate conclusion. The fill event reuses the published NestedSeatAlloc.ProbCond.nestEvent.

Pinned and corrected statements:

• A2 as printed is unsatisfiable. It asks for one $\delta$ valid for all $\theta_i$. Since $|h_i| \le 1$, the inequality fails once $|\theta_i - \theta^*_i| > 1/\delta$, so the printed Theorem 1 is vacuously true. The mission uses the windowed A2 above. That is what the proof applies to bounded iterates, and what the paper's sufficient condition (a density bounded below near $\theta^*$) yields. A formalization with the printed A2, or with any other hypothesis no instance satisfies, is ruled out.
• The base case. The page claims the rate "for all $0 < \beta < 1$". That is false in general. The milestone states "for some $\beta > 0$", which is all Theorem 1 uses.
• Integrability. Lemma 2 assumes each $Z_n$ integrable, so that the conditional expectation is genuine. Lemma 3's "bounded (a.s.)" is read as a bound uniform in $n$.
• Standing assumptions. Nonnegative demands and positive fares are added where the argument uses them.

Not posed: §4's modifications (projection onto the capacity, randomized integer levels, booking lead times) and the simulations of §5.

Contributions welcome: proofs of the general lemmas (Lemma 1, Robbins–Siegmund) stated for arbitrary probability spaces, measurability facts for the recursion, and the elementary Lemma 4.

Selected references

• G. van Ryzin and J. McGill, Revenue Management Without Forecasting or Optimization: An Adaptive Algorithm for Determining Airline Seat Protection Levels, Management Science 46(6), 760–775, 2000. https://doi.org/10.1287/mnsc.46.6.760.11936
• S. L. Brumelle and J. I. McGill, Airline Seat Allocation with Multiple Nested Fare Classes, Operations Research 41(1), 127–137, 1993. https://doi.org/10.1287/opre.41.1.127
• H. Robbins and D. Siegmund, A convergence theorem for non negative almost supermartingales and some applications, in Optimizing Methods in Statistics, Academic Press, 233–257, 1971. https://doi.org/10.1016/B978-0-12-604550-5.50015-8
• A. Benveniste, M. Métivier and P. Priouret, Adaptive Algorithms and Stochastic Approximations, Springer, 1990. https://doi.org/10.1007/978-3-642-75894-2
• E. Lukacs, Stochastic Convergence, 2nd ed., Academic Press, 1975.

Motivation

A simulation metamodel is a cheap surrogate for an expensive stochastic simulation: after running the simulation at a few design settings, one predicts the mean response at settings never simulated. Kriging, the interpolation method of geostatistics and of the design and analysis of computer experiments (DACE; Sacks, Welch, Mitchell and Wynn 1989; Santner, Williams and Notz 2003), treats the unknown response surface as a realization of a random field and predicts by the best linear predictor. Kriging was developed for deterministic computer codes, where an observation is the exact response. A stochastic simulation returns the response plus sampling noise, and the noise may be correlated across settings when the experimenter uses common random numbers (CRN).

Ankenman, Nelson and Staum, Stochastic Kriging for Simulation Metamodeling (Proc. 2008 Winter Simulation Conference, pp. 362–370), extend kriging to this setting. The mission formalizes §2 of that paper: the form of the optimal linear predictor, called stochastic kriging, and the formula for its mean squared error, which shows how much the simulation noise costs. The source used throughout is the WSC 2008 proceedings paper, not the later Operations Research 58(2) 2010 article, whose numbering differs.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and a set of design settings $\mathbf x$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation outputs (display (3), with constant trend)

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x).$$

The constant $\beta_0$ is the overall mean. $\mathsf M$ is a mean-zero random field, the extrinsic uncertainty imposed by the modeller. $\varepsilon_j(\mathbf x)$ is mean-zero sampling noise, the intrinsic uncertainty of the simulation; its variance may depend on $\mathbf x$, and noises at different settings may be correlated (CRN). All of these have finite second moments, and the field is uncorrelated with the noise.

An experiment design runs $n_i\ge1$ replications at each of $k$ design points $\mathbf x_1,\dots,\mathbf x_k$. The data are the sample means (4)

$$\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i),\qquad\bar{\mathcal Y}=(\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))^\top,$$

and the target at a point $\mathbf x_0$ is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

Write $\Sigma_{\mathsf M}(\mathbf x,\mathbf x')=\mathrm{Cov}[\mathsf M(\mathbf x),\mathsf M(\mathbf x')]$. Then $\Sigma_{\mathsf M}$ is the $k\times k$ matrix of these covariances at the design points, $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is the vector $(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$, and $\Sigma_\varepsilon$ is the $k\times k$ covariance matrix of the averaged noises $\sum_{j=1}^{n_i}\varepsilon_j(\mathbf x_i)/n_i$. A linear predictor (5) is $\lambda_0+\lambda^\top\bar{\mathcal Y}$ with arbitrary weights $(\lambda_0,\lambda)\in\mathbb R\times\mathbb R^k$; its mean squared error is $\mathrm E[(\lambda_0+\lambda^\top\bar{\mathcal Y}-\mathsf Y(\mathbf x_0))^2]$. The optimal MSE $\mathrm{MSE}^\star$ is the infimum of this quantity over all weights.

Formalization targets

Goal: display (7)

There is a rule $\Xi$, assigning a positive definite matrix $\Xi(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ to every pair of positive definite matrices, such that in every model as above with $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$, and at every $\mathbf x_0$,

$$\begin{aligned}\mathrm{MSE}^\star&=\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\\&=\Big[\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\Sigma_{\mathsf M}^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\Big]+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\,\Xi\,\Sigma_{\mathsf M}(\mathbf x_0,\cdot).\end{aligned}$$

The bracket is the usual kriging MSE. The goal also states that $\mathrm{MSE}^\star$ strictly exceeds it whenever $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\neq0$. The matrix $\Xi$ is existential and depends only on the two covariance matrices, not on $\mathbf x_0$.

Milestone: display (6)

If $\Sigma_{\mathsf M}+\Sigma_\varepsilon\succ0$, the stochastic kriging predictor

$$\widehat{\mathsf Y}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}(\bar{\mathcal Y}-\beta_0\mathbf 1_k)$$

is of the form (5), attains the smallest MSE among all predictors of that form, and is the only one that does.

Further items (not milestones)

• (9)–(10): the two-point example, with closed forms for the predictor and its MSE in terms of $\tau^2$, $r_{12}$, $r_0$, $\mathsf V$, $\rho$ and $n$.
• (14): the equicorrelated example, $\mathrm{MSE}^\star=\tau^2\big(1-k r_0^2/(1+(k-1)r+\gamma/n)\big)$.
• §3, p. 365: under the Gaussian Assumption 1, $\widehat{\mathsf Y}(\mathbf x_0)=\mathrm E[\mathsf Y(\mathbf x_0)\mid\bar{\mathcal Y}]$ almost surely.

Significance

Display (7) quantifies the price of simulation noise. With exact observations the best achievable error is the kriging MSE; with noisy observations it is larger by a positive definite form in the cross-covariance vector. Every point $\mathbf x_0$ that is correlated with the design therefore loses accuracy. The two examples make this concrete. In (10) the MSE increases with the intrinsic correlation $\rho$, which is why the paper concludes that common random numbers, a standard variance-reduction device for comparing systems, do not help prediction. Display (14) is the baseline against which the paper measures the cost of estimating the noise variance. The conditional-expectation statement shows that under Gaussian assumptions stochastic kriging is optimal among all predictors, not only linear ones.

The results are classical in substance: they are the best-linear-prediction formulas of second-order random-field theory, applied to data whose covariance is $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The paper states them without proof ("we can show"). No machine-checked version exists. Formalizing them yields a reusable, measure-theoretic best-linear-predictor theorem for random vectors defined from a model rather than postulated, together with the Loewner-order fact behind the positive definiteness of $\Xi$.

Difficulty

The algebra is short; the difficulty is in deriving the second-moment structure from the model, not assuming it. The covariance of $\bar{\mathcal Y}$ must be computed from the sample-mean definition, using bilinearity of covariance over finite sums of square-integrable variables, and shown to equal $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The cross-covariance with $\mathsf Y(\mathbf x_0)$ must be shown to equal $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$. The MSE of an arbitrary affine predictor then has to be written as a quadratic in $(\lambda_0,\lambda)$, and its infimum identified.

A first idea is to take "a random vector with mean $\beta_0\mathbf 1$ and covariance $\Sigma$" as the hypothesis. That proves a different, more abstract statement and leaves the model's content unproved.

For the second line of (7), the positive definiteness of $\Xi$ does not follow from a scalar argument. It is a matrix inequality, $(A+B)^{-1}\prec A^{-1}$ for $A,B\succ0$, and it fails without $\Sigma_\varepsilon\succ0$: when $\Sigma_\varepsilon=0$ the two lines coincide.

Formalization scope

• The model. Settings form an arbitrary type X; design points are x : Fin k → X, so the paper's index $i$ is i.val + 1. Replication $j\ge1$ is index j - 1 of a sequence ε : ℕ → X → Ω → ℝ.
• Moments and covariances. Field and noise are in $L^2$ and have mean zero. Covariances are Mathlib's covariance. The optimal MSE is an infimum over all of $\mathbb R\times\mathbb R^k$, with no unbiasedness constraint and no sign constraint.
• Hypotheses not printed in §2. Three are added and disclosed:
  1. uncorrelatedness of $\mathsf M$ and $\varepsilon$, the second-order content of Assumption 1's "independent of $\mathsf M$", without which the covariance of $\bar{\mathcal Y}$ is not $\Sigma_{\mathsf M}+\Sigma_\varepsilon$;
  2. square integrability;
  3. $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$ in the goal. Lean's matrix inverse returns $0$ for a singular matrix, and $\Xi\succ0$ is false when $\Sigma_\varepsilon$ is singular.
• Generality kept. No Gaussianity, no independence across replications and no "no CRN" restriction is imposed on (6) and (7).
• Examples. $|r_{12}|<1$ in (10) and $r\le1$ in (14) are added, since both are correlations.
• Non-trivializing. A positive semidefinite $\Xi$, or a $\Xi$ chosen after $\mathbf x_0$, would make the goal trivial or strictly weaker; the statement fixes $\Xi$ as a function of $(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ before the model and the point.

Infrastructure needed and reusable beyond this mission:

• covariance of averages of $L^2$ random variables;
• the MSE of an affine predictor as a quadratic form;
• minimization of a positive definite quadratic;
• the Loewner antitonicity of the inverse;
• for the conditional-expectation item, Gaussian conditioning via uncorrelated-hence-independent residuals.

Contributions to any of these are welcome.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2), 371–382, 2010 (later journal version, not used here). https://doi.org/10.1287/opre.1090.0754
• J. Sacks, W. J. Welch, T. J. Mitchell, H. P. Wynn, Design and Analysis of Computer Experiments, Statistical Science 4(4), 409–423, 1989. https://doi.org/10.1214/ss/1177012413
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8

Motivation

Stochastic simulation models of queues, supply chains and manufacturing systems are often too slow to run at every input setting of interest. A metamodel fitted to a moderate number of simulation runs then stands in for the simulation, for instance inside an optimization or a sensitivity analysis. Kriging, developed in geostatistics and in the design and analysis of deterministic computer experiments (DACE), predicts an unknown response surface by treating it as a realization of a Gaussian random field. Deterministic kriging interpolates its data exactly, which is wrong for stochastic simulation, where each run returns the response plus sampling noise.

Ankenman, Nelson and Staum's stochastic kriging separates the two sources of uncertainty: the extrinsic uncertainty of the unknown surface and the intrinsic uncertainty of the simulation output. Its optimal predictor requires the intrinsic noise variances, which are never known in practice and are replaced by sample variances computed from the same replications. This mission formalizes the paper's first key result: this plug-in step introduces no prediction bias.

The source is the Winter Simulation Conference 2008 proceedings version of the paper (pp. 362–370). The later Operations Research 58(2) (2010) article numbers its results differently and is not the reference here.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and design variables $\mathbf x\in\mathbb R^d$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation returns

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x),$$

where $\beta_0\in\mathbb R$ is a constant trend, $\mathsf M$ is a mean-zero random field on $\mathbb R^d$ (the unknown surface) and $\varepsilon_j(\mathbf x)$ is the noise of replication $j$. The quantity to predict at any point $\mathbf x_0$, simulated or not, is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

An experiment design is a list of pairs $(\mathbf x_i,n_i)$, $i=1,\dots,k$, with distinct design points $\mathbf x_i$ and $n_i$ replications at $\mathbf x_i$. The data are summarized by the sample means $\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i)$, collected in the vector $\bar{\mathcal Y}\in\mathbb R^k$, and the sample variances

$$\mathcal S^2(\mathbf x_i)=\frac1{n_i-1}\sum_{j=1}^{n_i}\big(\mathcal Y_j(\mathbf x_i)-\bar{\mathcal Y}(\mathbf x_i)\big)^2 .$$

The extrinsic covariances are the $k\times k$ matrix $\Sigma_{\mathsf M}=(\mathrm{Cov}[\mathsf M(\mathbf x_h),\mathsf M(\mathbf x_i)])_{h,i}$ and the vector $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)=(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$.

Assumption 1. $\mathsf M$ is a stationary Gaussian random field: all its finite-dimensional laws are multivariate normal with mean $0$, the covariance is $\tau^2R(\mathbf x-\mathbf x')$ with $\tau^2>0$ and $R(\mathbf 0)=1$, and the covariance matrix at finitely many distinct points is positive definite. At each design point the noises $\varepsilon_1(\mathbf x_i),\varepsilon_2(\mathbf x_i),\dots$ are i.i.d. $N(0,\mathsf V(\mathbf x_i))$ with $\mathsf V(\mathbf x_i)>0$; noises at different design points are independent (no common random numbers); and the whole noise family is independent of $\mathsf M$.

The intrinsic variance at a design point is estimated by $\widehat{\mathsf V}(\mathbf x_i)=\mathcal S^2(\mathbf x_i)$, giving $\widehat\Sigma_\varepsilon=\mathrm{Diag}\{\widehat{\mathsf V}(\mathbf x_1)/n_1,\dots,\widehat{\mathsf V}(\mathbf x_k)/n_k\}$ and the plug-in stochastic kriging predictor

$$\widehat{\widehat{\mathsf Y}}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\big[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon\big]^{-1}\big(\bar{\mathcal Y}-\beta_0\mathbf 1_k\big).\tag{13}$$

Formalization targets

Goal: Theorem 1 (p. 366)

Under Assumption 1, with $n_i\ge2$ replications at every design point, the prediction error is integrable and

$$\mathrm E\Big[\widehat{\widehat{\mathsf Y}}(\mathbf x_0)-\mathsf Y(\mathbf x_0)\Big]=0 .$$

Milestones (p. 365)

1. Under Assumption 1, $(\mathsf Y(\mathbf x_0),\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))$ is multivariate normal.
2. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ has a scaled chi-squared distribution: $(n_i-1)\mathcal S^2(\mathbf x_i)/\mathsf V(\mathbf x_i)\sim\chi^2_{n_i-1}$.
3. Under Assumption 1, $\mathcal S^2(\mathbf x_i)$ is strongly consistent for $\mathsf V(\mathbf x_i)$: the sample variance of the first $m$ replications converges to $\mathsf V(\mathbf x_i)$ almost surely as $m\to\infty$.

Significance

Theorem 1 says that the practitioner's shortcut, estimating the noise variances from the replications and plugging them into the optimal predictor, keeps the predictor unbiased; the cost of not knowing the intrinsic variance is paid entirely in mean squared error. This is what makes the paper's subsequent comparison of the plug-in MSE with the known-variance MSE the right measure of the penalty for estimation.

None of these statements has a machine-checked proof, and Mathlib has neither the chi-squared law of the normal sample variance nor the independence of the normal sample mean and sample variance. A complete development produces these classical facts of normal sampling theory for the first time in Lean, together with a reusable model of noisy observations of a Gaussian random field.

Difficulty

The predictor is a nonlinear function of the data: its weight vector $[\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is random and depends on the same replications as $\bar{\mathcal Y}$. Linearity of expectation therefore does not apply directly: the obvious argument, which treats the weights as constants, fails, and the dependence between the estimated variances, the sample means and the field has to be controlled exactly. The joint independence structure of Assumption 1 — across replications, across design points, and between noise and field — is load-bearing; a weaker pairwise or pointwise independence does not suffice. Integrability of the error must also be established, since the weights depend on random variances.

Formalization scope

Design points are indexed by Fin k (the paper's $i$ is i.val + 1) and lie in EuclideanSpace ℝ (Fin d); replication $j\ge1$ of the paper is index j - 1 of a sequence, so sums run over Finset.range. The field is M : EuclideanSpace ℝ (Fin d) → Ω → ℝ with Mathlib's IsGaussianProcess; covariances are ProbabilityTheory.covariance; stationarity is $\mathrm{Cov}[\mathsf M(\mathbf y),\mathsf M(\mathbf y')]=\tau^2R(\mathbf y-\mathbf y')$ with $R(\mathbf 0)=1$. The noise law is HasLaw … (gaussianReal 0 (V (x i))); independence is one iIndepFun over the family indexed by Fin k × ℕ plus one IndepFun between that whole family and the whole field. The variance function is ℝ≥0-valued with $\mathsf V(\mathbf x_i)>0$ assumed; replication counts satisfy $n_i\ge2$ wherever $\mathcal S^2$ appears and $n_i\ge1$ where only means appear. The chi-squared law $\chi^2_\nu$ is gammaMeasure (ν/2) (1/2).

The covariance quantities $\Sigma_{\mathsf M}$ and $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ in (13) are the true ones; only $\Sigma_\varepsilon$ is estimated, and $\widehat\Sigma_\varepsilon$ is built from the same replications as $\bar{\mathcal Y}$. A formalization in which $\widehat\Sigma_\varepsilon$ is deterministic or independent of the data, in which the matrix $\Sigma_{\mathsf M}+\widehat\Sigma_\varepsilon$ may be singular (Lean's Matrix.inv returns $0$ there, reducing the claim to $\mathrm E[\mathsf M(\mathbf x_0)]=0$), or in which the expectation is asserted without integrability, would trivialize the goal and is ruled out: positive definiteness of $\Sigma_{\mathsf M}$ at distinct design points is part of Assumption 1, and integrability is part of the conclusion. $\mathbf x_0$ may coincide with a design point.

A complete development needs: the law of linear images of independent Gaussian vectors; independence of the sample mean and the residual vector for i.i.d. normal samples (Cochran's theorem); the chi-squared law of the normal sample variance; the strong law of large numbers for the sample variance; and a measurability and boundedness argument for the random weight vector. The normal-sampling results are reusable well beyond this mission and are welcome as separate contributions.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2):371–382, 2010. https://doi.org/10.1287/opre.1090.0754
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8

Motivation

Many problems in statistics, signal processing and machine learning minimize the sum of a smooth function and a convex but nonsmooth one: $\ell_1$-regularized least squares (the Lasso), bound-constrained problems, and group-sparse regression are standard examples. Tseng and Yun (Math. Program. Ser. B 117 (2009) 387–423) proposed the coordinate gradient descent (CGD) method for this class and proved that it converges linearly under a local Lipschitzian error bound, Assumption 2(a) of their paper. An assumption is only as useful as the list of problems known to satisfy it. Section 6 of the paper supplies that list, and this mission formalizes it.

Error bounds of this kind have a history in smooth constrained optimization. Luo and Tseng proved them for quadratic objectives over polyhedral sets (Luo–Tseng 1992, SIAM J. Optim.), for strongly convex functions composed with a linear map (Luo–Tseng 1992, SIAM J. Control Optim.), and for dual functionals of linearly constrained strictly convex programs (Luo–Tseng 1993, Math. Oper. Res.), and used them to derive linear rates for feasible descent methods (Luo–Tseng 1993, Ann. Oper. Res.). Tseng and Yun carry these results over to the nonsmooth problem (1) by reformulating it as a smooth problem over the epigraph of the nonsmooth part.

Setting

Throughout, $\mathbb R^n$ carries the Euclidean norm $\|\cdot\|$. The problem is

$$\min_x\; F_c(x) = f(x) + cP(x), \qquad (1)$$

where $c > 0$, $P:\mathbb R^n \to (-\infty,\infty]$ is proper, convex and lower semicontinuous with effective domain $\operatorname{dom}P = \{x \mid P(x) < \infty\}$, and $f$ is continuously differentiable on an open set containing $\operatorname{dom}P$.

The residual at $x \in \operatorname{dom}P$ is

$$d_I(x) = \arg\min_d \Big\{\nabla f(x)^\top d + \tfrac12\|d\|^2 + cP(x+d)\Big\},$$

the unique minimizer of a strongly convex function. A point $x \in \operatorname{dom}P$ is stationary if the one-sided directional derivative satisfies $F_c'(x;d) \ge 0$ for every direction $d$; $\bar X$ denotes the set of stationary points and $\operatorname{dist}(x,\bar X)$ the distance to it. A point is stationary exactly when its residual vanishes.

Assumption 2(a) asks that $\bar X \neq \emptyset$ and that for every level $\zeta$ there be constants $\tau, \epsilon > 0$ with

$$\operatorname{dist}(x,\bar X) \le \tau\|d_I(x)\| \quad\text{whenever } F_c(x) \le \zeta,\ \|d_I(x)\| \le \epsilon.$$

Writing $\operatorname{epi}P = \{(x,\xi) \mid P(x) \le \xi\}$, problem (1) is equivalent to the smooth problem $\min\{f(x) + c\xi \mid (x,\xi) \in \operatorname{epi}P\}$ (41). Its projection residual at $(x,\xi)$ is the optimal solution $(\tilde d,\tilde\delta)$ of

$$\min_{(d,\delta)}\Big\{\nabla f(x)^\top d + \tfrac12\|d\|^2 + \tfrac12\delta^2 + c\delta \;\Big|\; (x+d,\xi+\delta) \in \operatorname{epi}P\Big\}. \qquad (42)$$

$P$ is polyhedral if $\operatorname{epi}P$ is the solution set of finitely many linear inequalities in $(x,\xi)$; $f$ is quadratic if $f(x) = \tfrac12 x^\top Ax + b^\top x + c_0$ with $A$ symmetric, not necessarily positive semidefinite. The four problem classes are:

• C1: $f$ quadratic, $P$ polyhedral;
• C2: $f(x) = g(Ex) + q^\top x$ with $g$ strongly convex and differentiable on $\mathbb R^m$, $\nabla g$ Lipschitz, $P$ polyhedral;
• C3: $f(x) = \max_{y\in Y}\{(Ex)^\top y - g(y)\} + q^\top x$ with $Y$ polyhedral and $g$ as in C2, $P$ polyhedral;
• C4: $f$ strongly convex with $\nabla f$ Lipschitz on $\operatorname{dom}P$ (condition (22)).

Formalization targets

Goal: Theorem 4 (p. 412)

$$\Big(\bar X \neq \emptyset \ \wedge\ (\mathrm{C1} \vee \mathrm{C2} \vee \mathrm{C3})\Big) \ \vee\ \mathrm{C4} \;\Longrightarrow\; \text{Assumption 2(a)}.$$

Under C4, nonemptiness of $\bar X$ is part of the conclusion. The constants $\tau, \epsilon$ are existential, so the goal does not depend on any particular estimate.

Milestone: Lemma 6 (p. 410)

If $P$ is Lipschitz on $\operatorname{dom}P$ with constant $K$, there is $\kappa > 0$ depending only on $K$ with

$$\|(\tilde d,\tilde\delta)\| \le \kappa\,\|d_I(x)\| \qquad (x \in \operatorname{dom}P,\ \xi = P(x)).$$

Milestone: Lemma 7 (p. 411)

If $\bar X \neq \emptyset$ and C1, C2 or C3 holds, then for every $\zeta$ there are $\tau', \epsilon' > 0$ with

$$\operatorname{dist}(x,\bar X) \le \tau'\|(\tilde d,\tilde\delta)\| \quad\text{whenever } F_c(x) \le \zeta,\ \|(\tilde d,\tilde\delta)\| \le \epsilon', \qquad (44)$$

where $(\tilde d,\tilde\delta)$ solves (42) with $\xi = P(x)$.

A further item, not a milestone, states the remark after Assumption 2 (p. 404) that Assumption 2(b) (stationary points with different objective values are uniformly separated) holds whenever $f$ is convex.

Significance

Theorem 4 is what makes the linear convergence theorem of the paper (Theorem 2, the subject of the second mission of this series) applicable. Through C1 it covers every problem with a quadratic $f$ and a polyhedral $P$, including the Lasso and $\ell_1$-regularized or bound-constrained quadratic programs, without convexity of $f$. Through C2 it covers losses of the form $g(Ex)$ with $g$ strongly convex, such as least squares $\tfrac12\|Ex - b\|^2$, where $f$ itself is not strongly convex because $E$ may have a nontrivial kernel. These error bounds became a standard tool for linear rates of proximal and coordinate methods without strong convexity.

The paper's result is proved, but Lemma 7 is proved by citation: it applies three error bounds of Luo and Tseng to the reformulation (41). A complete formalization therefore requires those error bounds for smooth problems over polyhedral sets, which are not in Mathlib. To our knowledge none of the results of this mission has a machine-checked proof.

Difficulty

Lemma 6 and the C4 case of Theorem 4 are short inequality arguments. The difficulty is Lemma 7. The natural first idea, proving an error bound from strong convexity, fails for C1–C3: $f$ may be nonconvex (C1) or have a degenerate Hessian (C2, C3), and the objective $f(x) + c\xi$ of (41) is never strongly convex in $(x,\xi)$. What replaces strong convexity is the polyhedral structure of $\operatorname{epi}P$, and the cited Luo–Tseng error bounds that exploit it are substantial results in their own right, each of which has to be established for a projection residual over a general polyhedron in $\mathbb R^{n+1}$.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n), with indices $0,\dots,n-1$. $P$ is the pair $(D, P)$, its effective domain and its finite values there, with "proper, convex, lsc" given by the published ProxNewton.Inexact.IsProperClosedConvex. The value $+\infty$ is never computed: every statement quantifies over $x \in D$, and membership in $\operatorname{epi}P$ is x + d ∈ D ∧ P (x + d) ≤ ξ + δ. $d_I(x)$ is a chosen minimizer of its subproblem, which at $x \in D$ is the paper's unique one. The solutions of (42) enter through a predicate, and every bound is asserted for every optimal solution. Stationarity is the liminf form of $F_c'(x;d) \ge 0$. $\operatorname{dist}$ is the infimum distance. Assumption 2(a) quantifies over every real $\zeta$, which agrees with the paper's "$\zeta \ge \min F_c$" because the condition is vacuous below $\inf F_c$. In Lemma 6 the constant $\kappa$ is quantified before the dimension and the data, so it depends on the Lipschitz constant only.

Explicit choices relative to the printed text:

• In C4, strong convexity and (22) are required on $\operatorname{dom}P$ only, since $f$ is only assumed smooth near $\operatorname{dom}P$. This hypothesis is weaker than the paper's, so the formal theorem is at least as strong.
• In C2 and C3, $E$ is a continuous linear map $\mathbb R^n \to \mathbb R^m$ (equivalently an $m \times n$ matrix). Strong convexity has a positive modulus and $\nabla g$ is globally Lipschitz.
• In C3, the maximum must be attained at every $x$, which forces $Y \neq \emptyset$, as the paper's formula presumes.
• Assumption 2(b) for convex $f$ is stated with convexity on $\operatorname{dom}P$.

A polyhedrality notion that admits only affine or constant $P$, an Assumption 2(a) whose constants are chosen after $x$, and a Lemma 7 whose residual is $d_I(x)$ instead of the solution of (42) would all trivialize or change the statements. The encoding rules out each of them: polyhedral means a finite system of linear inequalities in $(x,\xi)$, the constants precede $x$, and Lemma 7 is stated for (42).

Welcome contributions: proofs of Lemma 6 and of the C4 case; Lipschitz continuity of polyhedral functions on their domain (Rockafellar–Wets, Example 9.35); the Luo–Tseng error bounds for affine variational inequalities and for composite strongly convex objectives over polyhedra, which can be reused well beyond this mission.

Selected references

• P. Tseng, S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Math. Program. Ser. B 117 (2009) 387–423. https://doi.org/10.1007/s10107-007-0170-0
• Z.-Q. Luo, P. Tseng, Error bounds and the convergence analysis of matrix splitting algorithms for the affine variational inequality problem, SIAM J. Optim. 2 (1992) 43–54. https://doi.org/10.1137/0802004
• Z.-Q. Luo, P. Tseng, On the linear convergence of descent methods for convex essentially smooth minimization, SIAM J. Control Optim. 30 (1992) 408–425. https://doi.org/10.1137/0330025
• Z.-Q. Luo, P. Tseng, On the convergence rate of dual ascent methods for linearly constrained convex minimization, Math. Oper. Res. 18 (1993) 846–867. https://doi.org/10.1287/moor.18.4.846
• Z.-Q. Luo, P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Ann. Oper. Res. 46 (1993) 157–178. https://doi.org/10.1007/BF02096261
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3

Motivation

In the periodic-review inventory model with lost sales, demand that cannot be met from stock is lost rather than backordered. With a positive order lead time $L$ the state of the system is a vector of length $L$: on-hand stock plus the $L-1$ orders still in the pipeline. This model is standard in retail, where a customer facing an empty shelf buys elsewhere. Its optimal policy is not a base-stock policy and depends on the whole pipeline vector, which makes it much harder to analyse than the backorder model.

A basic question for any stochastic inventory model is how its optimal cost responds to uncertainty in demand. For the backorder model it is classical that more variable demand costs more (Song 1994). Zipkin's paper (Oper. Res. 56(4), 2008, 937–944) recasts the lost-sales model in a transformed state in which the optimal cost functions are L$^\natural$-convex. Its §5 uses this structure to prove the lost-sales analogue: if demand becomes more variable in the convex order, the optimal cost does not decrease. The paper presents this as new in the lost-sales setting.

Timeline of the structural results the mission rests on:

• 1958: Karlin and Scarf analyse the lost-sales model with lead time one and show that the optimal order decreases in the on-hand stock, with sensitivity less than one.
• 1969: Morton extends these monotonicity and bounded-sensitivity properties to general lead times and derives bounds on the optimal order.
• 2008: Zipkin reproves these properties through L$^\natural$-convexity in the transformed state (Theorem 4) and adds the parametric result of §5 (Theorem 11).

Setting

The order lead time is a positive integer $L$. The cost factors are the unit procurement cost $c$, the unit holding cost $\hat h$, the unit lost-sales penalty $p$ and the discount factor $\gamma$. Demands in different periods are independent and nonnegative, with a common law. The paper treats states and orders as continuous.

Transformed state. A state is a vector $v=(v_0,\dots,v_{L-1})$, where $v_l$ is the stock on hand plus all orders due to arrive $l$ or more periods hence, and $v_L=0$ by convention. So $v_0$ is the inventory position and $v_0-v_1$ is the stock on hand. The state space is

$$V=\{v\in\mathbb R^L:\ v_0\ge v_1\ge\cdots\ge v_{L-1}\ge0\}.$$

The action is $\zeta=-z\le0$, where $z\ge0$ is the order quantity, and $e$ is the all-ones vector. After demand $d$, the next state is

$$v_+=\big([v_0-v_1-d]^++v_1,\ v_2,\ \dots,\ v_{L-1},\ 0\big)-\zeta e.$$

Costs and recursion. The end-of-period holding and penalty cost is $\hat q(u)=\hat hu^++pu^-$, and its expectation given on-hand stock $y$ is $\hat q^0(y)=E[\hat q(y-d)]$. The optimal cost functions satisfy

$$\bar g_t(v,\zeta)=-\gamma^Lc\zeta+\hat q^0(v_0-v_1)+\gamma E[\bar f_{t+1}(v_+)],\qquad \bar f_t(v)=\min_{\zeta\le0}\bar g_t(v,\zeta),$$

with $\bar f_{T+L+1}=0$. This is the paper's recursion (1)–(2), written in the state $v$.

Program (4). Fix a demand $d$. Choosing the next stock level $v_+$, with $-d\le v_+-v_0\le0$ and $v_+\ge v_1$, gives

$$\bar\kappa_t(v,\zeta\mid d)=\min_{v_+}\Big\{\hat h(v_+-v_1)+p(v_+-v_0+d)+\gamma\bar f_{t+1}\big[(v_+,v_2,\dots,v_{L-1},0)-\zeta e\big]\Big\}.$$

Here $v_+-v_1$ is the leftover stock and $d-(v_0-v_+)$ the unmet demand.

Parametric demand. The demand $d(\phi)$ depends on a real parameter $\phi$, with law $\mu_\phi$, and $\bar f_t(v;\phi)$ is the optimal cost under $\mu_\phi$. A random variable $X$ is smaller than $Y$ in the convex order if $E\,g(X)\le E\,g(Y)$ for every convex $g:\mathbb R\to\mathbb R$ for which both expectations exist. This forces equal means and a variance that does not decrease.

Formalization targets

Goal: Theorem 11 (p. 941)

If $d(\phi)$ is increasing in $\phi$ with respect to the convex order, then for all $t$ and all $v\in V$,

$$\phi\le\phi'\ \Longrightarrow\ \bar f_t(v;\phi)\le\bar f_t(v;\phi').$$

Milestones, in attack order

1. Reformulation (proof of Theorem 4, p. 939). For $v\in V$, $\zeta\le0$, $d\ge0$, the stock level $v_+=[v_0-v_1-d]^++v_1$ is optimal in program (4):

$$\bar\kappa_t(v,\zeta\mid d)=\hat q(v_0-v_1-d)+\gamma\bar f_{t+1}(v_+).$$

2. Convexity of the optimal cost (Theorem 4, p. 939, with the remark on p. 938). $\bar f_t$ is convex on $V$ for every $t$.
3. The key step (proof of Theorem 11, p. 941). If $F$ is convex and nonnegative on $V$, then program (4) with continuation $F$ is jointly convex in $(v,\zeta,d)$ on $\{v\in V,\ \zeta\le0,\ d\ge0\}$.
4. The induction step (proof of Theorem 11, p. 941). If $\bar f_{t+1}(v;\phi)$ is nondecreasing in $\phi$ on $V$, then so is $\bar g_t(v,\zeta;\phi)$ for every $\zeta\le0$.

Significance

Theorem 11 is a comparative-statics statement: replacing demand by a mean-preserving spread cannot lower the optimal cost from any starting state, at any horizon. It justifies valuing variance reduction in lost-sales systems, for example forecasting or demand pooling, without solving the high-dimensional dynamic program. The induction pattern of the proof transfers to other parametric questions about the model; §6 of the paper notes that the parametric analysis remains valid in several extensions.

The paper proves Theorem 11 in a few lines on top of Theorem 4. No machine-checked proof of the result or of its structural prerequisites is known: neither the lost-sales recursion in the transformed state nor its convexity is formalized. A complete development checks the measure-theoretic steps the paper leaves implicit: that the optimal costs are finite, the existence of the expectations, and the passage from convexity in $d$ on $[0,\infty)$ to the convex order's test functions on $\mathbb R$.

Difficulty

The argument looks immediate: the costs are convex, so the convex order should apply directly. It does not. The function to which the convex order must be applied is the end-of-period cost as a function of demand. In recursion (1) that function is $\hat q(v_0-v_1-d)+\gamma\bar f_{t+1}(v_+(d))$, where the next state depends on $d$ through the kink $[v_0-v_1-d]^+$. Convexity of $\bar f_{t+1}$ alone does not make this composition convex in $d$. The paper avoids the issue by passing to program (4), where the amount of demand filled is a decision. That requires two facts: that selling as much as possible is optimal (milestone 1), and that the program is jointly convex in state, action and demand (milestone 3), which needs convexity of $\bar f_{t+1}$ on all of $V$ (milestone 2). A second difficulty is analytic: the optimal costs are defined by infima and expectations, and one must show that these are finite and integrable before any order comparison applies.

Formalization scope

• Model. States are Fin L → ℝ with 0-based indices matching the paper. $v_L=0$ is the helper vext. $V$ is the set of antitone, nonnegative vectors.
• Time. Data are stationary, so the optimal cost is indexed by the number of periods to go $k$, with $\bar f=0$ at $k=0$. The paper's "for all $t$" is "for all $k$".
• Demand. One-period demand has law $\mu$, a probability measure on $\mathbb R$ with $\mu((-\infty,0))=0$ and finite mean. The finite mean is not in the paper; it is added so that $\hat q^0$ is finite.
• Parametric family. A family $\mu_\phi$, $\phi\in\mathbb R$, with every hypothesis imposed for every $\phi$. Every object depending on demand takes the law as an argument, so $\bar f_t(v;\phi)$ is fbar … (μ φ) k v.
• Costs. "Unit cost" and "discount rate" are read as $c,\hat h,p\ge0$ and $0<\gamma\le1$.
• Minima. "Min" over orders is an infimum over $\zeta\le0$, and program (4) is an infimum over the interval $[\max(v_0-d,v_1),v_0]$. Lean's infimum is $0$ on sets that are unbounded below or empty. Nonnegative costs and demand keep the infima genuine, and the generic key-step milestone assumes the continuation is nonnegative.
• Expectations. Bochner integrals.
• Convex order. The platform definition StochasticOrders.Convex.ConvexOrder (Shaked–Shanthikumar (3.A.1)), applied to the identity on $(\mathbb R,\mu_\phi)$. Equal means are a consequence and are not assumed.
• Convexity milestone. Stated for the model's $\bar f$ only, not as "L$^\natural$-convex implies convex" in general, which fails without regularity.

Ruled out. Value functions that collapse to a junk constant would make Theorem 11 trivially true. This is excluded because the recursion is defined, not assumed, and every cost and demand hypothesis is imposed for every $\phi$. Replacing the convex order by "variance nondecreasing", by the increasing convex order, or by test functions convex only on $[0,\infty)$ changes the theorem and is not accepted.

Infrastructure and contributions. A complete development needs:

• finiteness and integrability of the optimal costs;
• convexity of partial infima over convex fibers;
• monotone extension of a function convex on $[0,\infty)$ to all of $\mathbb R$.

The model file is shared in shape with the series' missions 1 (L$^\natural$-convexity) and 2 (policy bounds), and its convexity lemmas are reusable for other parametric comparisons of the lost-sales model. Proofs of any milestone, and alternative arguments that bypass program (4), are welcome.

Selected references

• P. Zipkin, On the structure of lost-sales inventory models, Operations Research 56(4), 937–944, 2008. https://doi.org/10.1287/opre.1070.0482
• S. Karlin, H. Scarf, Inventory models of the Arrow–Harris–Marschak type with time lag, in Studies in the Mathematical Theory of Inventory and Production, Chapter 10, Stanford University Press, 1958.
• T. E. Morton, Bounds on the solution of the lagged optimal inventory equation with no demand backlogging and proportional costs, SIAM Review 11, 572–576, 1969 (as cited in Zipkin 2008).
• J.-S. Song, The effect of leadtime uncertainty in a simple stochastic inventory model, Management Science 40(5), 603–613, 1994. https://doi.org/10.1287/mnsc.40.5.603
• D. Stoyan, Comparison Methods for Queues and Other Stochastic Models, Wiley, 1983.
• M. Shaked, J. G. Shanthikumar, Stochastic Orders, Springer, 2007. https://doi.org/10.1007/978-0-387-34675-5

Motivation

Network design problems often ask for a cheap spanning tree in which no vertex is overloaded: in multicast and overlay networks the degree of a node bounds the number of copies it must forward, and in physical networks it bounds the number of ports. The minimum bounded degree spanning tree problem asks for a minimum-cost spanning tree whose degrees respect given bounds. Deciding whether a graph has a spanning tree of maximum degree 2 is already the Hamiltonian path problem, so exact solutions are out of reach, and the question becomes how little the bounds and the cost must be relaxed.

Some problems also impose lower bounds: a vertex that must serve as a hub, or a leaf-avoidance requirement, asks for degree at least $A_v$. Singh and Lau (STOC 2007) showed that the iterative rounding method they used for upper bounds extends to both kinds at once: a spanning tree of cost at most the LP optimum whose degree at every vertex lies in $[A_v-1,\,B_v+1]$.

Timeline. Fürer and Raghavachari (1994) found a spanning tree of maximum degree at most $\Delta^*+1$ in the unweighted case. For costs, Goemans (FOCS 2006) obtained cost at most OPT with degrees at most $B_v+2$. Singh and Lau (STOC 2007) improved this to $B_v+1$, which is the best possible additive violation unless P = NP, and gave the extension to lower and upper bounds formalized here. The journal version appeared in J. ACM 62(1), 2015.

Setting

Let $G=(V,E)$ be a finite simple graph with edge costs $c_e\in\mathbb R$; no sign or triangle inequality is assumed. Let $F$ be a forest on $V$ with no edge in common with $E$. A set $H\subseteq E$ is an $F$-tree if $H\cup F$ is a spanning tree of $V$, and $d_H(v)$ is the number of edges of $H$ at $v$. Integer lower bounds $A_v$ are given on $U\subseteq V$ and integer upper bounds $B_v$ on $W\subseteq V$. The minimum bounded degree connecting tree problem asks for a cheapest $F$-tree with $A_v\le d_H(v)$ on $U$ and $d_H(v)\le B_v$ on $W$; with $F=\varnothing$ it is the spanning-tree problem.

The supernodes are the vertex sets of the components of $(V,F)$, isolated vertices included, and $\mathcal I(F)$ is the family of unions of supernodes. For $S\subseteq V$, $E(S)$ and $F(S)$ are the edges with both endpoints in $S$, $\delta(S)$ the edges with exactly one endpoint in $S$, and $x(D)=\sum_{e\in D}x_e$. The linear relaxation LP-MBDCT$(G,\mathcal A,\mathcal B,U,W,F)$ minimizes $\sum_e c_ex_e$ subject to

$$x(E(V))=|V|-|F(V)|-1,\quad x(E(S))\le|S|-|F(S)|-1\ (S\in\mathcal I(F)),\quad A_v\le x(\delta(v))\ (v\in U),\quad x(\delta(v))\le B_v\ (v\in W),\quad x\ge0.$$

MBDCT Algorithm2 (Figure 5 of the paper) repeats: if $F$ is a spanning tree, stop; otherwise take a basic optimal solution $x^*$, delete the edges with $x^*_e=0$, move one edge with $x^*_e=1$ (if any) into $F$ and lower $A$ and $B$ by one at its endpoints, and otherwise remove from $U$ and $W$ one vertex whose support degree is at most two.

Formalization targets

Goal: Theorem 5.2

For every well-formed instance whose LP is feasible, MBDCT Algorithm2

1. has a terminating run;
2. makes at most $|V|-|F|-1+|U\cup W|$ iterations on every run;
3. returns only $F$-trees $H\subseteq E$ with

$$c(H)\le\sum_{e\in E}c_ex_e\ \text{for every feasible }x,\qquad A_v-1\le d_H(v)\ (v\in U),\qquad d_H(v)\le B_v+1\ (v\in W).$$

The statement fixes no constants beyond the paper's $\pm1$, and it is uniform over every choice of basic optimal solution, 1-edge and removed vertex.

Milestones

In the order the proof uses them: Lemma 5.3 (a basic solution is determined by a laminar family of tight sets and tight lower and upper degree rows, with $|E^*|=|\mathcal L|+|T_U|+|T_W|$); Claims 5.4, 5.7, 5.8 and 5.9 (special supernodes, cut sizes, the value $x^*(D(S))=r-1$ on the edges between the $r$ members of $S$, and when a set is special); Lemma 5.6 (every set of $\mathcal L$ keeps at least three tokens, exactly three only if it is special or $V$); and Lemma 5.1 (a basic solution has an edge with $x^*_e=1$ or a vertex of $U\cup W$ with support degree two).

Significance

The theorem gives, in one algorithm, a tree that is no more expensive than the LP lower bound and whose degrees miss both bounds by at most one. For the pure upper-bound problem it implies the $(1,B_v+1)$ guarantee, which cannot be improved to $B_v$ unless P = NP, since a Hamiltonian path is the case $B_v\equiv 2$. With lower bounds it also covers prescribed minimum degrees.

The result is proved in the paper; no machine-checked proof of it, or of any iterative rounding analysis, is known to exist. No platform item treats bounded-degree spanning trees or iterative rounding; the nearest related item is the open Williamson–Shmoys minimum-degree spanning tree local-search goal, which concerns a different, unweighted problem. Formalizing this mission requires the extreme-point structure of the spanning tree polytope with degree constraints (uncrossing to a laminar family), the token counting argument, and the induction over the algorithm's recursion, all of which are reusable for other iterative rounding results.

Difficulty

The obvious argument rounds an optimal LP solution once. It fails: an optimal vertex may be entirely fractional, and no single rounding keeps both the cost and the degrees. The algorithm instead re-solves the LP after each change, and its correctness rests on Lemma 5.1, that every extreme point of the current LP has an integral edge or a vertex with support degree two. That lemma is a counting statement about extreme points: the rank bound $|E^*|=|\mathcal L|+|T_U|+|T_W|$ from a laminar basis must be contradicted by a token distribution. With upper bounds only, every set of the laminar family can collect four tokens. With lower bounds a set may collect only three, and the proof must characterize those sets (special sets, Definition 5.5) and show by linear independence that the remaining cases cannot occur. A second difficulty is that the guarantee concerns the final tree after many iterations, so the degree accounting must survive the changes of $A$, $B$, $U$, $W$ and $F$ along the run.

Formalization scope

Vertices form a Fintype; edges are elements of Sym2 V with no loops, and $E$, $F$ are finite sets of edges with $F$ acyclic and disjoint from $E$. LP vectors are functions on Sym2 V that vanish off $E$. Right-hand sides are computed in $\mathbb R$, degree bounds are integers, and the subtour rows range over nonempty $S\in\mathcal I(F)$ (at $S=\varnothing$ the literal row $0\le-1$ would make the LP empty). A basic solution is a feasible $x$ that is the only vector on $E$ satisfying with equality every constraint tight at $x$, including the tight rows $x_e=0$. Supernodes are never contracted: the paper's contraction of a supernode with one active vertex is a proof device, and all of §5.1's notions are stated on the original instance.

The paper's loose phrases are made explicit as follows.

• "Polynomial time" becomes the iteration bound; the LP is solved by an oracle that returns any basic optimal solution. The ellipsoid method and the separation oracle are out of scope.
• "The algorithm returns" becomes three claims: a run exists, every run is bounded, and every returned set satisfies the guarantee. The guarantee alone would hold vacuously for a relation that can get stuck.
• "Cost at most the cost of the optimal LP solution" becomes $c(H)\le c\cdot x$ for every feasible $x$.
• "Distribute the tokens" (Lemma 5.6) becomes a counting inequality on the surplus of tokens.
• Claim 5.9's "contains exactly three special members" means that $S$ has exactly three members, all of them special.
• Figure 5's Step 4 is guarded by "no edge was picked in Step 3", as the paper's text under Lemma 5.1 states.
• "$F$ is not a spanning tree" and LP feasibility are explicit hypotheses where the paper leaves them implicit.

Two trivializations are ruled out. Existence of a cheap tree with degrees in $[A_v-1,B_v+1]$ is not the target: the statement is about the algorithm's outputs and quantifies over all of them. And a run that never terminates does not satisfy the goal, which requires a run to exist and every run to be bounded.

Contributions are welcome at every level: the laminar uncrossing argument (Lemma 5.3), which is reusable for any spanning-tree LP with degree rows; the counting claims; and the invariants of the recursion.

Selected references

• M. Singh, L. C. Lau, Approximating minimum bounded degree spanning trees to within one of optimal, STOC 2007, pp. 661–670. https://doi.org/10.1145/1250790.1250887
• M. Singh, L. C. Lau, Approximating minimum bounded degree spanning trees to within one of optimal, J. ACM 62(1), 2015. https://doi.org/10.1145/2629366
• M. X. Goemans, Minimum bounded degree spanning trees, FOCS 2006, pp. 273–282. https://doi.org/10.1109/FOCS.2006.48
• M. Fürer, B. Raghavachari, Approximating the minimum-degree Steiner tree to within one of optimal, J. Algorithms 17(3), 1994, pp. 409–423. https://doi.org/10.1006/jagm.1994.1042

Motivation

Many problems in statistics, signal processing and machine learning minimize the sum of a smooth loss and a convex, nonsmooth, separable regularizer: smooth optimization with $\ell_1$-regularization, bound-constrained optimization, the group Lasso and, through duality, support vector regression are special cases listed in §1 of the paper discussed below. For such problems, methods that update one coordinate, or one small block of coordinates, at a time are often the method of choice, because each update is cheap and exploits separability.

Paul Tseng and Sangwoon Yun (Math. Program. Ser. B 117 (2009) 387–423) introduced the coordinate gradient descent (CGD) method for this class: at each iteration it minimizes a quadratic model of the smooth part plus the exact nonsmooth part over a chosen block of coordinates, then takes an Armijo step. Their paper proves global convergence (Theorem 1) and, under a local error bound, linear convergence (Theorems 2 and 3). This mission is the second of a three-mission series on the paper; it targets the linear convergence theorem for the cyclic (Gauss–Seidel-type) choice of blocks.

The analysis follows Luo and Tseng's error-bound framework for smooth constrained problems (Ann. Oper. Res. 46 (1993) 157–178), which established linear rates for feasible descent methods without strong convexity. Tseng and Yun extend that framework to a nonsmooth convex term.

Setting

Let $\Re^n$ carry the Euclidean norm $\|\cdot\|$, let $\mathcal N=\{1,\dots,n\}$, and for $\mathcal J\subseteq\mathcal N$ let $x_{\mathcal J}$ be the subvector of coordinates in $\mathcal J$. The problem (1) is

$$\min_x\ F_c(x)=f(x)+cP(x),$$

where $c>0$, $P:\Re^n\to(-\infty,\infty]$ is proper, convex and lower semicontinuous, and $f$ is continuously differentiable on an open set containing $\operatorname{dom}P$.

Given $x\in\operatorname{dom}P$, a nonempty block $\mathcal J$ and a symmetric positive definite matrix $H$, the direction (6) is

$$d_H(x;\mathcal J)=\arg\min_d\Big\{\nabla f(x)^Td+\tfrac12d^THd+cP(x+d)\ \Big|\ d_j=0\ \forall j\notin\mathcal J\Big\}.$$

The CGD method starts at $x^0\in\operatorname{dom}P$, chooses at each iteration $k$ a block $\mathcal J^k$ and a matrix $H^k$, computes $d^k=d_{H^k}(x^k;\mathcal J^k)$, and sets $x^{k+1}=x^k+\alpha^kd^k$. The Armijo rule takes $\alpha^k$ as the largest element of $\{\alpha^k_{\rm init}\beta^j\}_{j\ge0}$ with

$$F_c(x^k+\alpha^kd^k)\le F_c(x^k)+\alpha^k\sigma\Delta^k,\qquad \Delta^k=\nabla f(x^k)^Td^k+\gamma\,d^{kT}H^kd^k+cP(x^k+d^k)-cP(x^k),$$

with $0<\beta,\sigma<1$ and $0\le\gamma<1$. The restricted Gauss–Seidel rule (12) asks for indices $0=t_0<t_1<\cdots$ (the set $\mathcal T$) such that the blocks used during each cycle $t_i\le k<t_{i+1}$ are pairwise disjoint and cover $\mathcal N$. $P$ is block-separable with respect to $\mathcal J$ if $P(x)=P_{\mathcal J}(x_{\mathcal J})+P_{\mathcal J^C}(x_{\mathcal J^C})$.

The analysis uses three assumptions: $\nabla f$ is $L$-Lipschitz on $\operatorname{dom}P$ (22); Assumption 1, $\bar\lambda I\succeq H^k\succeq\underline\lambda I$ with $0<\underline\lambda\le\bar\lambda$; and Assumption 2. Writing $\bar X$ for the set of stationary points of $F_c$ and $d_I(x)$ for the full direction with $H=I$, Assumption 2 says (a) $\bar X\ne\emptyset$ and the local Lipschitzian error bound $\operatorname{dist}(x,\bar X)\le\tau\|d_I(x)\|$ holds on every level set $\{F_c\le\zeta\}$ wherever $\|d_I(x)\|\le\epsilon$; (b) stationary points with different objective values are at least a fixed distance $\delta>0$ apart.

Formalization targets

Goal: Theorem 2(b)

Under (22), Assumptions 1 and 2, the restricted Gauss–Seidel rule, block-separability of $P$ with respect to every $\mathcal J^k$, and Armijo steps with $\sup_k\alpha^k_{\rm init}\le1$ and $\inf_k\alpha^k_{\rm init}>0$:

$$\{F_c(x^k)\}\downarrow-\infty\quad\text{or}\quad\big(\{F_c(x^k)\}_{\mathcal T}\ \text{Q-linear}\ \text{and}\ \{x^k\}_{\mathcal T}\ \text{R-linear}\big).$$

No rate constant is fixed: the goal asserts the shape of the convergence, so it is not invalidated by sharper constants.

Milestones

Lemma 4 (Hölder dependence of the subproblem's solution on its linear term), Lemma 5(a) (a per-block variational inequality), Lemma 5(b) (an explicit stepsize at which the Armijo test passes), Theorem 1(f) (Armijo stepsizes bounded away from zero, $\Delta^k\to0$, $d^k\to0$), and Theorem 2(a) (the residual $\|d_I(x^k)\|$ at the start of a cycle is bounded by the step lengths of the cycle).

Significance

Theorem 2(b) gives a linear rate for a block coordinate method on a nonsmooth, possibly nonconvex composite problem without strong convexity: the rate follows from the error bound, which holds for a large class of problems studied in §6 of the paper (the third mission of this series).

The theorem is proved in the paper; none of it has been machine-checked, as far as a search of the Prove2Me library shows. Formalizing it requires a reusable treatment of extended-valued convex functions through their effective domains, of the coordinate-restricted proximal subproblem, and of the Armijo rule for composite objectives. The milestones Lemma 4 and Lemma 5 are self-contained facts about this subproblem and are reusable for other proximal coordinate methods.

Difficulty

The classical route to a linear rate, as in Luo and Tseng's smooth analysis, derives $F_c(x^{k+1})-\bar\upsilon\le\tau'\|x^{k+1}-x^k\|^2$ from the error bound, where $\bar\upsilon$ is the limit value. With a nonsmooth $P$ this inequality is not available: the authors say so on p. 405, and work with $-\Delta^k$ in place of the quadratic term. A second obstacle is that the error bound controls the full residual $d_I(x^k)$, while each iteration only computes a block direction at a different point; relating the two over a cycle needs the disjointness of the blocks in the restricted rule and block-separability of $P$. Neither the global analysis nor the error bound alone gives the rate.

Formalization scope

$\Re^n$ is EuclideanSpace ℝ (Fin n), with 0-based coordinates. The extended-valued $P$ is encoded as a pair: its effective domain $D=\operatorname{dom}P$ and its finite values on $D$, with "proper convex lsc" given by the published definition ProxNewton.Inexact.IsProperClosedConvex; $F_c$ is never evaluated off $D$, and every statement carries membership in $D$ explicitly. The direction $d_H(x;\mathcal J)$ is the unique minimizer of (6) when $x\in D$ and $H\succ0$ (chosen by Classical.epsilon), and runs of the method carry the minimizer property itself. The Armijo rule takes the first admissible exponent $j$. $\mathcal T$ is a strictly increasing map $t$ with $t_0=0$. Stationarity is $F_c'(x;d)\ge0$ for all $d$, in liminf form. In Assumption 2, $\operatorname{dist}$ is the infimum distance and "for any $\zeta\ge\min F_c$" is "for every real $\zeta$" (equivalent, since the condition is empty below $\inf F_c$). Q-linear convergence includes convergence of the sequence; R-linear convergence is $\|x^{t_i}-\bar x\|\le Cq^i$ with $q<1$. The rates in the goal are along $\mathcal T$, as the paper states.

Explicit choices relative to the printed text:

• Theorem 2(a) is corrected. As printed, $\|d_I(x^k)\|\le\sup_j\alpha^j\,C\,r^k$ fails for small constant stepsizes and fails without block-separability (counterexamples in the item's statement). The milestone states $\|d_I(x^k)\|\le\max\{1,\sup_j\alpha^j\}\,C\,r^k$ under block-separability of $P$ with respect to every $\mathcal J^k$ (both are what the paper's proof gives and what Theorem 2(b) supplies), with the iterates in $\operatorname{dom}P$, and with $C$ quantified before the problem data, so it depends only on $n,L,\underline\lambda,\bar\lambda$.
• Lemma 5(b)'s range $0\le\alpha\le\min\{1,2\underline\lambda(1-\sigma+\sigma\gamma)/L\}$ is stated as $0\le\alpha\le1$ and $\alpha L\le2\underline\lambda(1-\sigma+\sigma\gamma)$, so that $L=0$ gives $[0,1]$ as on paper rather than $\{0\}$.
• Lemma 5(a) holds for every decomposition of $P$ in (20), and $\bar x$ ranges over $\operatorname{dom}P_{\mathcal J}$ (outside it the left side is $-\infty$).
• "$\lim F_c(x^k)>-\infty$" is "the values $F_c(x^k)$ are bounded below"; "$\{F_c(x^k)\}\downarrow-\infty$" is "nonincreasing and tending to $-\infty$"; suprema and infima of stepsizes are explicit bounds.

A statement that replaces Assumption 2 by strong convexity, claims a rate along the whole sequence rather than along $\mathcal T$, or asserts the Q-linear contraction without convergence to a limit is not this theorem and is out of scope.

Lemma 1, the claim after (10), Lemma 3 and Theorem 1(a), which the proof also uses, are milestones of the first mission of this series and are not restated here. Theorem 3 (the same conclusion along the whole sequence for the Gauss–Southwell-q rule) is not included: its proof is omitted in the paper. Contributions welcome: proofs of the milestones, and general facts about the coordinate-restricted proximal subproblem (existence and uniqueness of $d_H(x;\mathcal J)$, optimality conditions) that the milestones rely on.

Selected references

• P. Tseng and S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Mathematical Programming Ser. B 117 (2009) 387–423. https://doi.org/10.1007/s10107-007-0170-0
• Z.-Q. Luo and P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Annals of Operations Research 46 (1993) 157–178. https://doi.org/10.1007/BF02096261
• J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970, Chap. 9 (Q- and R-linear convergence). https://doi.org/10.1137/1.9780898719468
• R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3