Numerical Analysis

Motivation

Newton's method is the basic second-order method of continuous optimization: at the current point it replaces the objective by its second-order Taylor model and jumps to the stationary point of that model. Its defining property is speed near a nondegenerate minimum, where the error is squared at every step, so that the number of correct digits roughly doubles per iteration. This local behaviour is what makes Newton's method the inner engine of interior point methods, the polynomial-time algorithms for linear, conic and general convex programming (Nesterov and Nemirovski, 1994). In S. Bubeck's monograph Convex Optimization: Algorithms and Complexity (Foundations and Trends in Machine Learning, 2015; arXiv:1405.4980v2), §5.3.2 recalls the traditional local analysis of Newton's method, Theorem 5.3, before turning to the affine-invariant self-concordance analysis used for interior point methods. This mission formalizes that theorem and the four steps of its proof.

Setting

Let $\mathbb R^n$ carry the Euclidean norm $\|\cdot\|$, and write $\|A\|$ for the operator norm of a linear map $A:\mathbb R^n\to\mathbb R^n$, so that $\|Ax\|\le\|A\|\,\|x\|$. Let $f:\mathbb R^n\to\mathbb R$ be a $C^2$ function, with gradient $\nabla f(x)\in\mathbb R^n$ and Hessian $\nabla^2 f(x)$, a linear map $\mathbb R^n\to\mathbb R^n$ (the derivative of the gradient map). For a real number $c$, $A\succeq cI_n$ means $\langle Av,v\rangle\ge c\|v\|^2$ for all $v\in\mathbb R^n$.

The Hessian is $M$-Lipschitz if $\|\nabla^2 f(x)-\nabla^2 f(y)\|\le M\|x-y\|$ for all $x,y\in\mathbb R^n$.

Newton's method starts at $x_0\in\mathbb R^n$ and iterates, for $k\ge0$,

$$x_{k+1}=x_k-[\nabla^2 f(x_k)]^{-1}\nabla f(x_k).$$

A point $x^*$ is a local minimum of $f$ if $f(x^*)\le f(x)$ for all $x$ in a neighbourhood of $x^*$; it has strictly positive Hessian if $\nabla^2 f(x^*)\succeq\mu I_n$ for some $\mu>0$.

Formalization targets

Goal: Theorem 5.3 (p. 320)

Assume the Hessian of $f$ is $M$-Lipschitz, $M>0$, and $x^*$ is a local minimum with $\nabla^2 f(x^*)\succeq\mu I_n$, $\mu>0$. If $\|x_0-x^*\|\le\mu/(2M)$, then Newton's method from $x_0$ is well defined (every Hessian along the iterates is invertible, so the sequence exists and is unique) and

$$\|x_{k+1}-x^*\|\le\frac M\mu\,\|x_k-x^*\|^2\quad(k\ge0),\qquad x_k\to x^*.$$

Milestones (p. 321, the steps of the proof)

1. The integral formula $\int_0^1\nabla^2 f(x+sh)\,h\,ds=\nabla f(x+h)-\nabla f(x)$.
2. The error representation of one Newton step, $x_{k+1}-x^*=[\nabla^2 f(x_k)]^{-1}\int_0^1[\nabla^2 f(x_k)-\nabla^2 f(x^*+s(x_k-x^*))](x_k-x^*)\,ds$.
3. The Lipschitz bound $\int_0^1\|\nabla^2 f(x_k)-\nabla^2 f(x^*+s(x_k-x^*))\|\,ds\le\frac M2\|x_k-x^*\|$.
4. The Hessian lower bound $\nabla^2 f(x_k)\succeq(\mu-M\|x_k-x^*\|)I_n\succeq\frac\mu2I_n$ when $\|x_k-x^*\|\le\mu/(2M)$.

Significance

The theorem gives a quantitative basin of quadratic convergence: an explicit radius $\mu/(2M)$, depending only on the curvature at the minimum and the Lipschitz constant of the Hessian, inside which Newton's method needs only $O(\log\log(1/\varepsilon))$ iterations to reach accuracy $\varepsilon$. It is the classical statement whose shortcomings (dependence on a choice of norm, constants that change under linear changes of variables) motivate the self-concordance theory of the following subsections, and it is the local convergence result invoked whenever a damped or globalized Newton scheme is shown to enter its quadratic phase.

On the formal side, Mathlib has the calculus this needs (Fréchet derivatives, interval integrals of vector-valued maps, operator norms) but no convergence theorem for multivariate Newton's method for minimization. A formal proof produces reusable pieces: the integral form of the mean value theorem for gradients, the stability of a positive-definite lower bound under Lipschitz perturbations, and an inverse-operator norm bound from a quadratic-form lower bound. The result itself is classical and fully proved in the literature; what is open here is its machine-checked proof in this form.

Difficulty

The individual inequalities are short, but the argument is an induction in which well-definedness and the rate are proved together: the Hessian at $x_k$ is invertible only because $x_k$ is still in the ball of radius $\mu/(2M)$, and $x_{k+1}$ stays in that ball only because of the rate. A proof that first assumes the sequence exists and then bounds it is circular. The proof also passes between two kinds of control on the Hessian, a lower bound on its quadratic form and an operator-norm bound on its inverse, and the second is only meaningful once invertibility is established. Finally, the integral manipulations need integrability of the maps $s\mapsto\nabla^2 f(x^*+s(x_k-x^*))(x_k-x^*)$, which comes from the continuity of the Hessian.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n). The gradient and Hessian are explicit maps $g:\mathbb R^n\to\mathbb R^n$ and $H:\mathbb R^n\to(\mathbb R^n\to_L\mathbb R^n)$ with ContDiff ℝ 2 f, HasGradientAt f (g x) x and HasFDerivAt g (H x) x at every point; the norm on $H(x)$ is Mathlib's operator norm, as on the page. $A\succeq cI_n$ is the quadratic-form inequality. A Newton run is a sequence $x:\mathbb N\to\mathbb R^n$ indexed from $0$ satisfying the linear system $\nabla^2 f(x_k)(x_k-x_{k+1})=\nabla f(x_k)$; no inverse of a possibly singular operator appears in any hypothesis, and "well defined" is a conclusion: a unique run exists from $x_0$ and every Hessian along it is bijective. The rate and $x_k\to x^*$ are asserted for every run. The error representation is stated with both sides multiplied by $\nabla^2 f(x_k)$, which is equivalent to the printed form once the Hessian is invertible. Milestones 3 and 4 use only the Lipschitz property and are stated for any Lipschitz map $H$.

Added hypothesis: $M>0$ (the radius $\mu/(2M)$ divides by $M$; with $M=0$, Lean's convention $\mu/0=0$ would collapse the hypothesis to $x_0=x^*$). Convexity of $f$ is not assumed, as on the page; $x^*$ is a local minimum and $\nabla f(x^*)=0$ is derived, not assumed. Encoding the Newton step with Lean's inverse (which returns $0$ on singular maps), or replacing $\nabla^2 f(x^*)\succeq\mu I_n$ by mere invertibility, would change the theorem and is ruled out.

A complete development needs the fundamental theorem of calculus for $C^1$ vector-valued maps along segments, Hessian-based quadratic-form estimates, and operator-norm bounds for inverses; all are reusable for the analysis of damped Newton, cubic regularization and interior point methods. Proofs of the milestones independently of the goal are welcome.

Selected references

• S. Bubeck, Convex Optimization: Algorithms and Complexity, Foundations and Trends in Machine Learning 8(3–4):231–358, 2015. arXiv:1405.4980v2, §5.3.2, Theorem 5.3, pp. 320–321.
• Yu. Nesterov and A. Nemirovski, Interior-Point Polynomial Algorithms in Convex Programming, SIAM Studies in Applied Mathematics 13, 1994. doi:10.1137/1.9781611970791
• Yu. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer, 2004, Theorem 1.2.5. doi:10.1007/978-1-4419-8853-9

Motivation

Mean field games (MFG), introduced by Lasry and Lions and by Huang, Caines and Malhamé, describe Nash equilibria of very large populations of identical rational agents. The equilibrium is a coupled system: a backward Hamilton–Jacobi–Bellman equation for the value function $u$ of a representative agent and a forward Fokker–Planck equation for the density $m$ of the population. In the planning problem, proposed by P.-L. Lions, both the initial density $m_0$ and the final density $m_T$ are prescribed, and one asks for a cost structure under which the population moves from one to the other. This is a mean field analogue of optimal transport.

Achdou, Camilli and Capuzzo-Dolcetta (hal-00465404, SIAM J. Control Optim. 2012) propose finite-difference schemes for the planning problem. Mission I of this series treats the existence of a solution of the discrete planning scheme. Since the two boundary conditions on $m$ make the discrete system hard to solve directly, the paper also introduces a penalized scheme, in which the initial condition $M^0 = m_0$ is replaced by a penalty $U^0 = (M^0 - m_0)/\varepsilon$; for each $\varepsilon > 0$ this is a standard discrete MFG system with a unique solution. This mission formalizes the paper's §3.2: the penalized solutions converge, as $\varepsilon \to 0$, to a solution of the planning scheme.

Setting

Fix integers $N_h \ge 1$, $N_T \ge 1$, a horizon $T > 0$ and a viscosity $\nu \ge 0$; put $h = 1/N_h$ and $\Delta t = T/N_T$. The grid $\mathbb T^2_h$ consists of the points $x_{i,j}$, $(i,j) \in (\mathbb Z/N_h\mathbb Z)^2$ (periodic indices). A grid function is a real function on $\mathbb T^2_h$; time levels are $n = 0, \dots, N_T$.

• $(D_1^+U)_{i,j} = (U_{i+1,j} - U_{i,j})/h$, $(D_2^+U)_{i,j} = (U_{i,j+1} - U_{i,j})/h$, and the discrete gradient is $[D_hU]_{i,j} = ((D_1^+U)_{i,j}, (D_1^+U)_{i-1,j}, (D_2^+U)_{i,j}, (D_2^+U)_{i,j-1}) \in \mathbb R^4$.
• $\Delta_h$ is the five-point Laplacian.
• A numerical Hamiltonian $g(x_{i,j}, q)$, $q \in \mathbb R^4$, is monotone (G1), $C^1$ (G3), convex (G4) and coercive (G5) in $q$.
• $\mathcal B_{i,j}(U, M)$ is the discrete transport term $\mathrm{div}_h\big(M\,\nabla_q g(\cdot, [D_hU])\big)$.
• $W : \mathbb R \to \mathbb R$ is strictly convex, superlinear and $C^2$, with $V = W'$ (hypothesis (24)).
• $\mathcal K = \{M \ge 0 : h^2\sum_{i,j} M_{i,j} = 1\}$ is the set of discrete probability densities, and $m_0, m_T \in \mathcal K$ with $m_0 > 0$.

The planning scheme (18) asks for families $(U^n, M^n)$ with, for $n < N_T$,

$$\frac{U^{n+1} - U^n}{\Delta t} - \nu\Delta_hU^{n+1} + g(x, [D_hU^{n+1}]) = V(M^n), \qquad \frac{M^{n+1} - M^n}{\Delta t} + \nu\Delta_hM^n + \mathcal B(U^{n+1}, M^n) = 0,$$

$M^n \in \mathcal K$, $M^0 = m_0$ and $M^{N_T} = m_T$. The penalized scheme (20)–(23) has the same two equations, $M^n \in \mathcal K$ for $n < N_T$, $M^{N_T} = m_T$, and $U^0 = (M^0 - m_0)/\varepsilon$ in place of $M^0 = m_0$.

Formalization targets

Goal: Proposition 4

Under the hypotheses of Theorem 1, for every sequence $\varepsilon_k \downarrow 0$, every choice of penalized solutions $(U^{\varepsilon_k}, M^{\varepsilon_k})$, and every limit $M^{\varepsilon_k} \to M$ in $\mathcal K^{N_T+1}$,

$$\exists\, U,\ \text{a subsequence with } U^{\varepsilon_k} \to U, \quad (U, M) \text{ solves (43)–(46)}.$$

If $g$ is strictly convex, the whole sequence $(U^{\varepsilon_k}, M^{\varepsilon_k})$ converges to the unique solution of (43)–(46) with $\sum_{i,j} U^0_{i,j} = 0$.

Milestones

1. (13): (G5) implies $\max_{i,j} g(x_{i,j}, [D_hU]_{i,j}) / \|[D_hU]\|_\infty \to +\infty$.
2. Theorem 2: the penalized control problem (50), $\min \Theta^*(M,Z) + \frac{1}{2\varepsilon\Delta t}\sum(M^0 - m_0)^2$ under the discrete Fokker–Planck constraint, has a minimizer whose optimality system is (20)–(23).
3. Proposition 2: $\max_{i,j}|M^{\varepsilon,0}_{i,j} - (m_0)_{i,j}| \le C\varepsilon^{1/2}$.
4. Proposition 3: $\max_{n,i,j}|U^{\varepsilon,n}_{i,j}| \le C$.
5. Corollary 1: $\max_{i,j}|M^{\varepsilon,0}_{i,j} - (m_0)_{i,j}| \le C\varepsilon$.
6. Proposition 1: uniqueness of $M$ for (18), and of $U$ normalized by $\sum U^0 = 0$ when $g$ is strictly convex.

In 3–5 the constant $C$ depends on $h$, $\Delta t$ and the data, never on $\varepsilon$.

Significance

Proposition 4 justifies the penalized scheme as a way to compute solutions of the discrete planning problem. The penalized system has a standard initial–terminal structure that Newton's method (the paper's §4) handles well. The planning system does not, because it prescribes two conditions on $M$ and none on $U$. The estimates of Propositions 2–3 and Corollary 1 also give a rate: the initial mismatch of the penalized density is $O(\varepsilon)$.

The results are proved in the paper. None of them has been machine-checked: the platform holds no formalization of mean field games, of their finite-difference schemes, or of the Lasry–Lions monotonicity argument. A complete development would yield a reusable library of the monotone finite-difference operators on the discrete torus, a formal Lasry–Lions uniqueness argument for discrete MFG systems, and a template for compactness-plus-uniqueness convergence proofs in finite dimension.

Difficulty

Compactness of the densities is free, since $\mathcal K^{N_T+1}$ is compact. Passing to the limit in the scheme needs only continuity of $g$, $\nabla_q g$ and $V$. The substance lies in the two uniform estimates.

• Bounding $U^\varepsilon$ uniformly in $\varepsilon$ (Proposition 3) is the central difficulty. The initial condition $U^{\varepsilon,0} = (M^{\varepsilon,0} - m_0)/\varepsilon$ involves division by $\varepsilon$, so the obvious bound $|U^{\varepsilon,0}| \le 2/(h^2\varepsilon)$ blows up. The comparison principle for the discrete HJB equation, which is the natural first idea, cannot repair this: it propagates whatever bound $U^0$ has. The HJB equation alone does not bound $U^\varepsilon$; the coupling with the Fokker–Planck equation, the coercivity (13), and a lower bound on $M^{\varepsilon,0}$ that is uniform in $\varepsilon$ (which is where Proposition 2 enters) all have to be used.
• Proposition 2 compares the penalized problem with the planning problem through the convex duality of Theorem 2. It therefore needs a solution of the planning scheme, which is mission I's goal.
• Proposition 1 needs the discrete Lasry–Lions identity, a summation by parts across the coupled system.

Formalization scope

All objects live in the namespace MFGPlanning.Penalized.

• Grid points are ZMod Nh × ZMod Nh (periodicity is built in), time levels are Fin (NT + 1), and the four momentum components $q_1, \dots, q_4$ are Fin 4 indices 0, …, 3.
• The data structure carries $N_h, N_T \ge 1$, $T > 0$ and $\nu \ge 0$.
• $g$ is given only at the grid points. (G2), which only defines the continuous Hamiltonian, is not encoded.
• "Coercive" in (24) is read as superlinear.
• $V_h[M] = V(M_{i,j})$, as (24) prescribes.
• $\Theta^*$ and $(W+\chi)^*$ take values in EReal, so unbounded suprema are $+\infty$.
• Convergence is in the product topology, which on these finite-dimensional spaces is $\max_n \|\cdot\|_\infty$ convergence.
• "The unique solution" of (43)–(46) means unique under the normalization $\sum_{i,j} U^0_{i,j} = 0$, because $U$ is otherwise determined only up to a constant.

Three trivializing formalizations are ruled out. The constants in Propositions 2–3 and Corollary 1 are quantified before $\varepsilon$ and before the solution. The goal quantifies over every sequence $\varepsilon_k \to 0$ and assumes no convergence of $U$ and no limit solution. Its part 2 does not assume convergence of $M$.

The existence of a solution of the planning scheme (Theorem 1) is posed in mission I and not restated here. A solver of Proposition 2 will need it once mission I's goal is proved. The existence and uniqueness of solutions of (20)–(23) are quoted by the paper from Achdou–Capuzzo-Dolcetta (2010) and are not items. Welcome contributions include discrete summation-by-parts lemmas on the torus, the discrete comparison principle for monotone schemes, and the Poincaré-type inequality $\|W\|_\infty \le c\|[D_hW]\|_\infty$ on zero-mean grid functions.

Selected references

• Y. Achdou, F. Camilli, I. Capuzzo-Dolcetta, Mean field games: numerical methods for the planning problem, HAL preprint hal-00465404v1, 2010; SIAM J. Control Optim. 50 (2012). https://hal.science/hal-00465404v1, https://doi.org/10.1137/100790069
• Y. Achdou, I. Capuzzo-Dolcetta, Mean field games: numerical methods, SIAM J. Numer. Anal. 48 (2010) 1136–1162. https://doi.org/10.1137/090758477
• J.-M. Lasry, P.-L. Lions, Mean field games, Japanese Journal of Mathematics 2 (2007) 229–260. https://doi.org/10.1007/s11537-007-0657-8
• M. Huang, R. P. Malhamé, P. E. Caines, Large population stochastic dynamic games: closed-loop McKean–Vlasov systems and the Nash certainty equivalence principle, Communications in Information and Systems 6 (2006) 221–252. https://doi.org/10.4310/CIS.2006.v6.n3.a5

Motivation

Mean field games (Lasry and Lions, 2006–2007; Huang, Malhamé and Caines, 2006) model the limit of a large population of identical rational agents. Each agent solves an optimal control problem whose cost depends on the distribution $m$ of all agents, and the distribution is in turn transported by the agents' optimal feedback. In the continuous setting this gives a coupled system: a backward Hamilton–Jacobi equation for the value function $u$ and a forward Fokker–Planck equation for the density $m$.

In the usual formulation, $m$ is prescribed at the initial time and $u$ at the final time. The planning problem, introduced by P.-L. Lions in his Collège de France lectures, prescribes instead both the initial density $m_0$ and the final density $m_T$, and asks for a cost (through $u$) that steers the population from one to the other. According to the paper (§1, pp. 2–3), Lions proved existence for the continuous planning problem in mainly two cases: $\nu = 0$ with a smooth, strictly convex, superlinear Hamiltonian; and $\nu > 0$ with $H(p) = c|p|^2$ or close to it. In both cases the coupling is local and the densities are smooth and bounded away from $0$. Existence for $\nu > 0$ and more general Hamiltonians was then open, and for sublinear $H$, $m_0\ne m_T$ and short horizons there is no solution.

Achdou, Camilli and Capuzzo-Dolcetta (hal-00465404, 2010; SIAM J. Control Optim. 2012) introduced a finite-difference scheme for the planning problem and proved that the discrete system has a solution. The proof writes the scheme as the optimality system of a discrete optimal control problem, a Fokker–Planck equation driven by a control, and obtains a solution from a saddle point given by the Fenchel–Rockafellar duality theorem. This mission formalizes that existence result: Theorem 1 of §3.1 together with the lemmas on which its proof rests.

Setting

Fix integers $N_h, N_T \ge 1$, a horizon $T > 0$ and a viscosity $\nu \ge 0$; let $h = 1/N_h$ and $\Delta t = T/N_T$. The grid $\mathbb T^2_h$ is the periodic $N_h\times N_h$ grid on the two-dimensional torus, with points $x_{i,j}$, $(i,j)\in(\mathbb Z/N_h)^2$. On grid functions $U$ the scheme uses the forward differences $D_1^+$, $D_2^+$, the four-component discrete gradient $[D_hU]_{i,j} = ((D_1^+U)_{i,j}, (D_1^+U)_{i-1,j}, (D_2^+U)_{i,j}, (D_2^+U)_{i,j-1})$, the five-point Laplacian $\Delta_h$, a discrete divergence $\mathrm{div}_h$ of four-component fields, and the transport operator $\mathcal B(U, M) = \mathrm{div}_h(M\,\nabla_q g(\cdot, [D_hU]))$.

A numerical Hamiltonian $g(x_{i,j}, q_1,q_2,q_3,q_4)$ is monotone (nonincreasing in $q_1, q_3$, nondecreasing in $q_2, q_4$), $C^1$, convex, and superlinearly coercive in the directions where monotonicity does not bound it. The coupling is local: $V = W'$, with $W$ strictly convex, superlinear and $C^2$. The set $\mathcal K$ consists of the discrete probability densities, $h^2\sum_{i,j}M_{i,j} = 1$, $M\ge 0$. The discrete planning scheme (18) asks for $(U^n, M^n)_{0\le n\le N_T}$ with

$$\frac{U^{n+1}-U^n}{\Delta t} - \nu\Delta_hU^{n+1} + g(x,[D_hU^{n+1}]) = V(M^n),\qquad \frac{M^{n+1}-M^n}{\Delta t} + \nu\Delta_hM^n + \mathcal B(U^{n+1},M^n) = 0,$$

for $0\le n<N_T$, with $M^n\in\mathcal K$, $M^0 = m_0$ and $M^{N_T} = m_T$.

The duality is set up as follows. With $\chi$ the indicator of $\{m\ge0\}$, let $\Theta(\alpha,\beta) = \sum_{n,i,j}(W+\chi)^*(\alpha^n_{i,j} + g(x_{i,j},[\beta^n]_{i,j}))$ on dual variables $(\alpha^n,\beta^n)_{1\le n\le N_T}$. Let $\Lambda(\Psi)$ be the linear map sending $\Psi = (\Psi^n)_{0\le n\le N_T}$ to the discrete Hamilton–Jacobi operator and the discrete gradient of $\Psi^{n+1}$. Let $\Sigma(\alpha,\beta) = \mathcal F(\Psi)$ if $(\alpha,\beta) = \Lambda(\Psi)$ with $\sum\Psi^0 = 0$, and $+\infty$ otherwise, where $\mathcal F(\Psi) = \frac1{\Delta t}(\sum m_0\Psi^0 - \sum m_T\Psi^{N_T})$. The Legendre–Fenchel transforms $\Theta^*$, $\Sigma^*$ act on primal variables $(M^n, Z^n)_{0\le n<N_T}$, where $M^n$ is paired with $\alpha^{n+1}$ (Remark 2).

Formalization targets

Goal: Theorem 1

Under (G1), (G3)–(G5), (24), $m_0, m_T\in\mathcal K$, $m_0 > 0$, and either $\nu > 0$ or ($\nu = 0$ and $m_T > 0$):

$$\min_{M,Z}\ \Theta^*(M,Z) + \Sigma^*(-M,-Z) = -\min_{\alpha,\beta}\big(\Theta(\alpha,\beta) + \Sigma(\alpha,\beta)\big)$$

has a solution $(M,Z)$, $(\alpha,\beta)$ with a finite common value. Moreover $(\alpha,\beta) = \Lambda(U)$ for some $U$, and $(U, M)$, with $M^{N_T} = m_T$ appended, solves the scheme (18), with $Z^{k,n} = M^n\,\partial_{q_k}g(x,[D_hU^{n+1}])$.

Milestones, in attack order

1. §3.1, p. 7. $V$ maps $(0,\infty)$ onto $(\lambda,\infty)$; $(W+\chi)^*$ is finite, convex, continuous and nondecreasing, with explicit values on and off $\mathcal J_V$.
2. Lemma 1. $\Theta$ is convex and continuous, $\Sigma$ is convex and l.s.c., and $\Sigma$ and $\Theta$ are both finite at some point.
3. Lemma 2. $\Theta^*$ and $\Sigma^*$ are convex and l.s.c., with explicit formulas.
4. (30). $\Sigma^*(-M,-Z)$ is $0$ on the constraint set of the control problem (26) and $+\infty$ off it.
5. Lemma 3. If $m_0 > 0$, some $(M,Z)$ has $\Theta^*$, $\Sigma^*(-M,-Z)$ finite and $\Theta^*$ finite and continuous near it.
6. Positivity. A discrete strong maximum principle from the proof of Theorem 1: densities in $\mathcal K$ that solve the discrete Fokker–Planck equation (42) are strictly positive before the final time.

Significance

Theorem 1 is the existence result for the finite-difference planning problem. It is used in the paper's second part (§3.2), where solutions of a penalized scheme, in which the final condition is relaxed into a penalty, are shown to converge to a solution of (18) as the penalty parameter vanishes. The penalized scheme is what is solved numerically. The discrete existence result covers general convex, monotone, coercive numerical Hamiltonians and every $\nu\ge0$, which includes discretizations of the continuous cases that were still open when the paper was written.

The result has a published proof. What this mission adds is a machine-checked version of it: a complete discrete model of the planning problem (grid, operators, scheme, duality functionals) and the convex-analytic chain that connects it to the Fenchel–Rockafellar theorem. As far as a search of the platform shows, no part of this has been formalized; the series' second mission formalizes the convergence of the penalized scheme on the same model.

Difficulty

Existence for a coupled forward–backward nonlinear system with conditions at both ends of the time interval is not reachable by a fixed-point or time-marching argument: the Hamilton–Jacobi equation runs backward and the Fokker–Planck equation forward, and the final density is imposed rather than computed. The approach goes through duality, and three points are delicate.

1. All the functionals in the duality take the value $+\infty$, and Fenchel–Rockafellar needs a qualification condition on each side. Lemma 1 gives it for the dual problem; Lemma 3 gives it for the primal one, and needs $m_0 > 0$ and the coercivity (G5).
2. The optimality conditions only give a complementarity system: the Hamilton–Jacobi equation holds where $M^n > 0$ and becomes an inequality where $M^n = 0$. Recovering the scheme (18) needs the strict positivity of $M^n$, a discrete strong maximum principle. That principle uses the monotonicity (G1) and either diffusion or a positive final density.
3. The bookkeeping of the time lag between primal and dual variables and of the discrete integration by parts in $\Sigma^*$ must be exact.

Formalization scope

The Lean development lives in the namespace MFGPlanning.Existence. A structure Data carries $N_h, N_T \ge 1$, $T > 0$, $\nu\ge0$, $g$ at the grid points, $W$, $m_0$ and $m_T$. Committed conventions:

• grid indices are ZMod Nh × ZMod Nh (periodic), and $d = 2$ as in the paper;
• the components $q_1,\dots,q_4$ and $Z^1,\dots,Z^4$ are the Fin 4 indices 0..3;
• time levels are Fin (NT+1) for $U$ and the scheme, and Fin NT for the duality variables. M k is $M^k$ and α k is $\alpha^{k+1}$ (Remark 2), and Fin.snoc M mT appends $M^{N_T} = m_T$;
• $(W+\chi)^*$, $\Theta$, $\Theta^*$, $\Sigma$, $\Sigma^*$ are EReal-valued lattice suprema and infima, so unbounded suprema are $+\infty$ and not a junk value. Convexity of an extended-valued functional is convexity of its epigraph.

Disclosed encodings:

• (G2) is not encoded, since it only defines the continuous Hamiltonian;
• "coercive" in (24) is read as superlinear growth $W(m)/|m|\to\infty$, which the paper's own consequence $V((0,\infty)) = (\lambda,\infty)$ requires;
• $g$ is given only at grid points;
• the optimality conditions (32)–(33) are stated through what Theorem 1 says they are equivalent to: the scheme (18) and the relation (39).

Theorem 1 is not reducible to "the scheme (18) has a solution". The goal also asserts that the primal and dual problems attain their minima, that there is no duality gap with a finite value, and that $(\alpha,\beta) = \Lambda(U)$. The duality functionals are never real-valued suprema, which would make (30) and (31) hold or fail for junk reasons. A complete development needs finite-dimensional convex analysis on extended-valued functions: conjugates, the Fenchel–Rockafellar theorem with attainment, and subdifferential optimality conditions. These parts are reusable well beyond this mission, and proofs of them as separate theorems are welcome.

Selected references

• Y. Achdou, F. Camilli, I. Capuzzo-Dolcetta, Mean field games: numerical methods for the planning problem, preprint hal-00465404v1, 2010. https://hal.science/hal-00465404 ; published in SIAM J. Control Optim. 50(1), 2012. https://doi.org/10.1137/100790069
• Y. Achdou, I. Capuzzo-Dolcetta, Mean field games: numerical methods, SIAM J. Numer. Anal. 48(3), 2010. https://doi.org/10.1137/090758477
• J.-M. Lasry, P.-L. Lions, Mean field games, Japanese Journal of Mathematics 2(1), 2007. https://doi.org/10.1007/s11537-007-0657-8
• I. Ekeland, R. Temam, Convex Analysis and Variational Problems, North-Holland, 1976 (SIAM Classics reprint 1999). https://doi.org/10.1137/1.9781611971088

Motivation

Optimization problems with orthonormality constraints arise throughout numerical linear algebra and its applications: computing a few eigenvectors or singular vectors, Procrustes problems in statistics and shape analysis, orthogonal factor rotation, independent component analysis, and the orthogonality constraints of electronic-structure calculations. The feasible set of such a problem is the Stiefel manifold of orthonormal $m$-frames in $\mathbb R^n$. Riemannian optimization algorithms on this manifold (Riemannian gradient, Newton and trust-region methods; see Absil, Mahony and Sepulchre, Optimization Algorithms on Matrix Manifolds, 2008) take a step in a tangent direction and then need a retraction: a map that brings the updated point back onto the manifold while agreeing with the geometry to first order.

The most natural way to come back to a constraint set is to take the nearest point. P.-A. Absil and J. Malick, Projection-like retractions on matrix manifolds (SIAM J. Optim. 22(1), 2012; the mission follows the authors' version, HAL hal-00651608v2), show in their Proposition 3.2 that the metric projection onto a smooth submanifold yields a retraction, and then work out the projection explicitly for several matrix manifolds. For the Stiefel manifold, their Proposition 3.4 identifies the projection with a factor of the singular value decomposition, equivalently with the orthonormal factor of the polar decomposition. The paper notes that this result was mentioned without proof in Higham's survey of matrix nearness problems, and gives a proof along the lines of Horn and Johnson, Matrix Analysis, §7.4.

Section 4.5 of the same paper treats a second, "orthographic" retraction on the Stiefel manifold, which corrects a tangent step by a normal vector instead of projecting; for the orthogonal group $\mathbf O_n$ (the case $m=n$) it admits a closed form through a matrix square root (Proposition 4.12).

Setting

Fix natural numbers $1\le m\le n$. Matrices are real, and $\mathbb R^{n\times m}$ carries the Frobenius norm

$$\|X\|^2=\sum_{i,j}X_{ij}^2=\operatorname{trace}(X^\top X).$$

The Stiefel manifold is

$$V_{n,m}=\{X\in\mathbb R^{n\times m}:\ X^\top X=I_m\},$$

the set of matrices with orthonormal columns; $V_{n,n}$ is the orthogonal group $\mathbf O_n$.

The singular values of $X$ are $\sigma_1(X)\ge\sigma_2(X)\ge\dots\ge\sigma_{\min\{n,m\}}(X)\ge0$, the square roots of the eigenvalues of $X^\top X$. A singular value decomposition of $X$ is a factorization $X=U\Sigma V^\top$ with $U=[u_1,\dots,u_n]\in\mathbf O_n$, $V=[v_1,\dots,v_m]\in\mathbf O_m$, and $\Sigma\in\mathbb R^{n\times m}$ zero off its diagonal, with nonnegative nonincreasing diagonal entries. For $\bar X\in V_{n,m}$ every singular value equals $1$; in particular $\sigma_m(\bar X)=1$.

A projection of $X$ onto $V_{n,m}$ is a point $Y\in V_{n,m}$ with $\|X-Y\|\le\|X-Z\|$ for all $Z\in V_{n,m}$. A polar decomposition of $X$ is a factorization $X=WS$ with $W\in V_{n,m}$ and $S\in\mathbb R^{m\times m}$ symmetric positive definite.

Formalization targets

Goal: Proposition 3.4 (p. 10)

Let $\bar X\in V_{n,m}$ and let $X$ satisfy $\|X-\bar X\|<\sigma_m(\bar X)$. For every singular value decomposition $X=U\Sigma V^\top$,

$$\{\,Y:\ Y\text{ is a projection of }X\text{ onto }V_{n,m}\,\}=\Big\{\sum_{i=1}^m u_iv_i^\top\Big\},$$

and $\sum_{i=1}^m u_iv_i^\top$ is the $W$ of the polar decomposition $X=WS$: it is the orthonormal factor of every polar decomposition of $X$, and a polar decomposition with this factor exists.

The statement asserts existence, uniqueness and the closed form of the projection on the whole open ball of radius $\sigma_m(\bar X)$, for whichever singular value decomposition is supplied.

Steps of the proof (milestones)

1. For $Y\in V_{n,m}$: $\|X-Y\|^2=\|X\|^2+m-2\operatorname{trace}(Y^\top X)$.
2. For every $Y\in V_{n,m}$, $\operatorname{trace}(Y^\top X)\le\sum_{i=1}^m\sigma_i$, with equality at $Y=\sum_{i=1}^m u_iv_i^\top$.
3. If $\|X-\bar X\|<\sigma_m(\bar X)$ with $\bar X\in V_{n,m}$, then $X$ has full rank $m$.
4. The polar factor of a full-rank matrix is unique (Horn and Johnson, Theorem 7.3.2).

Further items (§4.5)

$$(S+I)^2=I-\Omega^\top\Omega \tag{4.14}$$

for $X\in\mathbf O_n$, $\Omega$ skew-symmetric and $S$ symmetric with $X+X\Omega+XS\in\mathbf O_n$; and Proposition 4.12,

$$R(X,X\Omega)=X\big(\Omega+\sqrt{I-\Omega^\top\Omega}\big),$$

stated as: $S_+=-I+\sqrt{I-\Omega^\top\Omega}$ is the unique symmetric correction of smallest Frobenius norm.

Significance

Proposition 3.4 makes the projective retraction on the Stiefel manifold computable from one singular value decomposition of an $n\times m$ matrix, and identifies it with the polar factor, which is also what many numerical codes already compute for re-orthonormalization. Combined with Proposition 3.2 of the paper, it yields a second-order-accurate retraction usable in any Riemannian algorithm on $V_{n,m}$. The same statement is the orthogonal Procrustes problem in the special case of a full-rank target: the nearest orthonormal frame to $X$.

The result is classical and proved; the paper's argument is short but cites two external facts (Weyl's perturbation bound for singular values and the uniqueness of the polar decomposition). To our knowledge no machine-checked proof of the nearest-orthonormal-frame property, of the uniqueness of the polar factor, or of the closed form of Proposition 4.12 exists in Mathlib or on this platform. A formalization produces reusable pieces: the trace inequality over the Stiefel manifold, the uniqueness of the polar decomposition, and the full-rank property near $V_{n,m}$, all of which recur in matrix analysis and in the analysis of Riemannian algorithms.

Difficulty

Existence of a nearest point is easy (the Stiefel manifold is compact), and attainment of the trace bound is a direct computation. The content is in two places. First, the bound $\operatorname{trace}(Y^\top X)\le\sum_i\sigma_i$ for every $Y\in V_{n,m}$ requires transporting $Y$ by the orthogonal factors of the singular value decomposition and bounding diagonal entries of a matrix with orthonormal columns. Second, uniqueness does not follow from the trace argument: when $X$ is rank deficient there are many maximizers. Uniqueness needs both a perturbation bound (singular values are 1-Lipschitz in the Frobenius norm, so $\sigma_m(X)>0$ near $\bar X$) and the uniqueness of the polar factor, which in turn rests on the uniqueness of the positive-semidefinite square root of $X^\top X$. Mathlib has singular values of linear maps and the spectral theorem but, at the pinned revision, neither Weyl's inequality for singular values nor the polar decomposition.

For Proposition 4.12, the printed proof compares only two solutions $S_\pm$ of (4.14), while (4.14) has other symmetric solutions (mixed signs of the square roots, or non-diagonal ones on repeated eigenvalues); minimality must be proved against all of them.

Formalization scope

Matrices are Matrix (Fin n) (Fin m) ℝ with the Frobenius norm brought in by open scoped Matrix.Norms.Frobenius. The Stiefel manifold is the set stiefel n m of matrices with Xᵀ * X = 1. Singular values are Mathlib's LinearMap.singularValues of Matrix.toEuclideanLin X, re-indexed to be 1-based as on the page (sv X m is $\sigma_m(X)$). A singular value decomposition is the predicate IsSVD X U S V of display (3.5), and $\sum_{i=1}^m u_iv_i^\top$ is U * E * Vᵀ with E the $n\times m$ rectangular identity (frameOfSVD U V). The projection is the set of nearest points, using the published predicate RandomGradFree.Nonsmooth.IsMetricProjection; "exists and is unique" is equality of that set with a singleton. Positive definiteness is Mathlib's Matrix.PosDef, which over $\mathbb R$ includes symmetry.

Standing assumptions and added hypotheses: $m\le n$ is the page's assumption of §3.3; $0<m$ is added so that $\sigma_m$ is meaningful. The radius is $\sigma_m(\bar X)$ with strict inequality, kept in that form although its value is $1$. Two misprints of the page are corrected in the Lean and kept in the verbatim milestone texts: the display of Proposition 3.4 reads $P_{\mathcal R_r}(X)$ for $P_{V_{n,m}}(X)$, and the distance identity reads $m^2$ for $m$. The uniqueness of the polar factor is stated for positive semidefinite factors under the rank hypothesis, which contains the positive definite case of Proposition 3.4. The matrix square root of §4.5 is defined through Mathlib's spectral theorem; Proposition 4.12 is stated for every admissible symmetric correction, and is vacuous only when no symmetric correction exists.

A trivializing formalization is ruled out: the projection is not taken as a hypothesis or chosen by definition; the goal asserts that the nearest-point set equals the singleton of the explicit matrix, for every singular value decomposition supplied, so neither existence nor uniqueness can be assumed away.

Contributions welcome beyond the stated items: Weyl's inequality $|\sigma_i(X)-\sigma_i(Y)|\le\|X-Y\|$ in Mathlib's singular-value API, existence of a singular value decomposition in the matrix form (3.5), and the polar decomposition with its uniqueness, all reusable well beyond this mission.

Selected references

• P.-A. Absil and J. Malick, Projection-like retractions on matrix manifolds, SIAM J. Optim. 22(1):135–158, 2012. https://doi.org/10.1137/100802529 ; authors' version https://hal.science/hal-00651608v2
• P.-A. Absil, R. Mahony and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244
• R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, 1985 (cited by the paper in its 1989 printing) (Theorem 7.3.2, §7.4). https://doi.org/10.1017/CBO9780511810817
• N. J. Higham, Matrix nearness problems and applications, in Applications of Matrix Theory (M. J. C. Gover and S. Barnett, eds.), Oxford University Press, 1989, pp. 1–27 (cited by the paper as [15, §4]; no DOI).
• A. Edelman, T. A. Arias and S. T. Smith, The geometry of algorithms with orthogonality constraints, SIAM J. Matrix Anal. Appl. 20(2):303–353, 1998. https://doi.org/10.1137/S0895479895290954

Motivation

An orthogonal projection onto a closed subspace is straightforward to describe when that subspace is given directly. It is less straightforward when the subspace is specified as the closed sum of two others: a vector can have contributions from both, and the two component projections generally do not commute. In §12 of Aronszajn's 1950 paper, this problem arises while expressing the reproducing kernel of a sum of two closed subspaces of a reproducing-kernel Hilbert space through their individual kernels. The projection formula is the analytic heart of that calculation. It also gives a series whose finite partial sums can be applied without first describing a basis for the closed sum.

Aronszajn states the result for closed subspaces of a complex Hilbert space of functions. The projection identity itself uses only Hilbert-space geometry; evaluation at a point and the reproducing property enter when the operator formula is translated into a kernel formula later in the section. This mission isolates the projection theorem so that the same formal result can be used both in that kernel calculation and in other settings with two closed subspaces.

Setting

Let $E$ be a complex Hilbert space and let $F_1,F_2$ be closed linear subspaces. The algebraic sum $F_1+F_2$ is the set of all $f_1+f_2$ with $f_i\in F_i$. It need not be closed. Write $F'=\overline{F_1+F_2}$ for its closure and $F_0=F_1\cap F_2$ for the intersection. All four subspaces are closed except possibly the algebraic sum itself. Let $P_1,P_2,P,P_0$ be the orthogonal projections onto $F_1,F_2,F',F_0$, respectively. Operator products denote composition, with the rightmost factor applied first; $(P_2P_1)^0$ is the identity.

The paper writes $F_1\oplus F_2$ for $F'$, even when $F_0$ is nonzero, and uses a distinct dotted plus for the algebraic sum. These symbols can look like a direct sum in newer notation. Here the definitions of $F'$ and $F_0$ remove that ambiguity. Aronszajn assumes complex scalars throughout Part I after §1, p. 343; the statements below keep that convention. Neither a topology on an underlying set of functions nor a measure is part of these operator assertions.

Formalization targets

The goal is the strong projection series, §12, Eq. (7). For each $f\in E$,

$$Pf=P_0f+\sum_{k=1}^{\infty}\left[P_1(P_2P_1)^{k-1}+P_2(P_1P_2)^{k-1}-(P_2P_1)^k-(P_1P_2)^k\right]f.$$

The sum means that its finite partial sums converge in the norm of $E$ for each fixed $f$. No operator-norm limit is claimed. The first term of the sum, at $k=1$, includes $P_1+P_2-P_2P_1-P_1P_2$ because both zero powers are identity operators.

Three source statements are milestones. Equation (1) is a finite identity that expresses a power of $(P-P_1)(P-P_2)$ through a partial sum and a power of $P_2P_1$. The paragraph following Eq. (4) identifies the strong limit of $(P_2P_1)^m$ as $P_0$. The paragraph before Eq. (7) asserts that $[(P-P_1)(P-P_2)]^m$ tends strongly to zero. Together these statements specify the finite and limiting parts of the displayed target. They require no assumption that $F_0$ is zero.

Significance

The formula describes the orthogonal projection onto a closed sum in terms of projections onto its two constituents. In Aronszajn's application, an orthogonal projection on a closed subspace determines that subspace's reproducing kernel; therefore the operator expansion supplies a kernel expansion without constructing a basis of the sum. It also makes explicit why an intersection term survives when the subspaces overlap. The result is proved in the 1950 paper; this mission asks for a Lean proof of that known result, with its convergence mode and closure convention made explicit.

A machine-checked development would provide a reusable complex-Hilbert-space statement for alternating projections and closed sums. Mathlib already provides closed submodules and their orthogonal projections, so the new contribution is the finite identity and the strong-limit claims connecting them. The resulting theorem can support the paper's later kernel expression, §12, Eq. (18), once the correspondence between bounded operators and kernels from §11 is available. Equation (18) is outside this mission's target list.

Difficulty

The algebraic sum of two closed subspaces can fail to be closed, so projecting onto it directly would leave the target undefined as a Hilbert-space orthogonal projection. Passing to its closure gives the correct target but does not imply the two projections commute or that their product has a uniform contraction factor. In particular, strong convergence of the powers and of the final series does not generally improve to convergence in operator norm. The later special case in §12, where the intersection is zero and the minimal angle is positive, has a stronger uniform-convergence conclusion; that extra angle condition is absent from Eq. (7).

The finite formula also has four terms per summand whose order matters. Reversing $P_1P_2$ and $P_2P_1$, beginning the power at $k$ rather than $k-1$, or omitting $P_0$ changes the statement. Even a proof of pointwise convergence of each selected term would not by itself establish the convergence of the bracketed partial sums in the target.

Formalization scope

The Lean representation is an arbitrary complex Hilbert space E with two ClosedSubmodule ℂ E objects. Mathlib's Submodule.starProjection supplies $P_1$ and $P_2$. Two small definitions name $P$ as projection onto the join of the closed submodules and $P_0$ as projection onto their meet. The closed-submodule join is the closure of the algebraic sum; the meet is the intersection. These definitions are computed from $F_1,F_2$, not independently quantified projections. The theorem quantifies over every vector $f$ and uses Filter.Tendsto along natural-number partial sums, giving norm convergence of the resulting vectors.

The source's ambient class has a reproducing kernel, but no step in Eqs. (1)–(7) uses point evaluations. The formal statements therefore apply to any complex Hilbert space, including the paper's RKHS. The stronger scope is explicit here and does not assert a new kernel identity. The $m$-th Lean partial sum uses indices $0,\ldots,m-1$ for the paper's $1,\ldots,m$; it has the same terms and is defined at $m=0$ as $P_0f$. The finite identity is stated only for $m\ge1$, as its $m=0$ version is false in general.

Eq. (1) on p. 375 is printed without the opening bracket of the summand (only the closing bracket appears); the formalization reads it as in Eq. (4) on p. 376, where both brackets are printed. There are two printing slips on p. 377: a running sentence drops brackets from the power of $(P-P_1)(P-P_2)$, and it prints $P\ominus P_2$ while discussing the projection $P-P_2$. The formalization follows the bracketed expression in Eqs. (1) and (4), using ordinary operator subtraction. The corresponding milestone preserves the printed wording. No assumption from the later angle analysis, especially $F_0=\{0\}$, is imported into these results. The series is stated through actual partial sums; it is not defined by an arbitrary RKHS chosen to have a desired kernel.

Selected references

• N. Aronszajn, Theory of Reproducing Kernels, Transactions of the American Mathematical Society 68 (1950), 337–404, §12, pp. 375–380. DOI: 10.1090/S0002-9947-1950-0051437-7.

Motivation

Most Markovian queueing models have no closed-form transient solution. The M/M/1 queue already needs modified Bessel functions (Chapter 2 of the book), and a finite-capacity or multi-class model with state-dependent rates has no closed form at all. What an analyst can always write down is the system of forward equations $p'(t)=p(t)Q$ for the state probabilities. Chapter 8 of Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory (4th ed., Wiley 2008, DOI 10.1002/9781118625651), presents two numerical techniques that turn such models into numbers: the randomization (or uniformization) method for the transient distribution of a finite continuous-time Markov chain, and the Fourier-series method for inverting a Laplace transform, as needed for the M/G/1 waiting-time transform (5.33) and the busy-period transform (5.37).

Uniformization goes back to Jensen (1953) and is the standard transient solver in performance-evaluation and reliability tools. Its appeal is that it replaces a matrix exponential, which is numerically delicate, by powers of a stochastic matrix weighted by Poisson probabilities, with an error bound that can be fixed before the computation starts (Grassmann 1977; Gross and Miller 1984). The Fourier-series method with Euler summation is due to Abate and Whitt (Abate and Whitt 1992; Abate, Choudhury and Whitt 1999).

Setting

A continuous-time Markov chain $X(t)$ on the states $\{0,1,\dots,N\}$ is described by its infinitesimal generator $Q=(q_{ij})$: for $i\ne j$, $q_{ij}\ge0$ is the rate of jumps from $i$ to $j$, and the diagonal entry is $-q_i$ with

$$q_i=\sum_{j\ne i}q_{ij},\qquad i=0,1,\dots,N.$$

The transient state-probability vector $p(t)=(p_0(t),\dots,p_N(t))$, $p_n(t)=\Pr\{X(t)=n\}$, is the solution of the forward equations

$$p'(t)=p(t)Q\quad(t\ge0),$$

started from a given probability vector $p(0)$. Fix a constant $\Lambda>0$ with $\Lambda\ge q_i$ for every $i$ (the book takes $\Lambda=\max_i q_i$) and define the uniformized matrix

$$\tilde P=\frac{Q}{\Lambda}+I,\qquad \tilde p_{in}=\begin{cases}q_{in}/\Lambda&(i\ne n),\\1-q_i/\Lambda&(i=n).\end{cases}$$

It is the transition matrix of a discrete-time chain $Y_k$: the state of $X$ after the $k$-th event of a Poisson process of rate $\Lambda$ that has been thinned. Write $\phi^{(k)}=p(0)\tilde P^{k}$ for its distribution after $k$ steps.

For the second half of the chapter, the Laplace transform of a real function $f$ on $[0,\infty)$ is $\bar f(s)=\int_0^\infty e^{-st}f(t)\,dt$, and the Fourier-series approximant with parameter $A$ is

$$f_{A,n}(t)=\frac{e^{A/2}}{2t}\Big[\bar f\Big(\frac{A}{2t}\Big)+2\sum_{k=1}^{n}(-1)^k\,\mathrm{Re}\,\bar f\Big(\frac{A+2k\pi i}{2t}\Big)\Big],$$

with $f_A(t)=\lim_{n\to\infty}f_{A,n}(t)$.

Formalization targets

Goal: the randomization formula with its truncation bound (Eqs. (8.9)–(8.12))

The forward equations have a solution, and every solution satisfies, for all $t\ge0$,

$$p(t)=\sum_{k=0}^{\infty}p(0)\tilde P^{(k)}\,\frac{e^{-\Lambda t}(\Lambda t)^k}{k!},$$

and whenever $\sum_{k=0}^{T}e^{-\Lambda t}(\Lambda t)^k/k!>1-\epsilon$, every component of the sum truncated at $k=T$ is within $\epsilon$ of $p_n(t)$.

Milestones

1. Eq. (8.12): $\tilde P$ has the entries above and is a stochastic matrix.
2. Eqs. (8.13)–(8.14): $\phi^{(k)}=\phi^{(k-1)}\tilde P$ and each $\phi^{(k)}$ is a probability vector.
3. p.385: $\phi=\phi\tilde P\iff0=\phi Q$.
4. Eqs. (8.27)–(8.28): for bounded Lipschitz $f$, $A>0$ and $t>0$,

$$f_A(t)-f(t)=\sum_{k=1}^{\infty}e^{-kA}f\big((2k+1)t\big),\qquad |f_A(t)-f(t)|\le\frac{Ce^{-A}}{1-e^{-A}}\ \text{ if } |f(x)|\le C \text{ for } x>3t.$$

The mission also contains Eqs. (8.7)–(8.8) as a further theorem, outside the milestone list: the transition probabilities satisfy $p_{in}(t)=\sum_k\tilde p^{(k)}_{in}e^{-\Lambda t}(\Lambda t)^k/k!$, and $p_n(t)=\sum_i p_i(0)p_{in}(t)$.

Significance

The randomization formula reduces the transient analysis of any finite Markovian queue (finite-buffer, multi-server, with balking, reneging or state-dependent rates) to repeated vector–matrix products with a sparse stochastic matrix. The truncation point is chosen from a Poisson tail alone, independently of $Q$. Milestone 3 shows that the same matrix gives the stationary equations, so one iteration serves both transient and steady-state computation. The discretization identity (8.27) is what justifies the parameter choice in Algorithm 8.1: the error decays like $e^{-A}$.

All of these results are classical and proved in the literature. None of them is formalized in Lean or Mathlib as far as a search of the platform and Mathlib shows. Mathlib has the matrix exponential and Poisson summation under decay hypotheses, but no continuous-time Markov chain generators, no uniformization, and no Laplace transform. This mission would add the finite-state link between generators, stochastic matrices and matrix exponentials that later chapters of queueing and reliability theory use, and a verified error formula for a numerical inversion method in wide use.

Difficulty

The book's derivation is probabilistic: it conditions on the number of events of the Poisson($\Lambda$) process and thins them. A formal statement cannot rest on that picture, because $p(t)$ is defined analytically, by the forward equations. The goal therefore contains a uniqueness statement for a linear ODE on $[0,\infty)$ with one-sided derivative at $0$, which the book never mentions. The componentwise bound then needs $\tilde P$ to be stochastic, so that every $\phi^{(k)}_n$ lies in $[0,1]$. That is exactly where $\Lambda\ge\max_i q_i$ is used; with a smaller $\Lambda$ the matrix $\tilde P$ has negative diagonal entries and the bound fails.

For (8.27), the book gives no proof. The identity is an aliasing (Poisson-summation) formula for a periodic function assembled from the values of $f$ at all odd multiples of $t$. The convergence of the conditionally summed series (8.24) is the delicate point: continuity of $f$ at $t$, the book's only hypothesis, does not guarantee convergence of a Fourier series. Mathlib's Poisson summation theorems require decay of the Fourier transform that the damped, reflected function built from $f$ does not have.

Formalization scope

• States are Fin (N+1); a row vector is Fin (N+1) → ℝ; $pQ$ is vecMul. A generator is a real matrix with nonnegative off-diagonal entries and diagonal $-\sum_{j\ne i}q_{ij}$.
• $p(t)$ is not defined as the series. It is any function with $p(0)=p_0$ and one-sided derivative $p(t)Q$ within $[0,\infty)$ at every $t\ge0$. The goal also asserts that such a function exists, so it cannot hold vacuously, and it asserts the series identity for every solution. Defining $p(t)$ as the series (8.9) would make the goal a tautology and is ruled out.
• $\Lambda$ is any real with $\Lambda>0$ and $\Lambda\ge q_i$ for all $i$ (the book takes equality with $\max_i q_i$).
• The truncation bound is stated componentwise, as on p.384 ("an error bound on $p_n(t)$ of $\epsilon$"), for an arbitrary real $\epsilon$ and truncation point $T$.
• The series (8.8), (8.9) are stated with HasSum, so convergence is part of the claim.
• The Laplace transform is the Lebesgue integral over $(0,\infty)$ at a complex argument. $f_A(t)$ is the limit of the partial sums $f_{A,n}(t)$, and the convergence is part of milestone 4.
• Strengthened hypotheses in milestone 4: $f$ bounded and Lipschitz on $[0,\infty)$ replaces "$t$ is a continuity point of $f$", which is not sufficient for convergence.
• Corrected misprints: $e^{-\lambda t}$ in (8.9) is $e^{-\Lambda t}$; $q_{ij}/\Lambda$ in (8.12) is $q_{in}/\Lambda$; $\phi(Q/\Lambda-I)$ on p.385 is $\phi(Q/\Lambda+I)$.
• Not formalized: Theorem 8.1 (Bromwich inversion) and the real form (8.21), which the book states without hypotheses on $f$; the limit claim $\lim_k\phi^{(k)}=\lim_t p(t)$ on p.385, which fails when $\tilde P$ is periodic; the Euler-summation approximation (8.26) and the round-off discussion, which are stated with "≈".

Useful infrastructure: the matrix exponential and its derivative (Matrix, NormedSpace.exp), uniqueness for linear ODEs (Grönwall), Fourier series on the circle, and a reusable Laplace transform file. Contributions of general lemmas on generators and stochastic matrices are welcome, as they apply to every finite Markovian model in the series.

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008, §§8.1.2–8.2. https://doi.org/10.1002/9781118625651
• A. Jensen, "Markoff chains as an aid in the study of Markoff processes", Skandinavisk Aktuarietidskrift 36 (1953) 87–91.
• W. K. Grassmann, "Transient solutions in Markovian queueing systems", Computers & Operations Research 4 (1977) 47–53.
• D. Gross, D. R. Miller, "The randomization technique as a modeling tool and solution procedure for transient Markov processes", Operations Research 32 (1984) 343–361. https://doi.org/10.1287/opre.32.2.343
• J. Abate, W. Whitt, "The Fourier-series method for inverting transforms of probability distributions", Queueing Systems 10 (1992) 5–87. https://doi.org/10.1007/BF01158520
• J. Abate, G. L. Choudhury, W. Whitt, "An introduction to numerical transform inversion and its application to probability models", in W. Grassmann (ed.), Computational Probability, Kluwer, 1999, 257–323.