Optimization

Where in-process inventories should sit

A product made by a chain of $n$ operations can carry stock in front of every operation, between every pair of them, and as finished goods. Each such inventory point buffers demand uncertainty, and each costs money to hold. Kenneth F. Simpson, Jr. asked which of these points should carry stock at all, and how much, in In-Process Inventories, Operations Research 6 (1958) 863–873. His answer is that the cheapest configuration is an extreme one: every inventory point is either kept full or eliminated.

The paper is the origin of what is now called the guaranteed-service model of multi-echelon inventory theory.

• 1958: Simpson proves that the optimal service times lie at a vertex of the feasible polyhedron, for a serial line, and notes extensions to branching chains and to correlated demand.
• 2000: Graves and Willems (Optimizing Strategic Safety Stock Placement in Supply Chains, M&SOM 2(1)) extend the service-time formulation to spanning-tree networks and solve it by dynamic programming over the vertex structure.
• Textbooks (Snyder and Shen, Fundamentals of Supply Chain Theory, Wiley) present the serial result as the "all-or-nothing" property of guaranteed-service systems.

The model

There are $n$ operations in series. Operation $i$ has a processing time $T_i$. Behind each operation is an inventory point that promises to fill any order within a service time $S_i \ge 0$. The incoming service time $S_0$ (from raw materials) and the outgoing one $S_n$ (to the customer) are fixed constants; the $n - 1$ intermediate service times $S_1, \dots, S_{n-1}$ are the decision variables. Inventory $i$ must protect against demand over the time $S_{i-1} - S_i + T_i$, and under the paper's demand model its holding cost is proportional to the square root of that time. With constants $c$ and $r_1, \dots, r_n > 0$ the total inventory cost is equation (9):

$$f(S_1, \dots, S_{n-1}) = c + \sum_{i=1}^{n} r_i \sqrt{S_{i-1} - S_i + T_i}.$$

Its domain $D \subseteq \mathbb R^{n-1}$ is the set of points with every variable nonnegative and every radicand nonnegative:

$$D = \{ (S_1, \dots, S_{n-1}) : S_i \ge 0 \ (1 \le i \le n - 1),\ \ S_{i-1} - S_i + T_i \ge 0 \ (1 \le i \le n) \}.$$

$D$ is a polyhedron cut out by $2n - 1$ affine inequalities. A vertex of $D$ is an extreme point: a point of $D$ that is not the midpoint of two distinct points of $D$. An inventory with $S_i = S_{i-1} + T_i$ has zero average stock; one with $S_i = 0$ is full.

Formalization targets

Goal: THEOREM I (p. 869)

For any constants $c, S_0, S_n, T_i$ and $r_i > 0$ with $D \ne \emptyset$,

$$\exists\, v \in \operatorname{ext} D \quad\text{with}\quad f(v) \le f(S)\ \ \text{for all } S \in D.$$

This is the paper's sentence "the minimum of the function $f$ occurs at one of the vertices of $D$", read as existence of a vertex minimizer. It fixes no constant and no particular vertex.

Milestones

1. $D$ is convex and compact (pp. 868, 872).
2. $f$ is concave on $D$ when $r_i \ge 0$ (p. 868, "concave downward").
3. $f$ attains its minimum on a nonempty $D$ (Appendix, p. 872).
4. LEMMA I (p. 872): for at least one variable, no interior point of $D$ minimizes $f$.
5. THEOREM II (p. 869): if $S_j = S_{j-1} + T_j$, the cost equals the cost of the model with operations $j$ and $j+1$ coupled into one with processing time $T_j + T_{j+1}$.
6. LEMMA II (pp. 872–873): for at least two variables, every minimizer has at least two independent active boundary normals.
7. Meaning of the theorem (p. 869): when $S_0 = S_n = 0$ and $T_i > 0$, the vertices of $D$ are exactly the $2^{n-1}$ points with each $S_i \in \{0,\ S_{i-1} + T_i\}$.
8. Extension (p. 871): THEOREM I holds with $\sqrt{t}$ replaced by $t^k$ for $\tfrac12 \le k \le 1$.

The p. 870 rephrasing, "for an optimal system, the service time for every nonvanishing inventory is zero", is already proved on the platform as SupplyChainTheory.gs_all_or_nothing (Snyder–Shen Theorem 6.5) and is included as a reference item. Under the dictionary stage $i \leftrightarrow n + 1 - i$, $r_i = h_i k$, $SI_N = S_0$, that statement is Simpson's model case $S_n = 0$.

Significance

THEOREM I turns a continuous optimization over the polyhedron $D$ into a comparison over its finitely many vertices, and in Simpson's model these are the $2^{n-1}$ all-or-nothing configurations. Each vertex has a physical reading as a mode of operation: which inventories are kept full and which operations are tightly coupled (THEOREM II). The property is what makes guaranteed-service models tractable by enumeration or dynamic programming in later work.

The result is proved in the paper, by an induction on faces. It has not been machine-checked. A formal development adds a checked statement of the theorem for general constants $S_0, S_n, T_i$ (the paper's model fixes $S_0 = S_n = 0$), an exact account of which hypotheses are needed (the positivity of the $r_i$ is essential, see below), and a checked vertex count for the model case. The existing platform theorem gs_all_or_nothing covers the model case in the all-or-nothing form; THEOREM I with general constants and the vertex characterization are not yet on the platform.

Difficulty

The domain $D$ is given by inequalities, while the conclusion is about its extreme points, which are never listed: with general constants $S_0, S_n, T_i$ the vertices of $D$ have no closed form, and $D$ may be lower-dimensional or have coinciding facets. Any argument therefore has to connect an inequality description of a polyhedron with its extreme-point structure, which in the paper's own text is treated informally ("faces", "edges", "an $n$-fold repetition"). LEMMA I and LEMMA II are stronger than the existence statement: they exclude minimizers in the interior of $D$ and in the relative interior of a face, and they become false as soon as some $r_i$ vanishes, so they cannot follow from concavity alone. THEOREM II, though described on the page as "examine the definitions", requires an exact reindexing of the operations, processing times and cost coefficients after one variable is removed.

The sign hypothesis cannot be dropped: with $n = 2$, $r_1 = r_2 = -1$, $S_0 = 0$, $S_2 = 1$, $T_1 = 2$, $T_2 = 1$, the function $f(S_1) = -\sqrt{2 - S_1} - \sqrt{S_1}$ on $D = [0, 2]$ is convex and its only minimizer is the interior point $S_1 = 1$.

Formalization scope

All declarations live in SimpsonInv.Vertex. With $m = n - 1$ variables a point of $D$ is S : Fin m → ℝ, and Lean index $j$ holds the paper's $S_{j+1}$; extVec S0 Sn S i returns the paper's $S_i$ for every $i$, including the constant ends. The constants $T, r$ are sequences ℕ → ℝ, read only at indices $1, \dots, n$. The space Fin m → ℝ carries the product topology, which is the Euclidean one, so interior is the paper's interior. A vertex is an element of Set.extremePoints ℝ D; a minimizer is IsMinOn f D v with $v \in D$.

Added hypotheses, each disclosed in its item: $r_i > 0$ for $1 \le i \le n$ (from equation (8); weakened to $r_i \ge 0$ for concavity and dropped for existence and THEOREM II); $D \ne \emptyset$ where a minimum is asserted; at least one variable for LEMMA I and at least two for LEMMA II; $S_0 = S_n = 0$ and $T_i > 0$ only for the vertex count, as in the paper's model. The page's "$n-1$-dimensional polyhedron" is not imposed: the goal is stated for every nonempty $D$. "On the boundary of the faces" (LEMMA II) is encoded by at least two independent active inequality normals, including when inequalities coincide.

Real.sqrt returns $0$ on negative inputs, so $f$ is never compared outside $D$. Minimizing over all of $\mathbb R^{n-1}$, or replacing the vertex set by a hand-listed finite set of points, would trivialize or falsify the goal and is ruled out: the goal quantifies over Set.extremePoints ℝ D and compares only on $D$. THEOREM II compares the original cost with the same cost function applied to new data ($n - 1$ operations, $T'$, $r'$), not with a substituted copy of the original.

A complete development needs compactness and extreme-point facts for polyhedra given by inequalities and concavity facts for square roots and powers of affine maps; both are reusable well beyond this mission. Proofs of any milestone, and reusable lemmas on polyhedra given by inequalities, are welcome.

Selected references

• Kenneth F. Simpson, Jr., In-Process Inventories, Operations Research 6(6):863–873, 1958. https://doi.org/10.1287/opre.6.6.863
• Stephen C. Graves and Sean P. Willems, Optimizing Strategic Safety Stock Placement in Supply Chains, Manufacturing & Service Operations Management 2(1):68–83, 2000. https://doi.org/10.1287/msom.2.1.68.23267
• Lawrence V. Snyder and Zuo-Jun Max Shen, Fundamentals of Supply Chain Theory, 2nd ed., Wiley, 2019 (Chapter 6, guaranteed-service serial systems).

Motivation

A loss network is a model of a circuit-switched network: telephone networks, and more generally any system in which a request seizes several resources at once and is lost if any of them is unavailable. Exact loss probabilities in such networks are given by a product-form distribution over a state space whose size grows exponentially with the number of routes, so practitioners compute an approximation instead: the Erlang fixed point, also called the reduced-load approximation. It pretends that links block independently, computes the traffic offered to each link after thinning by the other links, and applies Erlang's single-link formula link by link. Reduced-load approximations of this kind recur throughout §§3–4 of Kelly 1991, where they appear as limits of networks with fixed and with alternative routing.

An approximation defined by a fixed point equation raises an immediate question: does the equation have exactly one solution? Under fixed routing it does, and Kelly's Theorem 3.7 (Ann. Appl. Probab. 1991, p. 338, citing Kelly 1986, reference [33] of the paper) proves it by identifying the fixed point with the optimum of a strictly convex program, the revised dual problem. The convex characterization is specific to fixed routing; the paper notes (p. 337) that for networks with dynamic routing and control nonuniquely determined behaviour can occur.

Setting

There are $J$ links, link $j$ having $C_j \ge 1$ circuits, and a finite set of routes $r$. A call on route $r$ uses $A_{jr} \in \mathbb Z_+$ circuits of link $j$, and calls on route $r$ arrive as a Poisson stream of rate $\nu_r > 0$.

Erlang's formula (1.1) gives the blocking probability of a single link with $C$ circuits offered Poisson traffic of intensity $\nu$:

$$E(\nu, C) = \frac{\nu^C}{C!}\Big[\sum_{n=0}^{C} \frac{\nu^n}{n!}\Big]^{-1}.$$

The Erlang fixed point equations (3.1)–(3.2) for the vector of link blocking probabilities $(E_1, \dots, E_J)$ are

$$E_j = E(\rho_j, C_j), \qquad \rho_j = (1 - E_j)^{-1} \sum_r A_{jr}\,\nu_r \prod_i (1 - E_i)^{A_{ir}}, \qquad j = 1, \dots, J.$$

The utilization function $U(y, C)$ is defined by the implicit relation (3.4),

$$U\big(-\log(1 - E(\nu, C)),\, C\big) = \nu\,\big(1 - E(\nu, C)\big), \qquad \nu \ge 0:$$

the mean number of busy circuits on a single link whose blocking probability is $1 - e^{-y}$.

The revised dual problem (3.5) is

$$\text{minimize} \quad \sum_r \nu_r \exp\Big(-\sum_j y_j A_{jr}\Big) + \sum_j \int_0^{y_j} U(z, C_j)\,dz \qquad \text{subject to } y \ge 0,$$

and its stationarity conditions (3.6) are

$$\sum_r A_{jr}\,\nu_r \exp\Big(-\sum_i y_i A_{ir}\Big) = U(y_j, C_j), \qquad j = 1, \dots, J.$$

The revised dual differs from the dual problem (2.3) of §2, whose optimum gives the limiting blocking probabilities of a large network, only in its last term: $\sum_j y_j C_j$ becomes $\sum_j \int_0^{y_j} U(z, C_j)\,dz$.

Formalization targets

Goal: Theorem 3.7

The revised dual problem (3.5) has an optimum $y$, every optimum equals $y$, and

$$E \in [0,1]^J \text{ solves (3.1)–(3.2)} \iff E_j = 1 - \exp(-y_j) \text{ for every } j.$$

The goal states the characterization, not only the existence and uniqueness of the fixed point: the solution is $1 - e^{-y}$ for the optimum $y$ of (3.5).

Milestones (p. 338, in the order the argument uses them)

1. (3.4) defines a function: $\nu \mapsto -\log(1 - E(\nu, C))$ is a continuous strictly increasing bijection of $[0,\infty)$ onto itself for $C \ge 1$, $U$ satisfies (3.4), and $U \ge 0$.
2. $y \mapsto U(y, C)$ is strictly increasing on $[0, \infty)$.
3. $y \mapsto \int_0^y U(z, C)\,dz$ is strictly convex, and so is the objective of (3.5) on $y \ge 0$.
4. (3.5) has exactly one optimum.
5. A vector $y \ge 0$ is optimal for (3.5) if and only if it satisfies (3.6).
6. A solution of (3.1)–(3.2) in $[0,1]^J$ has all $E_j < 1$; and for $y \ge 0$, $E_j = 1 - e^{-y_j}$ solves (3.1)–(3.2) if and only if $y$ satisfies (3.6).

Significance

The theorem turns a nonlinear system of $J$ equations into a smooth strictly convex minimization over the nonnegative orthant. Uniqueness of the fixed point follows at once, and the program gives a computational handle: any convex-minimization method computes the Erlang fixed point, without relying on the convergence of repeated substitution. The comparison with the dual problem (2.3) also explains the approximation: (3.6) relaxes the conditions (2.7), which force a link's carried traffic to equal its capacity before it can block, into a smooth relation in which blocking grows with carried traffic as Erlang's formula prescribes (Remark 3.8).

On Prove2Me, existence and uniqueness of the fixed point alone is already the Proved theorem KellyStochasticNetworks.erlang_fixed_point_unique (Kelly–Yudovina, Stochastic Networks, Theorem 3.20), referenced by this mission. Also referenced: KellyStochasticNetworks.erlang_strictMono (Proved: $\nu \mapsto E(\nu, C)$ and $\nu \mapsto \nu(1 - E(\nu, C))$ are strictly increasing for $C \ge 1$) and KellyStochasticNetworks.erlang_mean_busy_circuits (Proved: the mean number of busy circuits is $\nu(1 - E(\nu, C))$, the reading of $U$ as a utilization). What this mission adds is the revised dual problem itself, the utilization function $U$, and the identification of the fixed point with the optimum of (3.5); none of these is formalized on the platform.

Difficulty

The function $U$ is defined only implicitly, through the inverse of $\nu \mapsto -\log(1 - E(\nu, C))$, so every property of the objective of (3.5) passes through properties of Erlang's formula: monotonicity, continuity, and the limit $E(\nu, C) \to 1$ as $\nu \to \infty$. Differentiating $\int_0^{y} U$ needs continuity of $U$, which must be derived from the inverse.

The paper's sentence "strictly convex: it thus has a unique minimum" covers only uniqueness. Existence of a minimizer on the unbounded set $y \ge 0$ is a separate fact, a growth estimate for the integral terms, and it is part of milestone 4. The first term of (3.5) is not strictly convex when $A$ has rank less than $J$, so strict convexity must come from the separable integral terms. Finally, optimality over $y \ge 0$ yields only one-sided conditions at coordinates with $y_j = 0$; that (3.6) is an equality there uses $U(0, C) = 0$.

Formalization scope

Lean namespace KellyLossNetworks.RevisedDual. Links are Fin J, routes Fin R, the incidence matrix is A : Fin J → Fin R → ℕ, rates ν : Fin R → ℝ, capacities C : Fin J → ℕ, following the published Kelly–Yudovina definitions. Every statement assumes $\nu_r > 0$ and $C_j \ge 1$; the paper uses both without stating them, and a link with $C_j = 0$ makes (3.4) meaningless, since $E(\nu, 0) = 1$.

Erlang's formula is the published KellyStochasticNetworks.erlang and the fixed point equations (3.1)–(3.2) are the published KellyStochasticNetworks.ErlangFixedPoint, whose factor $(1 - E_j)^{-1}$ is Lean's total inverse; milestone 6 shows that no solution in $[0,1]^J$ reaches $E_j = 1$. Solutions are sought in $[0,1]^J$, as on p. 338.

$U$ is defined by choice: $U(y, C) = \nu(1 - E(\nu, C))$ for a $\nu \ge 0$ with $-\log(1 - E(\nu, C)) = y$, and $0$ if none exists. For $C \ge 1$, $y \ge 0$ the $\nu$ is unique, and that (3.4) holds is milestone 1, a theorem, not an axiom built into the definition. Values at $y < 0$ or $C = 0$ are placeholders that no statement uses. The integral in (3.5) is the interval integral $\int_0^{y_j}$. An optimum of (3.5) is a $y \ge 0$ whose objective value is at most that of every $z \ge 0$.

A trivializing formalization is ruled out: replacing $U$ by any closed form other than the one fixed by (3.4) (for instance $U(z, C) = C$, the fluid utilization of Remark 3.8) would change the problem and reduce (3.5) to the dual (2.3); and the goal asserts the identification $E = 1 - e^{-y}$ rather than restating the existence and uniqueness that is already Proved.

A complete development needs: inverse-function arguments for a strictly increasing continuous map of $[0, \infty)$, derivatives of interval integrals with continuous integrands, strict convexity of separable sums, and first-order optimality conditions for convex functions on the nonnegative orthant. The facts about Erlang's formula and about $U$ are reusable in any reduced-load analysis. Proofs of individual milestones are welcome independently.

Selected references

• F. P. Kelly, Loss networks, The Annals of Applied Probability 1(3):319–378, 1991. https://doi.org/10.1214/aoap/1177005872
• F. P. Kelly, Blocking probabilities in large circuit-switched networks, Advances in Applied Probability 18:473–505, 1986 (reference [33] of Kelly 1991; DOI not verified offline).
• F. P. Kelly and E. Yudovina, Stochastic Networks, Cambridge University Press, 2014 (the source of the Proved platform items referenced here; DOI not verified offline).

Motivation

A linear program often has two objectives that pull in different directions, for example cost and risk, or two policy goals whose relative weight is a matter of judgment. Gass and Saaty (NRLQ 1955) put the question in its simplest form: let every cost coefficient be an affine function $d_j + \lambda d'_j$ of one real parameter $\lambda$, and find, for every value of $\lambda$, an optimal solution. Their paper turned this into a procedure that runs the simplex method once and then walks along the $\lambda$-axis, changing one basis column at a time. This parametric objective function method became a standard chapter of linear programming texts and is the prototype of parametric and sensitivity analysis in optimization software.

Timeline.

• 1951: Dantzig publishes the simplex method (Cowles Commission Monograph No. 13, Ch. XXI; the paper's reference [1]), with the optimality test $z_j - c_j \le 0$ for a minimization problem.
• 1954: Gass and Saaty summarize the parametric procedure (J. Oper. Res. Soc. Amer. 2 (1954), 316–319).
• 1955: the present paper gives the procedure in detail, proves the THEOREM that makes it work, and reports computations on the National Bureau of Standards SEAC computer.

Setting

For each real $\lambda$ the problem is

$$\text{minimize } \sum_{j=1}^n (d_j + \lambda d'_j)\,x_j \quad\text{subject to}\quad x_j \ge 0\ (j = 1,\dots,n),\qquad \sum_{j=1}^n a_{ij}x_j = a_{i0}\ (i = 1,\dots,m).$$

Write $A = (a_{ij})$, an $m\times n$ real matrix with columns $A_1,\dots,A_n$, and $b = (a_{i0})$. A basis is a list of $m$ linearly independent columns $A_{B(1)},\dots,A_{B(m)}$; its basic solution $x$ has $x_j = 0$ off the basis and $Ax = b$, and it is a basic feasible solution when also $x \ge 0$. Every column has coordinates $y_{ij}$ in the basis, $A_j = \sum_i y_{ij}A_{B(i)}$. For a cost vector $c$ the simplex method uses $z_j = \sum_i y_{ij}c_{B(i)}$; for the cost $c = d + \lambda d'$ the number $z_j - c_j$ is affine in $\lambda$:

$$z_j - c_j = \alpha_j + \lambda\beta_j .$$

The basis is optimal at $\lambda$ when $\alpha_j + \lambda\beta_j \le 0$ for all $j$ (inequalities (3)). The critical values of the basis are

$$\underline\lambda = \max_{\beta_j<0}\Big(-\frac{\alpha_j}{\beta_j}\Big)\ (-\infty \text{ if all }\beta_j\ge 0),\qquad \bar\lambda = \min_{\beta_j>0}\Big(-\frac{\alpha_j}{\beta_j}\Big)\ (+\infty \text{ if all }\beta_j\le 0).$$

The paper assumes throughout that the problem is nondegenerate: no basic feasible solution has more than $n - m$ zero components.

Case A is the situation where the current basis is optimal for some $\lambda$, so that (3) is consistent. If $\bar\lambda$ is finite, attained at an index $s$ with $\beta_s > 0$, and some $y_{is} > 0$, the method brings $A_s$ into the basis and removes the column $A_r$, $r = B(\ell)$, chosen by the simplex ratio test; $y_{rs} > 0$. Call the new basis $B'$ and its basic feasible solution $x'$.

Formalization targets

Goal: the THEOREM (p. 41)

$$\bar\lambda = \min\big\{\lambda\in\mathbb R : x' \text{ minimizes } (d+\lambda d')^{\top}x \text{ over } Ax = b,\ x \ge 0\big\}.$$

This contains both sentences of the paper's THEOREM: the new basis yields a minimum for at least one $\lambda$ (namely $\bar\lambda$), and the left end $\underline\lambda'$ of its interval of optimality equals $\bar\lambda$. The statement uses "least element" and so does not presuppose that the set is an interval.

Milestones

1. (3): under nondegeneracy, a basic feasible solution is a minimum at $\lambda$ if and only if $\alpha_j + \lambda\beta_j \le 0$ for all $j$.
2. (4)–(5): if (3) is consistent, it holds exactly for $\underline\lambda \le \lambda \le \bar\lambda$.
3. The reduced-cost optimality criterion and its nondegenerate converse, and feasibility of the basis after a simplex pivot (both already proved on the platform).
4. [1]: $\bar\lambda$ satisfies the inequalities (3') $\alpha'_j + \lambda\beta'_j \le 0$ of the new basis.
5. (7): $\alpha'_r = -\alpha_s/y_{rs}$ and $\beta'_r = -\beta_s/y_{rs}$.
6. (7)–(8): no $\lambda < \bar\lambda$ satisfies (3').
7. Case A stop: in the nondegenerate Case A setup, if every $y_{is} \le 0$, there is no minimum for any $\lambda > \bar\lambda$.
8. Summary (4): under the section's standing assumptions, the set of $\lambda$ for which a minimum exists is closed and connected.

Significance

The THEOREM is what makes the parametric procedure correct. The optimality intervals of successive bases meet end to end: the basis entered at $\bar\lambda$ takes over exactly where the old one stops, with no gap and no overlap. Together with the stopping rule this gives the paper's summary: starting from one optimal basis, the procedure produces a finite list of bases and characteristic values covering every $\lambda$ for which the problem has a minimum, and that set is a closed interval. The same structure underlies the piecewise linear, concave optimal value as a function of $\lambda$, the parametric right-hand-side method by duality, and the shadow-vertex simplex method.

The result is classical and its proof is short, but as printed it contains a sign slip: display (8) on p. 41 states the inequality in the wrong direction, so the contradiction the authors claim does not follow from what they wrote. The formal development states the corrected argument. To our knowledge none of these parametric statements has a machine-checked proof; the fixed-cost facts they rest on (the reduced-cost optimality criterion, Bertsimas–Tsitsiklis Theorem 3.1, and the feasibility of the pivoted basis, Theorem 3.2) are already proved on the platform in Shuze Chen's LinearOptimization library and are reused here as reference items.

Difficulty

The natural first idea is to observe that $\bar\lambda$ satisfies (3') and conclude. This proves only that the new basis is optimal at $\bar\lambda$. The statement that nothing below $\bar\lambda$ works is a statement about optimality, not about the inequalities (3'): one must know that failure of (3') at $\lambda$ means that $x'$ is not optimal at $\lambda$. That converse of the optimality test is false for a degenerate basic solution, which can be optimal while some $z_j - c_j$ is positive. The nondegeneracy assumption is therefore essential, and it must be transported to the new basic solution. The other ingredient is the pivot update of $z_j - c_j$ when one column leaves and another enters, which the paper uses through (7) without writing it down.

Formalization scope

The problem data are A : Matrix (Fin m) (Fin n) ℝ, b : Fin m → ℝ and cost vectors d d' : Fin n → ℝ; the parameter $\lambda$ is the Lean variable t (λ is a Lean keyword). Indices are 0-based. Feasible set, optimality, bases, reduced costs, basic directions, the coordinates $y_{ij}$ (pivotColumn) and basic feasible solutions (IsSimplexState) are the published LinearOptimization definitions. The commitments and readings of the paper's phrases are:

• Sign. The platform's reduced cost is $\bar c_j = c_j - c_B'B^{-1}A_j = -(z_j - c_j)$; $\alpha_j$ and $\beta_j$ are defined as the negatives of the reduced costs for $d$ and $d'$, and (3), (3'), (5), (7) are stated in the paper's sign.
• "Yields a minimum" means that the point is feasible and its cost is at most that of every feasible point (IsLpOptimal), never a comparison of values.
• $\bar\lambda$ is finite in the THEOREM; it is the real number $-\alpha_s/\beta_s$ for an index $s$ with $\beta_s > 0$ attaining the minimum. As a function, (5) is valued in WithTop ℝ and WithBot ℝ; no real infimum of a possibly empty set is used.
• Nondegeneracy is stated for every basic feasible solution, including in the Case A stopping rule and Summary (4). The stopping rule also carries Case A consistency, and Summary (4) carries the section's available basic feasible solution. Linear independence of the rows of $A$ is assumed where the platform's criteria need it; it follows from the existence of a basis.
• "By the simplex method" is the reduced-cost optimality criterion; "the simplex prescription" is the ratio test; "no minimum" in the stopping rule is the unbounded outcome, optimal value $-\infty$, with the absence of a minimizer stated explicitly.
• The paper's interval $\delta \le \lambda \le \phi$ of interest plays no role in the THEOREM and is dropped.

The goal is not satisfied by a degenerate reading: its set is defined by optimality, not by the inequalities (3'), and its hypotheses hold on an explicit instance, so the statement is neither the definition of the critical value restated nor vacuous.

Not formalized: Case B (the search for a first $\lambda$ with a minimum), Summary (1)–(2) on the whole procedure, the no-repetition remark, the computing experience (§3) and the worked example (§4). Contributions welcome include the pivot-update formula for reduced costs, which is reusable for any simplex-based argument, and the closedness of the set of cost vectors with a finite optimum.

Selected references

• S. Gass and T. Saaty, The computational algorithm for the parametric objective function, Naval Research Logistics Quarterly 2 (1955), 39–45. https://doi.org/10.1002/nav.3800020106
• T. Saaty and S. Gass, The parametric objective function, Journal of the Operations Research Society of America 2 (1954), 316–319. https://doi.org/10.1287/opre.2.3.316
• G. B. Dantzig, Maximization of a linear function of variables subject to linear inequalities, in T. C. Koopmans (ed.), Activity Analysis of Production and Allocation, Wiley, 1951, Ch. XXI (Cowles Commission Monograph No. 13).
• D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Optimization, Athena Scientific, 1997, §3.1–3.2 and §5.2.

Motivation

A Markov decision problem with the average reward criterion asks for a policy that maximises the long-run reward per period of a controlled finite Markov chain. It is the standard model for repeated operational decisions with no natural horizon: inventory replenishment, machine maintenance, queue admission and routing. For finitely many states and actions, Howard (1960) and Blackwell (1962) showed that a pure stationary average optimal policy exists, and policy iteration finds one.

Linear programming formulations are older for the discounted problem (d'Epenoux, 1960) and for chains with a single recurrent class (Manne, 1960). Without a structural assumption, a stationary policy can split the state space into several recurrent classes, the multichain case, and the gain becomes a vector rather than a number. Denardo and Fox (1968) gave a linear programming approach for this case. Hordijk and Kallenberg (1979) showed that a single pair of dual linear programs suffices: an optimal policy can be read off an extreme optimal solution of the dual program. Kallenberg's tract (1983, Chapter 4) is the systematic account of this theory, and this mission formalizes its §4.2–4.4.

Setting

Let $E$ be a finite set of states. Each state $i$ has a nonempty finite set $A(i)$ of actions, a reward $r_{ia}\in\mathbb R$ for each action, and transition probabilities $p_{iaj}\ge 0$ with $\sum_j p_{iaj}=1$. A policy $R$ chooses, at each epoch $t=1,2,\dots$, a probability distribution on $A(X_t)$ that may depend on the whole history. A stationary policy $\pi^\infty$ uses fixed weights $\pi_{ia}$ at every epoch. A pure stationary policy $f^\infty$ always plays $f(i)\in A(i)$. Write $P(\pi)_{ij}=\sum_a p_{iaj}\pi_{ia}$, $r(\pi)_i=\sum_a r_{ia}\pi_{ia}$, and $P(f)$, $r(f)$ in the pure case.

The average reward of $R$ from state $i$ is

$$\phi_i(R)=\liminf_{T\to\infty}\frac1T\sum_{t=1}^T\sum_j\sum_a\mathbb P_R(X_t=j,\ Y_t=a\mid X_1=i)\,r_{ja},$$

and $\hat\phi_i(R)$ is the same expression with $\limsup$. The AMD-value-vector is $\phi_i=\sup_R\phi_i(R)$, where $R$ ranges over all policies, and $R^*$ is average optimal if $\phi(R^*)=\phi$. A policy is Blackwell optimal if it maximises the discounted reward $v^\alpha(R)=\sum_{t\ge1}\alpha^{t-1}\mathbb E_R[r_{X_tY_t}]$ for all discount factors $\alpha$ in some interval $[\alpha_0,1)$.

For a stochastic matrix $P$, the stationary matrix $P^*$ is the Cesàro limit of $P^n$, and the deviation matrix is $D=(I-P+P^*)^{-1}-P^*$. For a pure policy, $\phi(f^\infty)=P^*(f)r(f)$ and $u(f^\infty)=D(f)r(f)$.

A vector $\tilde\phi$ is AMD-superharmonic if some $\tilde u$ satisfies $\tilde\phi_i\ge\sum_jp_{iaj}\tilde\phi_j$ and $\tilde\phi_i+\tilde u_i\ge r_{ia}+\sum_jp_{iaj}\tilde u_j$ for all $a\in A(i)$, $i\in E$. For weights $\beta_j>0$ with $\sum_j\beta_j=1$, the primal program (4.2.10) minimises $\sum_j\beta_j\tilde\phi_j$ over such pairs. Its dual is

$$\max\Big\{\sum_{i,a}r_{ia}x_{ia}\ \Big|\ \sum_{i,a}(\delta_{ij}-p_{iaj})x_{ia}=0,\ \ \sum_ax_{ja}+\sum_{i,a}(\delta_{ij}-p_{iaj})y_{ia}=\beta_j\ \ (j\in E),\ \ x,y\ge0\Big\}.\tag{4.2.11}$$

Write $E_x=\{i\mid\sum_ax_{ia}>0\}$.

Formalization targets

Goal: Theorem 4.2.4

If $(x^*,y^*)$ is an optimal solution of (4.2.11) and an extreme point of its feasible set, then a decision rule $f_*$ with

$$x^*_{if_*(i)}>0\ \ (i\in E_{x^*}),\qquad y^*_{if_*(i)}>0\ \ (i\in E\setminus E_{x^*})$$

exists, and every such $f_*^\infty$ is average optimal.

The goal is stated for every finite model, with no assumption on the chain structure, and for every rule-conforming $f_*$.

Milestones

1. Optimality equations (Theorem 4.2.1): a Blackwell optimal pure policy gives $(\phi^\circ,u^\circ)$ with $\phi^\circ_i=\max_{a\in A(i)}\sum_jp_{iaj}\phi^\circ_j$ and $\phi^\circ_i+u^\circ_i=\max_{a\in\bar A(i)}\{r_{ia}+\sum_jp_{iaj}u^\circ_j\}$.
2. Superharmonic characterisation (Theorem 4.2.2): $\phi$ is the smallest AMD-superharmonic vector.
3. Lim sup optimality (Theorem 4.2.3): a pure stationary average optimal $f^\infty$ satisfies $\hat\phi(f^\infty)\ge\hat\phi(R)$ for every $R$.
4. The three propositions inside the proof of Theorem 4.2.4: complementary slackness along $f_*$ (4.2.1), closedness of $E_{x^*}$ under $P(f_*)$ (4.2.2), and transience of $E\setminus E_{x^*}$ under $P(f_*)$ (4.2.3).
5. The policy–solution correspondence (§4.3): the representative $(x(\pi),y(\pi))$ is feasible (Theorem 4.3.1); the map is one-to-one with inverse (4.3.1) (Theorem 4.3.2); $E_x$ is closed under, and is exactly the recurrent set of, $P(\pi(x,y))$ (Propositions 4.3.2–4.3.3); the correspondence preserves optimality (Theorem 4.3.3); pure policies map to extreme points (Theorem 4.3.4).
6. Policy evaluation (Theorem 4.4.1): the system $(I-P(f))\tilde\phi=0$, $\tilde\phi+(I-P(f))\tilde u=r(f)$, $\tilde u+(I-P(f))\tilde z=0$ is solvable, and every solution has $\tilde\phi=\phi(f^\infty)$ and $\tilde u=u(f^\infty)$.

Significance

Theorem 4.2.4 reduces the multichain average reward problem to one linear program. Combined with Theorems 4.3.3 and 4.3.4, every pure stationary average optimal policy arises from some extreme optimal solution, so enumerating extreme optimal solutions enumerates all of them (Remark 4.2.5). The same program is the starting point for the constrained average reward problems of §4.7, where a policy must also satisfy linear constraints on its state–action frequencies. Theorem 4.4.1 is the evaluation step of multichain policy iteration.

The results are proved in the tract and, earlier, in Hordijk and Kallenberg (1979). None of them is formalized in a proof assistant. The platform has machine-checked average reward results only for unichain models, where the gain is a scalar. This mission adds the multichain theory, with general history-dependent randomized policies as competitors, and also yields reusable statements about stationary and deviation matrices of finite chains.

Difficulty

Under a unichain assumption, the dual variables $x$ form a stationary distribution and the $y$-variables are unnecessary, so the obvious argument reads a policy off $x$ alone. In the multichain case the $x$-part of an optimal solution can vanish on whole sets of states, and on those states the policy must be read off $y$. The argument then has to show that these states are transient under the selected policy, which is where extremality is essential. For a non-extreme optimal solution the book passes instead to the randomized policy (4.3.1), which needs the separate argument of Theorem 4.3.3. A second obstacle is that optimality is against all history-dependent randomized policies, so the comparison of $\phi(f_*^\infty)$ with $\phi$ must pass through the superharmonic characterisation (Theorem 4.2.2), whose proof uses the existence of a Blackwell optimal policy.

Formalization scope

The model is the published MarkovDecisionProcesses.StationaryMDP: finite state and action types, a nonempty Finset of admissible actions per state, and stochastic rows. General policies are the published AvgHRPolicy (decision rules depending on the epoch and the full history). The criteria are the published gainInf, gainSup and optGainInf; all values are real. Average optimality is the book's lim inf notion, $\phi_i(R^*)=\sup_R\phi_i(R)$ for every $i$. It is not Puterman's stronger IsAverageOptimal. The stationary and deviation matrices are the published Cesàro limitMatrix and deviationMatrix. Their existence and invertibility (Theorem 2.4.1) are never assumed as hypotheses. LP variables are indexed by admissible state–action pairs, so Set.extremePoints is the book's notion of extreme feasible solution. Recurrence is the classical finite-chain notion: every state accessible from $i$ leads back to $i$. The ergodic set of a recurrent state is the set of states accessible from it.

A formalization that assumes a unichain or ergodic model trivializes the mission (it reduces Theorem 4.2.4 to known unichain results). No statement here makes such an assumption, and the weights $\beta_j$ are strictly positive and sum to one, as in (4.2.10).

The proofs need finite Markov chain theory: the Cesàro limit, $P^*P=PP^*=P^*$, $P^*D=0$, transience versus $P^*_{\cdot j}=0$, and positivity of the stationary distribution on ergodic sets. They also need LP duality with complementary slackness and the existence of a Blackwell optimal pure policy. The last result is posed separately on the platform (Blackwell 1962, Theorem 5). LP duality is available as proved platform theorems (LinearOptimization.*), whose sign conventions must be matched. Contributions on the Markov chain layer are reusable across the whole series.

Selected references

• L. C. M. Kallenberg, Linear Programming and Finite Markovian Control Problems, Mathematical Centre Tracts 148, Mathematisch Centrum, Amsterdam, 1983. https://ir.cwi.nl/pub/13008
• A. Hordijk and L. C. M. Kallenberg, Linear programming and Markov decision chains, Management Science 25(4):352–362, 1979. https://doi.org/10.1287/mnsc.25.4.352
• D. Blackwell, Discrete dynamic programming, Annals of Mathematical Statistics 33(2):719–726, 1962. https://doi.org/10.1214/aoms/1177704593
• E. V. Denardo and B. L. Fox, Multichain Markov renewal programs, SIAM Journal on Applied Mathematics 16(3):468–487, 1968.
• J. G. Kemeny and J. L. Snell, Finite Markov Chains, Van Nostrand, 1960.

Motivation

Fitting a probability distribution to prescribed margins or moments is a basic task in statistics and information theory. Starting from a reference distribution, one can enforce one constraint at a time by choosing the closest distribution in I-divergence, then cycle through the constraints. In a contingency table this procedure is iterative proportional fitting. A single projection onto the intersection is the desired answer, but each step of the iteration sees only one constraint set. Csiszár's 1975 paper proves that cyclic I-projections converge to the joint I-projection when the underlying set is finite and the constraint sets are linear. This mission formalizes that result together with the geometric statements needed to identify its limit.

Setting

Let $(X,\mathcal X)$ be a measurable space. A probability distribution (PD) is a probability measure on it. For PDs $P$ and $Q$, the I-divergence $I(P\|Q)$ is the Kullback–Leibler divergence: integrate the logarithm of the density of $P$ relative to $Q$ against $P$, with value $+\infty$ when $P$ is not absolutely continuous with respect to $Q$ or the integral diverges. It takes values in $[0,+\infty]$ and vanishes exactly when $P=Q$. The logarithmic conventions are $\log0=-\infty$ and $0\cdot(\pm\infty)=0$ (Csiszár 1975, (1.1)–(1.3)).

Given a PD $R$ and a set $\mathcal E$ of PDs that contains at least one $P$ with finite $I(P\|R)$, an I-projection $Q$ of $R$ onto $\mathcal E$ minimizes $I(P\|R)$ over $P\in\mathcal E$ and itself has finite divergence. The variation distance is the $L^1$ distance between the densities of two PDs relative to a common dominating measure. A set of PDs is linear here when every real affine combination of two members that is still a PD also belongs to the set. This includes distributions satisfying fixed linear moment constraints (Csiszár 1975, pp. 147 and 150–151).

For the goal, $X$ is finite. Choose $k\ge1$ linear sets $\mathcal E_1,\ldots,\mathcal E_k$, and write $\mathcal E=\bigcap_{i=1}^k\mathcal E_i$. Beginning with $Q_0=R$, let $Q_n$ be the I-projection of $Q_{n-1}$ onto $\mathcal E_i$ when $n=mk+i$ and $1\le i\le k$. The common feasible set is assumed nonempty and contains a PD absolutely continuous with respect to $R$. These assumptions ensure that the intended finite-divergence projection problem is present (Csiszár 1975, Theorem 3.2).

Formalization targets

The goal is the convergence statement of Theorem 3.2:

$$\forall x\in X,\qquad Q_n(\{x\})\longrightarrow Q(\{x\}), \quad Q=\operatorname{proj}^{I}_{\mathcal E}(R).$$

The goal also asserts existence of the cyclic sequence. On finite $X$, convergence of every point mass is equivalent to convergence in variation distance. The target fixes no rate. Its milestones include existence on variation-closed convex sets (Theorem 2.1), the identity (2.6) and the tangent-hyperplane Lemma 2.1, the Pythagorean characterization and transitivity of I-projections (Theorems 2.2–2.3), the density form under general linear constraints (Theorem 3.1) and the finite-moment exponential-density characterization (Corollary 3.1(A)), and the finite-space identity used in Theorem 3.2's proof. The two displayed identities in that proof are also targets:

$$I(Q\|R)=I(Q\|Q_n)+\sum_{i=1}^{n}I(Q_i\|Q_{i-1}),\qquad I(P\|Q_n)=I(P\|Q)+I(Q\|Q_n)\quad(P\in\mathcal E).$$

They specify the divergence relationships that the eventual convergence theorem must preserve (Csiszár 1975, (3.14)–(3.15)).

Significance

The theorem justifies enforcing finitely many linear distribution constraints one at a time while still obtaining the same I-projection as a simultaneous minimization. It gives a precise limit for the iterative proportional fitting procedure discussed by Csiszár, including cases in which the initial distribution has zero probabilities; strict positivity of $R$ is not required (Csiszár 1975, p. 155). The accompanying identities make the relation between a partial iterate and the final solution explicit. The result was proved in 1975. The work here is to provide machine-checked proofs of that known result and its supporting geometry, while keeping the source's extended-value and boundary conventions visible.

Difficulty

Successive projections minimize divergence over different sets. Minimization for one set by itself does not identify the minimizer over their intersection, and ordinary metric projection facts do not directly apply because I-divergence is asymmetric and can be infinite. The paper's statements need to handle zero densities and feasible distributions at infinite divergence without turning extended-real expressions into default values. On the finite space, convergence of a subsequence is available, but the formal claim must still identify every possible limit as the joint I-projection. These are the mathematical and representation difficulties faced by a complete proof (Csiszár 1975, §§2–3).

Formalization scope

Lean uses Measure X and Mathlib's InformationTheory.klDiv : ℝ≥0∞ for $I(P\|Q)$. A projection is a finite-divergence minimizer over the stated feasible set; no real-valued infimum or implicit arbitrary witness replaces it. Variation distance uses the densities relative to $P+Q$, which dominates both PDs and gives the value of (1.6). Variation closure is expressed by sequential convergence in this distance. The integral $\int\log q_R\,dP$ is represented in extended reals using positive and negative logarithmic parts and ENNReal.log; the theorem hypotheses keep the indeterminate $+\infty-(+\infty)$ case out. Density equations are $R$-almost everywhere, and a claim that $Q$ has an $R$-density explicitly includes $Q\ll R$.

For moment constraints, every integral is accompanied by integrability. Real affine combinations are tested through their values on measurable sets; the candidate combination must be a PD. The goal assumes a finite carrier with measurable singleton sets and $k\ge1$. Lean indexes the sets by Fin k from zero, so Qs (n+1) is projected onto the set at n % k; this is the paper's rule (3.13) shifted by one. The goal asserts existence of a recursive sequence before speaking about every such sequence, and the required feasible witness has $P\ll R$. These choices rule out an empty cyclic schedule, junk integrals, and a vacuous convergence assertion. Reusable contributions include the variation-distance and I-projection interfaces, their Pythagorean geometry, and the finite-moment density characterization; the mission's main closure remains Theorem 3.2.

Selected references

• I. Csiszár, I-Divergence Geometry of Probability Distributions and Minimization Problems, The Annals of Probability 3(1), 146–158 (1975). DOI: 10.1214/aop/1176996454.

Motivation

Delivering goods from one depot to many customers with a fleet of trucks is the vehicle routing problem. Dantzig and Ramser formulated it in 1959 as the "truck dispatching problem" (Management Sci. 6 (1959), doi:10.1287/mnsc.6.1.80). Five years later Clarke and Wright proposed the savings procedure (Oper. Res. 12 (1964), 568–581, doi:10.1287/opre.12.4.568): start with one truck per customer and repeatedly join two routes where joining saves the most distance, as long as the trucks can still carry the loads. On the twelve-customer example of Dantzig and Ramser it reached 290 units against their 294.

The savings procedure became the standard construction heuristic for vehicle routing. It appears in every textbook treatment of the problem, is the initial solution of many metaheuristics, and is the benchmark against which later construction methods are compared. The paper states the procedure and its bookkeeping (the half matrix, the vector $Q$, Table II) but proves nothing about it. That every run of the procedure ends, and ends in a feasible allocation, is left to the reader.

Setting

A depot $P_0$ serves customers $P_1,\dots,P_M$. The distance $d_{y,z}$ between every two points is given, and customer $P_j$ requires a load $q_j$. Trucks come in $n$ capacity classes: $x_i$ trucks of capacity $C_i$ are available, with $C_1<\dots<C_n$, and the trucks of the smallest capacity are unlimited, $x_1=\infty$ (p. 569).

A run is the ordered list of customers $P_{a_1},\dots,P_{a_k}$ served by one truck, which drives $P_0P_{a_1}\cdots P_{a_k}P_0$. Its load is $\sum_i q_{a_i}$. A state of the procedure is a list of runs; its mileage is the total length of its runs. A list of runs is an allocation if every customer lies on exactly one run, and it is feasible if each run can be given a truck of capacity at least its load without using more trucks of any class than are available.

The saving of the cell $(y:z)$ is $d_{0,y}+d_{0,z}-d_{y,z}$. The half matrix records $t_{y,z}=1$ if $P_y$ and $P_z$ are adjacent on a run, and $t_{y,0}\in\{0,1,2\}$ counts how many ends of runs lie at $P_y$. Table II counts, for each capacity level $C_i$, the runs with load above $C_i$ and the trucks with capacity above $C_i$.

The procedure starts with the runs $[P_1],\dots,[P_M]$, so $t_{y,0}=2$ for every customer. A cell $(y:z)$ is admissible if (I) $t_{y,0}>0$ and $t_{z,0}>0$; (II) $P_y$ and $P_z$ are on different runs; (III) after replacing those two runs by one run of load $Q_y+Q_z$, no column of Table II has more runs than trucks. A step links an admissible cell of maximum saving, ties broken arbitrarily, by joining the two runs end to end at $P_y$ and $P_z$. The procedure stops when no cell is admissible.

Formalization targets

Goal: correctness of the procedure

Assume symmetric distances, $C_1<\dots<C_n$, $x_1=\infty$, and that the initial one-truck-per-customer allocation passes the Table II test. Then every run of the procedure is finite whatever the tie-breaks; at every reachable state from which a link is possible a step exists; and every final state $S$ is a feasible allocation satisfying relation (A), $\sum_{z\ne y}t_{y,z}=2$ for every customer, with

$$\text{mileage}(S)=2\sum_{j=1}^M d_{0,j}-\sum_{\substack{1\le y<z\le M\\ t_{y,z}=1}}\bigl(d_{0,y}+d_{0,z}-d_{y,z}\bigr).$$

Milestones

$$\text{mileage}(\text{link}(S,y,z))=\text{mileage}(S)-(d_{0,y}+d_{0,z}-d_{y,z})\quad\text{for admissible }(y:z),$$

relation (A) at every reachable state, the permanence of links between customers, and

$$\#\{\text{runs with load}>C_i\}\le\sum_{k>i}x_k\ \ (i=1,\dots,n)\iff\text{the runs can be allocated to trucks},$$

which together give feasibility of every reachable state. Two further milestones: when $C_n\ge\sum_j q_j$ and the distances are a metric, the optimum equals the traveling salesman optimum (p. 569); and on the data of Table I every run of the procedure ends with total distance 290.

Significance

The goal is the specification the paper's procedure meets: it always terminates, never gets stuck, and outputs routes the fleet can actually drive, with the mileage the half-matrix bookkeeping predicts. The equivalence of the Table II test with truck feasibility is the reason the procedure can check capacities by counting columns instead of solving an assignment problem; it is a nested-class instance of Hall's marriage condition and is reusable for any routing or bin-assignment problem with ordered vehicle classes.

None of these statements is formalized anywhere, and the paper does not prove them. The paper makes no claim of optimality ("near-optimal", p. 568, is not quantified), and the mission does not either. A formal model of the procedure as a nondeterministic transition system is a base on which later results about savings heuristics (worst-case ratios, parallel and sequential variants) can be stated.

Difficulty

The obvious argument for termination, "each link removes a run", needs the invariant that every reachable state is a partition of the customers into runs; the procedure manipulates ordered lists and reverses runs, so the invariant must be carried through every step. Feasibility is not preserved by an arbitrary link but only by one that passes condition (III), and turning the column test into an actual assignment of runs to trucks is a matching argument that fails without both $x_1=\infty$ and the ordering of the capacities. Because ties are broken arbitrarily, every statement must hold for all runs of a nondeterministic procedure, not for one canonical execution. The worked example requires following every tie-break.

Formalization scope

Points are Fin (M+1) with depot 0, and a run's length is SupplyChainTheory.routeCost from the referenced module SupplyChainTheory_vrp. Capacity class i : Fin (n+1) is the paper's $C_{i+1}$; availabilities are in ℕ∞; loads and capacities are real. A state is a list of runs, each an ordered list of customers; the matrix $t$ and the vector $Q$ are computed from the state. The procedure is a step relation, and every invariant is stated for reachable states. Termination is Acc of the step relation at the initial state.

Choices committed to:

• distances are arbitrary real, symmetric numbers (the half matrix, p. 573); no nonnegativity or triangle inequality is assumed in the goal, and savings may be negative;
• every run: ties are free, as the paper suggests choosing randomly;
• reachable states: invariants are not claimed for arbitrary lists of runs;
• Table II is read cumulatively: column "Over $C_i$" compares runs with load $>C_i$ against all trucks of capacity $>C_i$, as the printed Tables II, IV and VI show;
• $x_1=\infty$ and $C_1<\dots<C_n$ are hypotheses; $q_j\le C_n$ is not separately assumed, since it follows from the initial Table II test;
• the initial one-truck-per-customer allocation is assumed feasible (p. 572);
• the metric hypotheses enter only the traveling-salesman milestone.

A correctness claim only about final states, without termination and progress, would hold for a procedure with no moves; a step relation with a fixed tie-break would prove less than the paper; invariants for arbitrary states are false; and a mileage identity for an arbitrary partition says nothing about the procedure's output. None of these is the target.

Not formalized: the informal $C_1\ll\sum q_j$, the shadow-cost reading of p. 571, the decomposition savings (2)–(5) of the general scheme, the load-splitting reduction of pp. 572–573, the Dantzig–Ramser methods, and the empirical comparisons and appendices.

Contributions welcome: proofs of the milestones, in particular the Table II equivalence, which needs a Hall-type argument for nested classes, and a decision procedure for the worked example.

Selected references

• G. Clarke and J. W. Wright, Scheduling of vehicles from a central depot to a number of delivery points, Operations Research 12(4) (1964), 568–581. https://doi.org/10.1287/opre.12.4.568
• G. B. Dantzig and J. H. Ramser, The truck dispatching problem, Management Science 6(1) (1959), 80–91. https://doi.org/10.1287/mnsc.6.1.80
• P. Hall, On representatives of subsets, Journal of the London Mathematical Society 10 (1935), 26–30. https://doi.org/10.1112/jlms/s1-10.37.26

Why natural policy gradient

Policy gradient methods optimize a parameterized policy directly by gradient ascent on its expected discounted return. They underlie most practical reinforcement-learning systems, yet for a long time their convergence was understood only locally: the objective $\theta\mapsto V^{\pi_\theta}(\rho)$ is not concave, and plain gradient ascent can stall on plateaus. The natural policy gradient (NPG) of Kakade (2001) preconditions the gradient by the pseudoinverse of the Fisher information matrix of the policy, so that steps are measured in the geometry of the induced action distributions rather than in parameter space. NPG is the conceptual ancestor of trust-region methods such as TRPO and PPO.

Agarwal, Kakade, Lee and Mahajan (arXiv:1908.00261v5, JMLR 2021) prove that exact NPG with softmax policies converges to the optimal value at rate $O(1/T)$, with constants that depend on neither the number of states, nor the number of actions beyond $\log|\mathcal A|$, nor any distribution mismatch coefficient. This mission formalizes that result, Theorem 5.3 of the paper, and the lemmas its proof rests on.

Timeline. Kakade (2001) introduced NPG and observed that, for softmax policies, it moves toward the greedy policy. Kakade and Langford (2002) proved the performance difference lemma. Even-Dar, Kakade and Mansour (2009) analyzed experts algorithms in changing MDPs, whose regret argument this proof adapts; Geist, Scherrer and Pietquin (2019) studied regularized policy iteration with related rates. Agarwal et al. (preprint 2019, journal 2021) gave the dimension-free $\log|\mathcal A|/(\eta T)+1/((1-\gamma)^2T)$ bound formalized here.

Setting

A finite discounted MDP consists of finite sets $\mathcal S$ (states) and $\mathcal A$ (actions, nonempty), a transition kernel $P(s'\mid s,a)$, rewards $r(s,a)\in[0,1]$ and a discount factor $\gamma\in[0,1)$. A policy $\pi$ assigns to each state a distribution $\pi(\cdot\mid s)$ over actions. Its value $V^\pi(s)$ is the expected discounted sum of rewards started from $s$; for a start distribution $\rho$, $V^\pi(\rho)=\sum_s\rho(s)V^\pi(s)$. The state–action value is $Q^\pi(s,a)=r(s,a)+\gamma\sum_{s'}P(s'\mid s,a)V^\pi(s')$ and the advantage is $A^\pi(s,a)=Q^\pi(s,a)-V^\pi(s)$. The discounted state visitation distribution is

$$d^\pi_\rho(s)=(1-\gamma)\sum_{t\ge0}\gamma^t\,\Pr{}^\pi(s_t=s\mid s_0\sim\rho).$$

An optimal policy $\pi^\star$ has $V^{\pi^\star}(s)\ge V^\pi(s)$ for every policy $\pi$ and every state $s$; write $V^\star=V^{\pi^\star}$.

The softmax policy with parameter $\theta\in\mathbb R^{|\mathcal S||\mathcal A|}$ is $\pi_\theta(a\mid s)=\exp(\theta_{s,a})/\sum_{a'}\exp(\theta_{s,a'})$. The Fisher information matrix with respect to $\rho$ is

$$F_\rho(\theta)=\mathbb E_{s\sim d^{\pi_\theta}_\rho}\mathbb E_{a\sim\pi_\theta(\cdot\mid s)}\big[\nabla_\theta\log\pi_\theta(a\mid s)\,\nabla_\theta\log\pi_\theta(a\mid s)^\top\big],$$

and the NPG update with step size $\eta$ is

$$\theta^{(t+1)}=\theta^{(t)}+\eta\,F_\rho(\theta^{(t)})^\dagger\,\nabla_\theta V^{(t)}(\rho),\tag{15}$$

where $M^\dagger$ is the Moore–Penrose pseudoinverse and $\pi^{(t)},V^{(t)},A^{(t)}$ are the quantities of $\theta^{(t)}$. Following §5.3, every state is assumed reachable from $\rho$: for each $s$ some policy $\pi$ has $d^\pi_\rho(s)>0$. The normalizer $Z_t(s)=\sum_a\pi^{(t)}(a\mid s)\exp(\eta A^{(t)}(s,a)/(1-\gamma))$ appears in the lemmas.

Formalization targets

Goal: Theorem 5.3

For NPG (15) from $\theta^{(0)}=0$ with any fixed $\eta>0$, for all $T>0$,

$$V^{(T)}(\rho)\ \ge\ V^\star(\rho)-\frac{\log|\mathcal A|}{\eta T}-\frac{1}{(1-\gamma)^2T}.$$

Companion: iteration complexity

If moreover $\eta\ge(1-\gamma)^2\log|\mathcal A|$ and $\epsilon>0$, then $V^\star(\rho)-V^{(T)}(\rho)\le\epsilon$ for every $T\ge 2/((1-\gamma)^2\epsilon)$.

Milestones

1. Lemma 3.2 (performance difference): $V^\pi(s_0)-V^{\pi'}(s_0)=\frac{1}{1-\gamma}\mathbb E_{s\sim d^\pi_{s_0}}\mathbb E_{a\sim\pi(\cdot\mid s)}[A^{\pi'}(s,a)]$.
2. Lemma 5.1 (soft policy iteration, policy form): $\pi^{(t+1)}(a\mid s)=\pi^{(t)}(a\mid s)\exp(\eta A^{(t)}(s,a)/(1-\gamma))/Z_t(s)$.
3. Lemma 5.2 (improvement): $V^{(t+1)}(\mu)-V^{(t)}(\mu)\ge\frac{1-\gamma}{\eta}\mathbb E_{s\sim\mu}\log Z_t(s)\ge0$ for every $\mu\in\Delta(\mathcal S)$.
4. Monotonicity (proof of Theorem 5.3): $V^{(t+1)}(s)\ge V^{(t)}(s)$ for all $s$.

Significance

The result. Theorem 5.3 is a global, non-asymptotic guarantee for a non-concave policy optimization problem, with no dependence on $|\mathcal S|$ and no distribution mismatch coefficient: the same analysis does not go through for plain softmax policy gradient, whose rates in the same paper carry such coefficients or are only asymptotic. It explains why Fisher preconditioning helps, and it is the template for the paper's later analysis of NPG with function approximation (§6) and for subsequent work on linear convergence of NPG and on entropy-regularized variants.

Formalizing it. The theorem is proved on paper; no machine-checked proof is known. A formal development needs: a policy-gradient identity for softmax policies, the computation of the Moore–Penrose pseudoinverse of a singular, block-structured Fisher matrix acting on the policy gradient, the performance difference lemma in the unnormalized convention, and a telescoping regret argument with Kullback–Leibler divergences between action distributions. Each of these is reusable: the performance difference lemma and the visitation-distribution calculus serve every result of the paper and much of the reinforcement-learning literature.

Difficulty

The obvious approach, treating NPG as gradient ascent on a smooth objective, fails: $V^{\pi_\theta}(\rho)$ is non-concave and its smoothness constant grows with $1/(1-\gamma)^3$, which would give a rate depending on the state space and on distribution mismatch. The dimension-free rate comes only after the NPG step is rewritten exactly as a multiplicative-weights update on each state's action distribution, and this rewriting is itself delicate. The Fisher matrix is always singular (each state's block annihilates the all-ones vector), the pseudoinverse returns the minimum-norm solution, which is not $A^{(t)}/(1-\gamma)$ but differs from it by a state-dependent offset; and the identity holds at a state only if that state carries positive visitation weight, which is why reachability from $\rho$ is assumed. A newcomer's first idea, replacing $F^\dagger$ by "any solution of $Fu=\nabla V$" or by an ordinary inverse, does not define the algorithm.

Formalization scope

Lean conventions (Mathlib, namespace PolicyGradTheory.NPG):

• Policies are functions S → A → ℝ (π s a $=\pi(a\mid s)$) satisfying the published IsPolicy; transitions P s a s' $=P(s'\mid s,a)$; values, $Q$-functions and state occupation probabilities are the published FoundationsML definitions PolicyValue, QFunction, OccupationDist, and optimality is IsOptimalPolicy. $V^\star(\rho)$ is the value of a policy assumed optimal, never a supremum over arbitrary functions.
• Parameters live in EuclideanSpace ℝ (S × A), so norms are Euclidean and gradient is $\nabla_\theta$.
• The pseudoinverse is the published predicate VectorSpaceOpt.IsPseudoinverse (minimum-norm least-squares inverse, which on Euclidean space is the Moore–Penrose pseudoinverse); Mathlib has none.
• Standing assumptions carried by every statement: rewards in $[0,1]$, $\gamma\in[0,1)$, $\mathcal A$ nonempty, $\rho\in\Delta(\mathcal S)$, every state reachable from $\rho$ (§5.3, footnote 6). The goal also assumes $\theta^{(0)}=0$ and $\eta>0$, as the page does. Lemmas 5.1 and 5.2 hold from any start.
• $\log$ is the natural logarithm; $T$ counts updates and $T>0$ is a natural number.

The goal does not assume the soft-policy-iteration form, the improvement bound or monotonicity; those are milestones to be proved, and a statement taking any of them as a hypothesis would not be Theorem 5.3. The parameter form "$\theta^{(t+1)}=\theta^{(t)}+\frac{\eta}{1-\gamma}A^{(t)}$" of Lemma 5.1 is not stated, since the pseudoinverse step differs from it by a state-dependent offset (as the lemma's own proof on p. 70 shows); only the policy form, which is invariant to that offset, is a target.

Contributions welcome: the softmax policy-gradient identity and score-function formula, existence and uniqueness of the minimum-norm least-squares inverse on Euclidean space, Kullback–Leibler identities for finite distributions, and the summability facts behind the visitation distribution.

Selected references

• A. Agarwal, S. M. Kakade, J. D. Lee, G. Mahajan, On the Theory of Policy Gradient Methods: Optimality, Approximation, and Distribution Shift, JMLR 22(98), 2021; preprint arXiv:1908.00261v5, 2020. https://arxiv.org/abs/1908.00261v5
• S. M. Kakade, A Natural Policy Gradient, NeurIPS 2001. https://proceedings.neurips.cc/paper/2001/hash/4b86abe48d358ecf194c56c69108433e-Abstract.html
• S. Kakade, J. Langford, Approximately Optimal Approximate Reinforcement Learning, ICML 2002. https://dl.acm.org/doi/10.5555/645531.656005
• E. Even-Dar, S. M. Kakade, Y. Mansour, Online Markov Decision Processes, Mathematics of Operations Research 34(3), 2009. https://doi.org/10.1287/moor.1090.0396
• M. Geist, B. Scherrer, O. Pietquin, A Theory of Regularized Markov Decision Processes, ICML 2019. https://arxiv.org/abs/1901.11275