All missions

Motivation

Primal-dual interior-point methods are the standard polynomial-time algorithms for linear programming and the basis of most commercial LP solvers. They do not move along the boundary of the feasible polyhedron, as the simplex method does. Instead they follow a curve through its interior, the central path, toward the solution set. Every convergence and complexity result about these methods assumes that this curve exists. Stephen J. Wright's monograph Primal-Dual Interior-Point Methods (SIAM, 1997, doi:10.1137/1.9781611971453) proves existence in its Chapter 2 from the standard LP background: existence of solutions, boundedness of solution sets, the Goldman–Tucker partition, and a primal-dual logarithmic barrier.

This mission is the first in a series that formalizes the book chapter by chapter. Later missions (potential reduction, path following, infeasible methods, finite termination, the homogeneous self-dual form) take the objects defined here as given.

Setting

Let $A \in \mathbb{R}^{m\times n}$, $b \in \mathbb{R}^m$ and $c \in \mathbb{R}^n$. Nothing is assumed about the rank of $A$. The primal and dual linear programs in standard form are

$$\text{(2.1)}\quad \min_x\ c^{\mathsf T}x \ \ \text{s.t.}\ Ax = b,\ x \ge 0, \qquad\qquad \text{(2.2)}\quad \max_{\lambda,s}\ b^{\mathsf T}\lambda \ \ \text{s.t.}\ A^{\mathsf T}\lambda + s = c,\ s \ge 0,$$

with all vector inequalities componentwise. The primal-dual feasible set $\mathcal F$ consists of the triples $(x,\lambda,s)$ with $Ax = b$, $A^{\mathsf T}\lambda + s = c$ and $(x,s) \ge 0$. The strictly feasible set $\mathcal F^\circ$ is the same set with $(x,s) > 0$. The solution sets are $\Omega_P$, the minimizers of (2.1); $\Omega_D$, the maximizers $(\lambda^*, s^*)$ of (2.2); and $\Omega = \Omega_P \times \Omega_D$. The index sets $\mathcal B$ and $\mathcal N$ collect the coordinates $j$ for which some $x^* \in \Omega_P$ has $x^*_j \neq 0$, and some $(\lambda^*, s^*) \in \Omega_D$ has $s^*_j \ne 0$, respectively.

For $\tau > 0$, the central path conditions (2.25) ask for $(x,\lambda,s)$ with

$$A^{\mathsf T}\lambda + s = c, \qquad Ax = b, \qquad x_j s_j = \tau \ \ (j = 1,\dots,n), \qquad (x,s) > 0 .$$

The reduced strictly feasible set $\mathcal H^\circ$ is the set of pairs $(x,s)$ such that $(x,\lambda,s) \in \mathcal F^\circ$ for some $\lambda \in \mathbb{R}^m$. The primal-dual log barrier on it is

$$f_\tau(x,s) = \frac1\tau\, x^{\mathsf T}s - \sum_{j=1}^n \log(x_j s_j).$$

In Lean all of these objects live in the single definition module WrightPDIPM.CentralPath.LP, under the names F, Fo, OmegaP, OmegaD, Omega, indexB, indexN, Hcirc, ftau and CentralPathConditions.

Formalization targets

Goal: Theorem 2.8 (p. 39)

If $\mathcal F^\circ \ne \emptyset$ and $\tau > 0$, then two things hold. First, $f_\tau$ has exactly one local minimizer relative to $\mathcal H^\circ$:

$$\exists!\,(x,s) \in \mathcal H^\circ \ \text{such that}\ (x,s) \text{ is a local minimizer of } f_\tau \text{ on } \mathcal H^\circ .$$

Second, the central path conditions (2.25) have a solution $(x,\lambda,s)$.

Milestones

1. Theorem 2.1 (p. 25). If $\mathcal F \neq \emptyset$, then $\Omega \ne \emptyset$. If either problem has an optimal solution, so does the other, and the optimal values agree.
2. Theorem 2.3 (p. 26). Suppose $\mathcal F \neq \emptyset$. If the dual is strictly feasible, then $\Omega_P$ is nonempty and bounded. If the primal is strictly feasible, then $\{s^* : (\lambda^*, s^*) \in \Omega_D\}$ is nonempty and bounded.
3. Theorem 2.4 (Goldman–Tucker, p. 28). $\mathcal B \cup \mathcal N = \{1,\dots,n\}$; equivalently, some $x^* \in \Omega_P$ and $(\lambda^*, s^*) \in \Omega_D$ satisfy $x^* + s^* > 0$. It is stated here under $\mathcal F \ne \emptyset$.
4. Lemma 2.5 (p. 30). If $\mathcal F^\circ \ne \emptyset$, then for every $K \ge 0$ the set $\{(x,s) : (x,\lambda,s) \in \mathcal F \text{ for some } \lambda,\ x^{\mathsf T}s \le K\}$ is bounded.

Significance

The result itself. Theorem 2.8 shows that the central path $\{(x_\tau, \lambda_\tau, s_\tau) : \tau > 0\}$ exists as soon as the problem has a strictly feasible primal-dual point. Every path-following and potential-reduction method in the book either tracks this curve or measures its iterates against neighbourhoods of it, so the theorem is the foundation the book's convergence and complexity results stand on. Theorem 2.4 provides the strictly complementary partition $\mathcal B \cup \mathcal N$, which the local and superlinear convergence analysis and the finite-termination results use. Theorems 2.1 and 2.3 are the duality and boundedness facts beneath both.

Formalizing it. These are classical results with complete textbook proofs. The work left is machine-checking them in the book's own form: standard-form LP, no rank assumption, the primal-dual barrier rather than the primal barrier, and minimality relative to $\mathcal H^\circ$. Related statements are already on the platform: Farkas's lemma and LP strong duality in Bertsimas–Tsitsiklis form, a central-path theorem in Matoušek–Gärtner form under the assumption $\operatorname{rank} A = m$ with the primal barrier, and Goldman–Tucker for the inequality-form LP. None of them is any of the statements above, but their proofs may be reusable.

Difficulty

The natural first step is to minimize $f_\tau$ by compactness. This does not work directly, because $\mathcal H^\circ$ is neither closed nor bounded and $f_\tau$ is not coercive on all of $\mathbb{R}^{2n}$. The quadratic term $x^{\mathsf T}s$ is indefinite in general, and only its restriction to the affine set of feasible pairs is linear. To obtain a minimizer one must show that the level sets of $f_\tau$ inside $\mathcal H^\circ$ are compact. This needs a boundedness statement (Lemma 2.5) that uses a strictly feasible point in an essential way, together with lower bounds that keep $x_j s_j$ away from $0$. A second difficulty is that $\mathcal H^\circ$ has empty interior in $\mathbb{R}^{2n}$ whenever $A \ne 0$. Passing from the minimizer to the conditions (2.25) is therefore a constrained first-order argument with multipliers, and it must work without $A$ having full row rank, so the multiplier $\lambda$ is not unique. For Goldman–Tucker, the obstacle is the step from "every index is positive at some optimum" to a single optimum that is positive on every index of $\mathcal B$ and of $\mathcal N$ at once.

Formalization scope

Vectors are Fin n → ℝ and matrices are Matrix (Fin m) (Fin n) ℝ. Inequalities are componentwise. Indices are 0-based in Lean and $\{1,\dots,n\}$ is Set.univ. Primal-dual points are triples (x, λ, s). "Bounded" means Bornology.IsBounded. $\Omega_P$ and $\Omega_D$ are defined as sets of optimal points (feasible, and at least as good as every feasible point), so no infimum or supremum is taken and an unattained optimum gives the empty set. No rank assumption on $A$ is made anywhere. No result in this chapter has an $O(\cdot)$ or an unspecified constant, so no constants are instantiated.

"Local minimizer in $\mathcal H^\circ$" means IsLocalMinOn (ftau τ) (Hcirc A b c), that is, minimality among nearby points of $\mathcal H^\circ$. Local minimality in all of $\mathbb{R}^{2n}$ would make the claim false. Uniqueness is asserted for the minimizer $(x,s)$ only, never for $\lambda$. $f_\tau$ uses Real.log and is only ever evaluated on $\mathcal H^\circ$, where every $x_j s_j > 0$.

The hypothesis $\mathcal F^\circ \neq \emptyset$ is part of the goal and cannot be weakened. If $\mathcal F^\circ = \emptyset$, then $\mathcal H^\circ = \emptyset$ and the first conclusion becomes false. A formalization that drops it, or one that replaces "unique local minimizer on $\mathcal H^\circ$" with "a global minimizer exists", is not this theorem.

The book prints Theorem 2.4 without a hypothesis. The faithful item carries the section's standing case $\mathcal F \ne \emptyset$, which by Theorem 2.1 gives $\Omega \neq \emptyset$.

Beyond the paper. One sharpness probe, theorem_2_4_probe, states Theorem 2.4 exactly as printed, with no hypothesis. It is expected to be false. For $m = n = 1$, $A = (1)$, $b = -1$, $c = 0$ the primal is infeasible and the dual objective is unbounded, so $\mathcal B \cup \mathcal N = \emptyset$. A disproof is the sought outcome.

Infrastructure that is reusable beyond this mission includes the standard-form LP objects, the solution sets and the $\mathcal B$/$\mathcal N$ partition, all of which later chapters use; KKT conditions for linearly constrained convex minimization; and strict convexity of $\sum_j -\log t_j$. Contributions of general lemmas toward any of these are welcome.

Selected references

• S. J. Wright, Primal-Dual Interior-Point Methods, SIAM, Philadelphia, 1997, Chapter 2. doi:10.1137/1.9781611971453
• A. J. Goldman and A. W. Tucker, "Theory of linear programming," in Linear Inequalities and Related Systems, Annals of Mathematics Studies 38, Princeton University Press, 1956, pp. 53–97.
• N. Megiddo, "Pathways to the optimal set in linear programming," in Progress in Mathematical Programming, Springer, 1989, pp. 131–158.
• D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Optimization, Athena Scientific, 1997.

Motivation

Statistical learning theory asks when minimizing the error on a finite sample also controls the error on unseen data. For a set of classifiers, the answer given by Vapnik and Chervonenkis is governed by a single combinatorial quantity: how many different ways the set can label $\ell$ points. Bounds on the generalization error of empirical risk minimization are stated in terms of the logarithm of this number, the growth function, and they are informative only if the growth function grows sublinearly in $\ell$.

The theorem at the centre of this mission says that there is no middle ground: the growth function is either exactly $\ell \ln 2$ for every $\ell$, or, from some point on, it is bounded by a logarithmic function $h(\ln(\ell/h) + 1)$. The integer $h$ is the VC dimension. This dichotomy is what turns the abstract bounds of Chapter 3 of V. N. Vapnik, The Nature of Statistical Learning Theory (2nd ed., Springer 2000, doi:10.1007/978-1-4757-3264-1) into constructive ones (Eq. (3.24) of the book).

Timeline. Vapnik and Chervonenkis announced the connection between the growth function and the shattering number in 1968 and proved it in 1971 (Theory Probab. Appl. 16 (1971)), in the form $m^S(r) = 2^r$ for all $r$ or $m^S(r) \le r^n + 1$. Sauer (J. Combin. Theory A 13 (1972)) and Shelah (Pacific J. Math. 41 (1972)) independently proved the sharp polynomial bound $\sum_{i=0}^{h}\binom{\ell}{i}$. The logarithmic form used in the book follows from that sum together with the estimate $\sum_{i \le h}\binom{\ell}{i} \le (e\ell/h)^h$. The book states its results without proofs (Preface, p. xi).

Setting

Let $Z$ be a nonempty set and $\Lambda$ a nonempty set of parameters, and let $Q(z, \alpha)$, $\alpha \in \Lambda$, be a set of indicator functions on $Z$: each takes only the values $0$ and $1$. For a sample $z_1, \dots, z_\ell \in Z$ (repetitions allowed), let

$$N^\Lambda(z_1, \dots, z_\ell)$$

be the number of different binary vectors $(Q(z_1, \alpha), \dots, Q(z_\ell, \alpha))$ obtained as $\alpha$ ranges over $\Lambda$ (Section 2.3.2, p. 42). The growth function is

$$G^\Lambda(\ell) = \ln \sup_{z_1, \dots, z_\ell} N^\Lambda(z_1, \dots, z_\ell)$$

(Section 2.7, p. 55), with $\ln$ the natural logarithm. Since $N^\Lambda \le 2^\ell$, always $G^\Lambda(\ell) \le \ell \ln 2$.

Points $z_1, \dots, z_h$ are shattered if for every $\delta_1, \dots, \delta_h \in \{0, 1\}$ some $\alpha$ has $Q(z_i, \alpha) = \delta_i$ for all $i$. The VC dimension is the maximum number of shattered points, or $\infty$ if arbitrarily many points are shattered (Section 3.6, p. 80). For a set of real functions $A \le Q(z, \alpha) \le B$, the VC dimension is that of the indicators $\theta\{Q(z, \alpha) - \beta\}$, $\beta \in (A, B)$ (Eq. (3.22)), where the step function is $\theta(u) = 0$ for $u < 0$ and $\theta(u) = 1$ for $u \ge 0$.

In the Lean development these are numSeparations, growthFunction, Shatters, vcDim, stepFunction and vcDimReal in the namespace NatureSLT.GrowthFunction.

Formalization targets

Goal: Theorem 3.3 (p. 79)

Either

$$G^\Lambda(\ell) = \ell \ln 2 \quad \text{for every } \ell,$$

or there is an integer $h \ge 0$ with $G^\Lambda(h) = h \ln 2$ and $G^\Lambda(h+1) < (h+1)\ln 2$ such that

$$G^\Lambda(\ell) \le h\left(\ln\frac{\ell}{h} + 1\right) \qquad \text{for every } \ell > h .$$

The range $\ell > h$ is the one the book gives in both of its restatements of the theorem: the display after the Definition on p. 79 and Eq. (3.23) on p. 83.

Milestones

1. Equivalent definition of the VC dimension (Section 3.6, p. 80): the VC dimension is infinite if and only if $G^\Lambda(\ell) = \ell\ln 2$ for all $\ell$, and equals $h$ if and only if $G^\Lambda(h) = h\ln 2$ and $G^\Lambda(h+1) < (h+1)\ln 2$.
2. Eq. (3.23) (p. 83): if the VC dimension is a finite $h$, then $G^\Lambda(\ell) \le h(\ln(\ell/h) + 1)$ for $\ell > h$.
3. Example 1(i) (p. 81): the linear indicator functions $\theta\{\sum_{p=1}^n \alpha_p z_p + \alpha_0\}$ on $\mathbb R^n$ have VC dimension $n + 1$.
4. Example 1(ii) (p. 81): the real linear functions $\sum_{p=1}^n \alpha_p z_p + \alpha_0$ have VC dimension $n + 1$.
5. Example 2(i) (p. 82): $\theta(\sin \alpha z)$, $\alpha \in \mathbb R$, has infinite VC dimension; the points $10^{-1}, \dots, 10^{-\ell}$ are shattered by the explicit parameter $\alpha = \pi(\sum_{i=1}^\ell (1 - \delta_i)10^i + 1)$.

Significance

The dichotomy is what makes the VC dimension a usable complexity measure. Combined with the bounds of Section 3.3, it gives distribution-free generalization bounds whose confidence term depends on $\ell$ and $h$ only (Eq. (3.24)). It also explains why finiteness of the VC dimension is a sufficient condition for consistency of empirical risk minimization for every probability measure (p. 80): the growth function, divided by $\ell$, tends to zero exactly when it is not linear. Examples 1 and 2 show that the VC dimension is not the number of parameters: linear functions have as many as their parameters, while a one-parameter family of sines has infinitely many.

The theorem is classical and has been proved many times. Related statements are already machine-checked on the platform in other forms: Sauer's lemma for concept classes (ComputationalLearning.sauer_lemma, FoundationsML.RademacherVC.sauer_lemma, HighDimProb.Chaining.sauer_shelah), the bound $\Phi_d(m) \le (em/d)^d$ (ComputationalLearning.phi_polynomial_bound), the 1971 dichotomy with the bound $r^n + 1$ for sets of events (VapnikChervonenkis.GrowthFunction.theorem1_growth_dichotomy), and the VC dimension of open halfspaces with $\pm$ labels (UnderstandingML.vcDim_halfspaces). What is not on the platform is the book's statement: the logarithmic bound for the growth function of an indexed set of $\{0,1\}$-valued functions, with $h$ characterized through the growth function itself, and the book's examples with the closed half-space convention $\theta(0) = 1$. This mission produces that statement and connects it to the shattering definition.

Difficulty

The obvious argument proves only half. That $G^\Lambda(h) = h\ln 2$ and $G^\Lambda(h+1) < (h+1)\ln 2$ pin down the largest shattered set is elementary; the content lies in the bound for $\ell > h$: a class that shatters no $h + 1$ points must realize only polynomially many labelings of every $\ell$ points, with the explicit constant of the logarithmic form, and nothing in the definition of the growth function gives this directly. A tempting shortcut is to read the bound as valid for every $\ell$, as the theorem's display does in isolation: this is false whenever $h \ge 2$, since $G^\Lambda(1) = \ln 2$ exceeds $h(\ln(1/h) + 1)$. The range $\ell > h$ is essential. For the examples, the upper bound $n + 1$ in Example 1 is a Radon-type statement about $n + 2$ points in $\mathbb R^n$, and Example 2 requires controlling the sign of $\sin(\alpha\, 10^{-i})$ simultaneously for all $i$ with one parameter.

Formalization scope

Indicator functions are boolean-valued, true standing for $1$, and written Q α z for the book's $Q(z, \alpha)$; samples are maps Fin ℓ → Z, so indices run from $0$. $N^\Lambda$ is the cardinality of the finite set of realized vectors in Fin ℓ → Bool; the supremum over samples is a supremum of natural numbers bounded by $2^\ell$, and $G^\Lambda$ is Real.log of its cast. The VC dimension is valued in ℕ∞. $Z$ and $\Lambda$ are assumed nonempty in every theorem, so the supremum is over a nonempty set of samples and is at least $1$; the junk value $\ln 0 = 0$ never enters. At $h = 0$ the book's expression $h(\ln(\ell/h)+1)$ is undefined; the formal statements evaluate it as $0$, the limit as $h \to 0^+$, and the correct bound in that case, where $G^\Lambda \equiv 0$. The step function has $\theta(0) = 1$. Points of $\mathbb R^n$ are EuclideanSpace ℝ (Fin n); in Example 1(ii) the bounds of the real-function definition are $A = -\infty$, $B = +\infty$.

A trivializing formalization is ruled out: the goal is the disjunction of a linear growth function and the explicit bound $h(\ln(\ell/h) + 1)$ with $h$ fixed by $G^\Lambda(h) = h \ln 2$ and $G^\Lambda(h+1) < (h+1)\ln 2$, not an existence statement over an arbitrary $h$ and not a bound in which an empty parameter set makes the growth function vanish.

Beyond the paper. One variant is posed beside the faithful goal: a sharpness probe (growth_function_dichotomy_probe) claiming the logarithmic bound for every $\ell$, i.e. without $\ell > h$. Its evidence: for $Z = \{1, \dots, 10\}$ and all subsets, $h = 10$ and $G^\Lambda(1) = \ln 2 > 10(\ln(1/10) + 1)$. A disproof is the sought outcome.

A complete development needs: finite cardinality arguments for sets of binary vectors, the Sauer–Shelah lemma for an indexed family (Mathlib has Finset.card_shatterer_le_sum_vcDim for finite set families), the estimate $\sum_{i\le h}\binom{\ell}{i} \le (e\ell/h)^h$, Radon's theorem or an affine-independence argument for Example 1, and elementary trigonometry for Example 2. The bridge between indexed families of indicator functions and set families is reusable by every later mission of this series. Contributions of any milestone, of bridge lemmas to the existing Sauer-type theorems, and of a disproof of the probe are welcome.

Selected references

• V. N. Vapnik, The Nature of Statistical Learning Theory, 2nd ed., Springer, 2000. https://doi.org/10.1007/978-1-4757-3264-1
• V. N. Vapnik and A. Ya. Chervonenkis, On the uniform convergence of relative frequencies of events to their probabilities, Theory of Probability and Its Applications 16(2), 1971. https://doi.org/10.1137/1116025
• N. Sauer, On the density of families of sets, Journal of Combinatorial Theory, Series A 13(1), 1972. https://doi.org/10.1016/0097-3165(72)90019-2
• S. Shelah, A combinatorial problem; stability and order for models and theories in infinitary languages, Pacific Journal of Mathematics 41(1), 1972. https://doi.org/10.2140/pjm.1972.41.247
• V. N. Vapnik, Statistical Learning Theory, Wiley, 1998. ISBN 978-0-471-03003-4

Motivation

A Jackson network of exponential single-server queues has a stationary distribution of product form: the numbers of customers at the nodes are independent in steady state, each geometric (Jackson, 1957; Jackson, 1963). The same property holds for a large family of networks (Kelly, Reversibility and Stochastic Networks, 1979), and it is what makes performance formulas for such networks explicit. A network can also be read as a quasi-birth-and-death process (QBD): take the number of customers at one node as the level and the vector of queue lengths at the other nodes as the phase. The phase set is then infinite, $\mathbb N^{K-1}$, and product form becomes the statement that level and phase are independent in steady state. Ramaswami and Taylor (Stochastic Models 12, 1996, pp. 121–142 and 143–164) characterized when a QBD has this independence property and used the characterization to find new product-form networks. Chapter 15 of Latouche and Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling (SIAM, 1999), gives a simplified version of these results, and this mission formalizes it.

Setting

Fix a countable set $\iota$ of phases, finite or infinite, and four real matrices $B, A_0, A_1, A_2$ indexed by $\iota$. A continuous-time QBD is the Markov process on the states $(n,i)$, with level $n\ge 0$ and phase $i\in\iota$, whose infinitesimal generator is block-tridiagonal:

$$Q=\begin{bmatrix} B & A_0 & & \\ A_2 & A_1 & A_0 & \\ & A_2 & A_1 & \ddots \\ & & \ddots & \ddots\end{bmatrix}.$$

$A_0$ and $A_2$ are nonnegative, $A_1$ and $B$ have nonnegative off-diagonal and strictly negative diagonal entries, and the rows of $Q$ sum to zero. The level $\ell(n)$ is the set $\{(n,i): i\in\iota\}$. The process stays in a state $x$ for an exponential time of rate $q_x=-Q_{xx}$ and then jumps to $y\ne x$ with probability $Q_{xy}/q_x$; this jump chain together with the mean sojourns $1/q_x$ describes it. It is assumed irreducible and regular (finitely many transitions in every bounded time interval). It is positive recurrent if it returns to every state with probability one and the expected return time is finite. A stationary probability vector $\pi=(\pi_0,\pi_1,\dots)$, partitioned by levels, satisfies $\pi\ge 0$, $\pi\mathbf 1=1$ and $\pi Q=0$.

Two matrices indexed by phases are defined by what they mean for the process:

• $N_{ij}$, the expected time spent in $(n,j)$, starting from $(n,i)$, before the first visit to $\ell(n-1)$ ($n\ge 1$);
• the rate matrix $R_{ij}$, the rate of sojourn in $(n+1,j)$ per unit of local time in $(n,i)$: $q_{(n,i)}$ times the expected time spent in $(n+1,j)$, starting from $(n,i)$, before the first return to $\ell(0)\cup\dots\cup\ell(n)$.

The level is independent of the phase in steady state if $\pi_n=\alpha_nv$ for some probability vector $v$ on $\iota$ and some probability density $\{\alpha_n\}$ on the levels. A vector $\boldsymbol\theta$ is positive if $\theta_i>0$ for every phase.

Formalization targets

Goal: the independence criterion (Theorem 15.1.3)

The process is positive recurrent with level and phase independent in steady state if and only if there are a positive vector $\boldsymbol\theta$ with $\boldsymbol\theta\mathbf 1=1$ and a scalar $0<\eta<1$ with

$$\boldsymbol\theta(A_0+\eta A_1+\eta^2A_2)=\mathbf 0,\qquad \boldsymbol\theta(B+\eta A_2)=\mathbf 0;$$

then $\pi_n=(1-\eta)\eta^n\boldsymbol\theta$ for all $n\ge 0$. The criterion uses only the blocks of $Q$, not $R$.

Milestones

1. Theorem 6.4.1. For a positive recurrent continuous-time QBD, $R$ and $N$ are finite and the same at every level, $R=A_0N$, and $\pi_n=\pi_0R^n$ for $n\ge 0$.
2. Theorem 15.1.1. The same equivalence as the goal, with $\boldsymbol\theta R=\eta\boldsymbol\theta$ in place of the first equation.
3. Theorem 15.2.1. For a Jackson network written as a QBD (level $n_1$, phase $(n_2,\dots,n_K)$), with $\boldsymbol y$ the solution of the traffic equations $\boldsymbol y=\boldsymbol\lambda+\boldsymbol yP$ and $y_i/\mu_i<1$ for all nodes, the scalar $\eta=y_1/\mu_1$ and the vector $\theta_{\boldsymbol r}=c\prod_{2\le j\le K}(y_j/\mu_j)^{n_j}$, $c=\prod_{2\le j\le K}(1-y_j/\mu_j)$, satisfy the conditions of the goal.

Significance

The goal says that independence of level and phase forces the level to be geometric, and it reduces the question to a linear system in $\boldsymbol\theta$ for a fixed $\eta$. That system involves only the generator's blocks, so it can be checked for models whose rate matrix has no closed form. Together with Theorem 15.2.1 it gives a QBD proof of the product form of open Jackson networks, and the same approach led Ramaswami and Taylor to product-form networks beyond Jackson's. The level-dependent version (§15.3) is described in the book but not proved there.

Formalizing the chapter adds the continuous-time layer of the matrix-analytic series: a generator on a countable state space, regularity, and positive recurrence defined by return times, with expected sojourn times as sums over trajectories of the jump chain. The product form of open migration processes, which contains the Jackson case, is on the platform as KellyStochasticNetworks.open_migration_equilibrium. It is a balance-equation statement about a network and does not concern QBDs. No QBD statement is on the platform.

Difficulty

The "if" direction of the goal looks like a computation: $(1-\eta)\eta^n\boldsymbol\theta$ satisfies $\pi Q=0$. But a solution of $\pi Q=0$ shows positive recurrence only through regularity, and on an infinite phase set that is a real hypothesis. Relating the generator condition $\boldsymbol\theta(A_0+\eta A_1+\eta^2A_2)=\mathbf 0$ to the eigenvector condition $\boldsymbol\theta R=\eta\boldsymbol\theta$ is the hard part. The book notes that the finite-dimensional argument (Perron–Frobenius theory for $R$) "fails in the case of an infinite number of phases": spectral arguments for $R$ are not available when $\iota$ is infinite, and every product of infinite matrices must be justified by nonnegativity. For Theorem 15.2.1, positivity of $\boldsymbol\theta$ and $\eta$ is not visible in the hypotheses on the load. It comes from irreducibility of the network.

Formalization scope

All declarations are in MatrixAnalytic.ProductForm. States are ℕ × ι with ι countable and arbitrary (levels from $0$; no finiteness of phases is assumed anywhere). The generator is qbdQ B A₀ A₁ A₂, and the hypothesis IsGeneratorC encodes the sign and zero-row-sum conditions. The continuous-time process is encoded through its jump chain and mean holding times, without measure theory. Irreducibility is that of the jump chain. Regularity (RegularCT) says that the sum of mean holding times along the jump chain diverges almost surely. Positive recurrence (PositiveRecurrentCT) says the return probability is one and the expected return time is finite. $R$ and $N$ are sums over taboo trajectories valued in $[0,+\infty]$, never solutions of $A_0+RA_1+R^2A_2=0$. Defining $R$ by that equation, or positive recurrence by the existence of a stationary vector, would make the milestones circular and is ruled out. Positive recurrence and stationarity are separate objects. "Level independent of phase" is the factorization $\pi_{(n,i)}=\alpha_nv_i$, not the geometric form, which is the conclusion. Vector equations over infinite phase sets are convergent sums; products of nonnegative matrices are computed in $[0,+\infty]$. The Jackson network indexes nodes by Option J with J finite, none being node 1. Its blocks are written out rate by rate, and the openness of the network ($I-P$ invertible) is carried as the meaning of "the solution" of the traffic equations.

A complete development needs the theory of continuous-time chains on countable spaces through the jump chain (invariant measures, positive recurrence, the regular case), the continuous-time matrix-geometric property, and the series representation of $R$. These are reusable well beyond this mission, and contributions of any of them are welcome.

Selected references

• G. Latouche and V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM Series on Statistics and Applied Probability, SIAM, 1999. doi:10.1137/1.9780898719734
• V. Ramaswami and P. G. Taylor, An operator-analytic approach to product-form networks, Communications in Statistics. Stochastic Models 12(1):121–142, 1996. (Reference [102] of Latouche & Ramaswami.)
• V. Ramaswami and P. G. Taylor, Some properties of the rate operators in level dependent quasi-birth-and-death processes with a countable number of phases, Communications in Statistics. Stochastic Models 12(1):143–164, 1996. (Reference [103] of Latouche & Ramaswami.)
• J. R. Jackson, Networks of waiting lines, Operations Research 5(4), 1957. doi:10.1287/opre.5.4.518
• J. R. Jackson, Jobshop-like queueing systems, Management Science 10(1), 1963. doi:10.1287/mnsc.10.1.131
• F. P. Kelly, Reversibility and Stochastic Networks, Wiley, 1979. statslab.cam.ac.uk
• J. R. Norris, Markov Chains, Cambridge University Press, 1997. doi:10.1017/CBO9780511810633

Motivation

A retailer of a seasonal or perishable product must commit to a stock level before demand is known, and must also choose the price at which to sell. The two decisions interact: the price shifts the demand distribution, and the stock level determines how much of that demand can be served. The price-setting newsvendor is the single-period model of this interaction. It goes back to Whitin (1955) and Mills (1959), and it is the building block of the joint pricing and inventory literature in operations research and economics.

Petruzzi and Dada, Pricing and the Newsvendor Problem: A Review with Extensions, Operations Research 47(2), 1999 (doi:10.1287/opre.47.2.183), unify the additive and multiplicative versions of the model and give conditions under which the joint optimum is characterized analytically. Their Theorem 1 treats additive demand. Earlier results required the demand noise to be a Pólya frequency function of order 2 (Ernst 1970) or extended this to the log-normal (Young 1978). Theorem 1 replaces those assumptions by a condition on the hazard rate that every distribution with a nondecreasing hazard rate satisfies.

Setting

Demand at price $p$ is $D(p,\varepsilon) = y(p) + \varepsilon$ with $y(p) = a - bp$, where $a > 0$, $b > 0$. The noise $\varepsilon$ is a random variable on the range $[A,B]$, $A < B$, with probability density $f$ (nonnegative, zero off $[A,B]$, integrating to $1$), distribution function $F(z) = \int_A^z f$ and mean $\mu$. The paper requires $A > -a$. Stocking costs $c$ per unit; each leftover unit is disposed of at cost $h$ with $h \ge -c$ (a negative $h$ is a salvage value); each unit of unmet demand costs a penalty $s$. No sign is imposed on $c$ or $s$.

The firm stocks $q = y(p) + z$ units, where $z$ is the stocking factor. Leftovers occur when $\varepsilon < z$, shortages when $\varepsilon > z$. The expected profit is

$$E[\Pi(z,p)] = \int_A^z \big(p[y(p)+u] - h[z-u]\big) f(u)\,du + \int_z^B \big(p[y(p)+z] - s[u-z]\big) f(u)\,du - c[y(p)+z].$$

With expected leftovers $\Lambda(z) = \int_A^z (z-u) f(u)\,du$ and expected shortages $\Theta(z) = \int_z^B (u-z) f(u)\,du$, it splits as $E[\Pi(z,p)] = \Psi(p) - L(z,p)$, a riskless profit $\Psi(p) = (p-c)[y(p)+\mu]$ less a loss function $L(z,p) = (c+h)\Lambda(z) + (p+s-c)\Theta(z)$ (eqs. (1)–(3)). The riskless price is $p^0 = (a + bc + \mu)/2b$. The firm solves

$$\max_{z \in [A,B],\ p \in \mathbb R} E[\Pi(z,p)]. \qquad (4)$$

For each fixed $z$ the best price is $p(z) = p^0 - \Theta(z)/2b$ (Lemma 1), so (4) reduces to maximizing $E[\Pi(z,p(z))]$ over $z$. The hazard rate is $r(z) = f(z)/[1-F(z)]$.

The Lean development lives in the namespace PetruzziDada.Additive. The definition Model carries the data $a,b,c,h,s,A,B,f$ with the standing assumptions above, and defines y, F, μ, Λ, Θ, expProfit (the integral above), Ψ, L, p0, price, reducedProfit ($z \mapsto E[\Pi(z,p(z))]$), R (its derivative), hazard and HazardCondition.

Formalization targets

Goal: Theorem 1(c)

If $2r(z)^2 + r'(z) > 0$ for $A \le z < B$ and $a - b(c-2s) + A > 0$, then there is $z^* \in [A,B]$ with

$$\Big\{ z \in [A,B] : \tfrac{d}{dz} E[\Pi(z,p(z))] = 0 \Big\} = \{z^*\}, \qquad E[\Pi(z,p)] \le E[\Pi(z^*,p(z^*))] \ \ \forall z \in [A,B],\ \forall p \in \mathbb R.$$

The optimal policy is then $q^* = y(p^*) + z^*$, $p^* = p(z^*)$. The goal leaves the density arbitrary beyond the stated condition.

Milestones, in the order of the paper's proof

1. Eqs. (1)–(3): $E[\Pi(z,p)] = \Psi(p) - L(z,p)$.
2. Eq. (5): $\partial E[\Pi(z,p)]/\partial z = -(c+h) + (p+s+h)[1-F(z)]$.
3. Lemma 1: $p(z)$ is the unique optimal price for each $z$.
4. Theorem 1(a): an exhaustive search over $z \in [A,B]$ determines the optimum.
5. Appendix, p. 192: $dE[\Pi(z,p(z))]/dz = -(c+h) + [p^0+s+h-\Theta(z)/2b][1-F(z)] =: R(z)$, and the formula for $dR/dz$.
6. Appendix, p. 193: at a critical point of $R$, $R'' = -\frac{f(1-F)}{2br^2}\{2r^2 + r'\}$.
7. Appendix, p. 193: under the hazard condition, $R$ is monotone or unimodal on $[A,B]$ and has at most two roots.
8. Appendix, p. 193: $R(B) = -(c+h)$ and $2bR(A) = a - b(c-2s) + A$.

Theorem 1(b) is also stated, as printed, but it is not a milestone (see Difficulty).

Significance

Theorem 1 tells a practitioner when the joint stocking and pricing problem needs only its first-order condition: one root-finding problem in $z$ replaces a two-dimensional search. The hazard condition $2r^2 + r' > 0$ holds for every density with a nondecreasing hazard rate, such as the uniform density on $[A,B]$, and it tolerates a hazard rate that decreases, as long as it does not decrease faster than $2r^2$. The result is the analytical core of much later work on joint pricing and inventory decisions.

The theorem is proved in the paper; it has not been machine-checked. Formalizing it puts the classical price-setting newsvendor on the platform: an objective defined by an integral whose limits move with the decision, a reduction by optimizing out the price, and a unimodality argument for a derivative. The definitions (the density model, $\Lambda$, $\Theta$, the hazard rate) are reusable for the multiplicative case of the same paper and for other newsvendor results.

Difficulty

The objective is an integral with variable limits and a density that is assumed only integrable, so eq. (5) is a Leibniz rule in which pointwise continuity of $f$ cannot be used. The derivative $R$ of the reduced objective is itself differentiable only where $f$ is, and its second derivative involves $f'$. The hazard condition, however, controls the sign of $R''$ only at critical points of $R$, and turning that local information into a global shape of $R$ on a closed interval, with a density that may jump at $A$, is a real-analysis lemma that Mathlib does not provide in this form.

Two published steps need care. The printed proof says "$R(B) = -(c+h) < 0$", but the model allows $c + h = 0$; Theorem 1(c) still holds then, by a different argument (finding F3). Theorem 1(b) claims that $z^*$ is the largest root of $R$ under the hazard condition alone, and this is false as printed: when $a - b(c-2s) + A < 0$ the maximum can sit at $z = A$, which need not be a root (finding F1; for $\varepsilon$ uniform on $[0,1]$, $a = 0.5$, $b = 1$, $c = 1$, $h = s = 0$, no root exists). Theorem 1(c)'s extra hypothesis is exactly what excludes this.

Formalization scope

The distribution of $\varepsilon$ is given by its density $f : \mathbb R \to \mathbb R$, with $f \ge 0$, $f = 0$ off $[A,B]$, $f$ interval-integrable on $[A,B]$ and $\int_A^B f = 1$; the bounded range $A < B$ is part of the model. All integrals are interval integrals. The expected profit is the integral of p. 184, not $\Psi - L$, and optimality ranges over $z \in [A,B]$ and every real price $p$. The hazard rate is $f(z)/[1-F(z)]$, which Lean evaluates to $0$ where $F(z) = 1$; the hazard condition is therefore required on $[A,B)$ instead of the printed $A \le z \le B$, since $r(B)$ is undefined. At each $z \in [A,B)$ the condition asks for a genuine derivative $r'$ of $r$ relative to $[A,B]$ (one-sided at $A$) with $2r(z)^2 + r' > 0$, never a junk derivative. The derivative $dE[\Pi(z,p(z))]/dz$ is Mathlib's deriv; the reduced objective is differentiable everywhere, which is part of a milestone. The milestones on $dR/dz$ and $d^2R/dz^2$ add continuity, respectively differentiability, of $f$ and $f(z) > 0$, which the paper uses silently; the goal adds nothing. The milestone on the endpoint values states $R(B) = -(c+h)$ without the printed "$< 0$".

These choices rule out the trivializing formalizations: an objective defined as $\Psi - L$ (which would make eq. (1) definitional and Lemma 1 a quadratic identity), optimality over $p(z)$ alone or over all real $z$, a hazard condition read through a junk derivative (which would reduce it to $2r^2 > 0$), and a goal restricted to uniform, exponential or IFR densities.

Beyond the paper. Two variants are posed beside the faithful statements. theorem_1b_corrected (corrected, finding F1) replaces the conclusion of Theorem 1(b) by "$z^*$ is $A$ or the largest root"; the uniform counterexample above shows the printed conclusion fails while the corrected one holds. theorem_1c_probe (sharpness probe, finding F2) is Theorem 1(c) with optimality over all real $z$; for $z < A$ the reduced objective is $b d^2 - (d+s)w + w^2/4b$ with $w = \mu - z$, which is unbounded, so a disproof is expected and shows that the restriction to $[A,B]$ cannot be dropped.

Contributions welcome: proofs of the milestones, general lemmas on interval integrals with moving limits against an integrable density, and a reusable unimodality lemma ("a $C^1$ function whose critical points are all strict local maxima rises, then falls").

Selected references

• N. C. Petruzzi, M. Dada, Pricing and the Newsvendor Problem: A Review with Extensions, Operations Research 47(2):183–194, 1999. doi:10.1287/opre.47.2.183
• E. S. Mills, Uncertainty and Price Theory, Quarterly Journal of Economics 73(1):116–130, 1959.
• T. M. Whitin, Inventory Control and Price Theory, Management Science 2(1):61–68, 1955.
• R. E. Barlow, F. Proschan, Statistical Theory of Reliability and Life Testing, Holt, Rinehart and Winston, 1975.

Motivation

Many queueing, inventory and storage models are Markov chains on a two-dimensional state space $(n, i)$: a level $n \ge 0$ (a queue length, a stock level) and a phase $i$ (the state of an arrival or service environment). When customers arrive in batches, the level can jump up by several units in one step but comes down one unit at a time; when services are grouped, it is the other way round. These two families, the M/G/1-type and GI/M/1-type chains, were developed by Neuts along separate lines (Neuts 1981; Neuts 1989, references 2 and 3), and they are the embedded chains of the classical M/G/1 queue at departures and the GI/M/1 queue at arrivals.

Chapter 13 of Latouche and Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling (SIAM 1999, doi:10.1137/1.9780898719734), shows that both families are embedded in quasi-birth-and-death processes (QBDs), and derives from the QBD theory of the earlier chapters the two central results: the first-passage matrix $H$ of an M/G/1-type chain is the minimal nonnegative solution of $H = \sum_n C_n H^n$, and the stationary distribution of a GI/M/1-type chain is matrix-geometric, $\tilde\pi_n = \tilde\pi_0 K^n$. This mission formalizes those statements.

Setting

Fix a countable set of phases $\iota$ and two sequences of real $\iota \times \iota$ matrices $D_0, D_1, \dots$ and $C_0, C_1, \dots$. States are pairs $(n, i)$, $n \ge 0$, $i \in \iota$; the level $\ell(n)$ is the set of states with first coordinate $n$. Row vectors act on the left: $\tilde\pi_n$ is the row vector $(\tilde\pi_{(n,i)})_{i}$.

The M/G/1-type chain (skip-free to the left) has transition matrix (13.1): from $(0,i)$ it moves to $(k,j)$ with probability $(D_k)_{ij}$; from $(n,i)$, $n \ge 1$, it moves to $(n-1+k, j)$ with probability $(C_k)_{ij}$, $k \ge 0$. The GI/M/1-type chain (skip-free to the right) has transition matrix (13.2): from $(n,i)$ it moves to $(0,j)$ with probability $(D_n)_{ij}$ and to $(n+1-k, j)$ with probability $(C_k)_{ij}$, $0 \le k \le n$. Both matrices are assumed stochastic: nonnegative entries, rows summing to one.

The M/G/1-type chain is embedded in a QBD (13.3) on the states $(n, d, j)$: the level is $n$ and the phase is the pair $(d, j)$, where $d \ge 0$ counts the levels still to be climbed. A jump from $(n,i)$ to $(n+k,j)$ becomes the sequence $(n,0,i) \to (n,k,j) \to (n+1,k-1,j) \to \dots \to (n+k,0,j)$, with probabilities $(C_{k+1})_{ij}, 1, \dots, 1$ (Table 13.1, blocks (13.4)–(13.7)). In this QBD,

$$H_{ij} = \mathrm P[\tau < \infty \ \&\ (X_\tau, D_\tau, J_\tau) = (n-1, 0, j) \mid (X_0, D_0, J_0) = (n, 0, i)], \tag{13.9}$$

with $\tau$ the first passage time from $\ell(n)$ to $\ell(n-1)$, and $G$ is the QBD's matrix of first-passage probabilities from $\ell(n)$ to $\ell(n-1)$.

A state is positive recurrent when its expected return time is finite, and a stationary probability vector $\tilde\pi$ is nonnegative, sums to one and satisfies $\tilde\pi \tilde P = \tilde\pi$.

Formalization targets

Goal: $H$ is the minimal solution of (13.10) (Theorem 13.1.4)

For a stochastic M/G/1-type chain and every level $n \ge 1$, the matrix $H$ of (13.9) has finite entries and is the minimal nonnegative solution of

$$H = \sum_{n \ge 0} C_n H^n, \tag{13.10}$$

that is, it solves (13.10) and lies entrywise below every nonnegative real matrix $Y$ with $Y = \sum_n C_n Y^n$.

Milestones

1. Theorem 13.1.3. The QBD's first-passage matrix has the block structure (13.8): $G_{(k,i),(d,j)} = (H^{k+1})_{ij}$ when $d = 0$, and $0$ otherwise.
2. Corollary 13.1.2. If $\tilde\pi$ is the stationary probability vector of an irreducible, positive recurrent M/G/1-type chain and $\pi$ is a stationary probability vector of the QBD, then $\tilde\pi_i = c\,\pi_{i,0}$ for a constant $c$.
3. Theorem 13.1.5. With $\bar C_n = \sum_{i \ge 0} C_{n+i} H^i$ and $\bar D_n = \sum_{i\ge 0} D_{n+i}H^i$, the matrix $(I - \bar C_1)^{-1} = \sum_\nu \bar C_1^{\,\nu}$ is finite and, for $n \ge 1$,

$$\tilde\pi_n = \Big(\tilde\pi_0 \bar D_n + \sum_{1 \le j \le n-1} \tilde\pi_j \bar C_{n-j+1}\Big)(I - \bar C_1)^{-1}. \tag{13.12}$$

4. Theorem 13.2.3. For an irreducible, positive recurrent GI/M/1-type chain with stationary probability vector $\tilde\pi$, there is a nonnegative $K$ with

$$\tilde\pi_n = \tilde\pi_0 K^n, \tag{13.26}$$

$K$ is the minimal nonnegative solution of $K = C_0 + KC_1 + K^2C_2 + \cdots$ (13.27), and, as printed, $K_{ij}$ is for every $n\ge1$ the expected number of visits to $(n+1, j)$ before the first return to $\ell(n-1)$, starting from $(n,i)$.

Significance

Theorem 13.1.4 turns a probabilistic object, the matrix of first-passage probabilities one level down, into the minimal solution of a matrix equation of the size of one level. Every numerical method for M/G/1-type chains, including the quadratically convergent algorithms of Bini and Meini, computes this solution. Theorem 13.1.5 (Ramaswami's formula) then gives the stationary distribution recursively from $H$. Theorem 13.2.3 is the matrix-geometric property of GI/M/1-type chains, the origin of the matrix-analytic approach. Together these results cover the stationary analysis of batch-arrival and batch-service queues.

All four results are classical and proved in the book. None of them is formalized. Lean's Mathlib has no theory of first-passage matrices or of block-structured Markov chains. On Prove2Me the closest results are scalar ones: the geometric law of the G/M/1 arrival chain (QueueingFundamentals.GM1.arrival_point_geometric), the M/G/1 departure chain (QueueingFundamentals.MG1.*) and Foster-type classification of GI/M/1 chains (FosterQueues.GIM1.*). The matrix statements here contain those as the case of one phase.

Difficulty

The equation (13.10) follows formally from a first-step decomposition once $H$ is known to be well defined independent of the level. The hard part is the passage between the chain and its QBD embedding. A first passage of the QBD from $\ell(n+k)$ down to $\ell(n-1)$ splits into $k+1$ successive one-level passages, a strong-Markov decomposition that must here be carried out as an identity of sums over trajectories. Minimality cannot be read off the equation: it needs the comparison of $H$ with the iterates of the map $Y \mapsto \sum_n C_n Y^n$ started from $0$, with every series over a possibly infinite phase set. For Theorem 13.1.5 the inverse $(I - \bar C_1)^{-1}$ has to be shown finite, which uses recurrence. Reading the interpretive clause of Theorem 13.2.3 literally is a trap (see below).

Formalization scope

Everything lives in the namespace MatrixAnalytic.SkipFree. Chains are real matrices ℕ × ι → ℕ × ι → ℝ on a countable phase type ι, and nothing requires finitely many phases. The QBD embedding has phase type ℕ × ι. Levels are counted from $0$. Stochasticity, irreducibility, positive recurrence (finite expected return time) and stationary probability vectors restate the platform definition mm_countable. Probabilities and expected numbers of visits ($H$, $G$, and the visits in (13.26)) are sums over finite trajectories of products of one-step transition probabilities, valued in $[0,\infty]$. Products, powers and series of nonnegative matrices are also computed in $[0,\infty]$. The book's $(I - M)^{-1}$ is $\sum_\nu M^\nu$, and finiteness is a conclusion wherever the book asserts it. "Minimal" is entrywise, over nonnegative real matrices.

$H$ and $K$ are never defined as solutions of (13.10) or (13.27), nor as "the minimal nonnegative solution". That definition would make Theorems 13.1.4 and 13.2.3 true by construction. $H$ is the first-passage probability (13.9), and $K$ is existentially quantified with its probabilistic meaning stated.

Beyond the paper. Finding F1/F2: the interpretive clause of Theorem 13.2.3 names the wrong taboo set. With one phase, $C_0 = 0.2$, $C_1 = 0.3$, $C_2 = 0.5$, the minimal solution of (13.27) is $K = 0.4$, but the expected number of visits to $n+1$ from $n$ before reaching $n-1$ is $0.8$. The faithful item states the page and is expected to be disproved. The corrected variant gm1_matrix_geometric_corrected counts visits before the return to $\ell(0) \cup \dots \cup \ell(n)$, which is what the proof's $K = R^{(1)}_{00}$ records.

Infrastructure needed: path-sum decompositions (first step, last exit, concatenation at a first-passage time) for countable chains, monotone iteration for matrix power series in $[0,\infty]$, and censoring of a chain on a subset of states. All of these are reusable for the other chapters of this series. Proofs that go directly through the skip-free chain, without the QBD embedding, are welcome.

Selected references

1. G. Latouche and V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, ASA-SIAM Series on Statistics and Applied Probability, SIAM, 1999. doi:10.1137/1.9780898719734
2. M. F. Neuts, Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach, Johns Hopkins University Press, 1981; Dover reprint 1994.
3. M. F. Neuts, Structured Stochastic Matrices of M/G/1 Type and Their Applications, Marcel Dekker, 1989.
4. V. Ramaswami, "Stable recursion for the steady state vector for Markov chains of M/G/1 type", Communications in Statistics. Stochastic Models 4 (1988) 183–188 (reference [94] of the book; cited there in Remark 13.1.6).
5. D. A. Bini and B. Meini, "On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems", in W. J. Stewart (ed.), Computations with Markov Chains, Kluwer, 1995, pp. 21–38 (reference [6] of the book).