All missions

Motivation

In the local model of differential privacy, no trusted curator holds the data: each individual randomizes their own record before handing it to the analyst. This is the model of randomized response in survey statistics and of the privacy-preserving telemetry deployed by large software vendors. Kasiviswanathan, Lee, Nissim, Raskhodnikova and Smith, What Can We Learn Privately? (arXiv:0803.0924v3, published in SIAM J. Comput. 40(3) (2011) 793–826, DOI 10.1137/090756090), asked which learning tasks remain possible in this model. Their answer (§5) is that local algorithms are exactly as powerful as statistical query (SQ) algorithms in the sense of Kearns (J. ACM 1998), up to polynomial factors. This mission formalizes one direction of that equivalence: every local algorithm run on i.i.d. data can be simulated by an SQ algorithm.

This is the fourth of five missions on the paper. Mission III formalizes the converse direction (Theorem 5.7). All theorem numbers refer to arXiv:0803.0924v3.

Setting

A database is $z=(z_1,\dots,z_n)\in D^n$. Here its entries are drawn i.i.d. from a probability distribution $P$ on $D$.

• An $\varepsilon$-local randomizer (Definition 5.1) is a randomized map $R:D\to W$ to a discrete set $W$ such that $\Pr[R(u)=w]\le e^{\varepsilon}\Pr[R(u')=w]$ for all $u,u'\in D$ and $w\in W$.
• An $\varepsilon$-local algorithm $A$ making $t$ queries (Definitions 5.2, 5.3) sees the database only through an LR oracle. At its $k$-th call it names an index $i_k$ and an $\varepsilon_k$-local randomizer $R_k$, both possibly depending on its earlier answers, and receives a fresh sample of $R_k(z_{i_k})$. For each index $i$ the budgets of the calls on $i$ sum to at most $\varepsilon$. $A$ is noninteractive if its calls do not depend on earlier answers. Its output is a function of the $t$ answers. Its output distribution is taken over $z\sim P^n$ and the randomizers' coins.
• An SQ oracle for $P$ (Definition 5.4) answers a query $(g,\tau)$, with $g:D\to[-1,1]$ and tolerance $\tau$, with any number $v$ such that $|v-\mathbb E_{u\sim P}[g(u)]|\le\tau$. It may choose its answers adversarially and adaptively.
• An SQ algorithm $B$ (Definition 5.5) accesses $P$ only through such an oracle. It may use its own coins, and the number of queries it makes may be random and unbounded.
• The statistical difference of two distributions on a discrete space (p. 8) is $\max_S|\mu(S)-\nu(S)|$.

Formalization targets

Goal: Lemma 5.8 (p. 21)

There is an absolute constant $C$ such that for every $\varepsilon$-local algorithm $A$ making $t$ queries there is an SQ algorithm $B$ with the following properties against every distribution $P$ and every valid oracle:

$$\tau=\frac{\beta}{3e^{2\varepsilon}t},\qquad \mathbb E[\#\text{queries of }B]\le C\,t\,e^{\varepsilon},\qquad \mathrm{SD}\bigl(B,\ A(z),\ z\sim P^n\bigr)\le\beta .$$

The constant $C$ is the paper's $O(\cdot)$. The tolerance is the proof's own choice.

Milestones

1. Display (5) (p. 22): a single query estimates $p(w)=\Pr_{z_i\sim P}[R(z_i)=w]$ within a factor $1\pm\beta/(3t)$, whatever valid answer the oracle gives.
2. One rejection-sampling iteration (p. 23): every iteration terminates with probability at least $\frac{1-\varphi}{1+\varphi}e^{-\varepsilon}$. Conditioned on terminating, it outputs $w$ with probability in $(1\pm3\varphi)p(w)$.
3. The interactive estimate (p. 24): two queries estimate the conditional probability of the next answer given earlier answers on the same entry, within a factor $1\pm3e^{2\varepsilon}\tau$.
4. Claim 5.9 (p. 22): the noninteractive case, with the explicit bound $2t\,e^{\varepsilon}$ on the expected number of queries.

Significance

The result. Lemma 5.8, combined with Theorem 5.7, shows that a concept class is learnable by a locally private algorithm if and only if it is learnable with statistical queries (Theorem 5.14). Every SQ lower bound therefore becomes a lower bound for local privacy. For instance, parity functions are privately learnable by a centralized algorithm (Theorem 4.4) but not by a local one, because parities need exponentially many statistical queries (Corollary 5.15, which also uses the SQ lower bound of Blum et al., STOC 1994).

Formalizing it. The result is proved on paper. As far as we know it has no machine-checked proof. The paper's proof is a few paragraphs and leaves the model implicit: what an SQ algorithm with a random number of queries is, how an adversarial oracle interacts with it, and in what sense the per-randomizer errors add up over $t$ adaptively chosen randomizers. The formal statement makes each of these explicit. The milestones isolate the estimates the proof uses, so they can be proved independently of the composition argument.

Difficulty

The estimates in the milestones are elementary inequalities. The difficulty is in the goal. First, the oracle's answers, and therefore the estimates $\tilde p(w)$, change from iteration to iteration, so the simulated distribution is not a fixed rejection sampler: the output law must be controlled iteration by iteration against an adversary. Second, the SQ algorithm has no bound on its number of steps. Its output law is a limit over an unbounded run, and its query count is an expectation that is finite only because every iteration terminates with probability bounded away from zero. Third, in the interactive case the randomizers applied to one entry are correlated through the entry. The simulation must sample from the conditional law given earlier answers on that entry, while the errors of all $t$ steps combine along adaptively chosen histories.

Formalization scope

• Local randomizers have discrete output: R u : PMF W. The simulation uses the point probabilities $\Pr[R(z_i)=w]$. Every map $u\mapsto\Pr[R(u)=w]$ is required to be measurable.
• A local algorithm makes exactly $t$ calls. Its index, randomizer and budget at call $k$ are functions of the earlier answers. The budget is required along every answer sequence. Its own coins are calls to $0$-local randomizers.
• An SQ algorithm is a state machine with a random start, random transitions depending on the answer, and an output on stopping. Its output law is a sub-probability mass function. Its expected query count lies in $[0,\infty]$ and counts every query, including those in rejected iterations.
• The oracle is an arbitrary deterministic function of the whole state trajectory, constrained only by Definition 5.4. The conclusions must hold for every valid oracle, not for the exact-mean oracle alone.
• $B$ is quantified before $P$ and sees $P$ only through answers, and every query of $B$ is a measurable $[-1,1]$-valued function with tolerance exactly $\beta/(3e^{2\varepsilon}t)$. Without these two constraints the statement would be trivial: $B$ could hard-code $P$, or rescale a query to shrink its effective tolerance.
• The constant $C$ is quantified outside every other object.
• The implicit hypotheses $\varepsilon>0$ (the queries divide by $e^{\varepsilon}-e^{-\varepsilon}$), $0<\beta\le1$ and $t\ge1$ (so that $\varphi=\beta/(3t)\le1/3$ and $\tau$ is defined) are stated. The paper's reference input $\mathbf 0$ is an arbitrary point $u_0\in D$.
• Claim 5.9 states "$t\cdot e^{\varepsilon}$ queries". Its proof gives at most $2t\,e^{\varepsilon}$, which is the bound stated.
• The paper's refinement "noninteractive $A$ yields nonadaptive $B$" is not formalized. The simulation decides from each answer whether to stop, and so which randomizer the next query belongs to. It therefore does not prepare its queries before receiving answers in the sense of Definition 5.5. The goal is the adaptive statement for all local algorithms, which contains the noninteractive case. Claim 5.10, whose statement coincides with this goal, is represented by its estimate (milestone 3) rather than restated.
• The statistical difference is valued in $[0,\infty]$, which avoids junk values of a real supremum.

Reusable beyond this mission: the SQ-algorithm state machine with expected query count, the statistical difference of sub-probability mass functions, and the model of interactive local algorithms. Welcome contributions include proofs of the estimate milestones, a general lemma bounding the statistical difference of sequentially composed approximate samplers, and the termination and expected-runtime analysis of rejection sampling with varying acceptance probabilities.

Selected references

• S. P. Kasiviswanathan, H. K. Lee, K. Nissim, S. Raskhodnikova, A. Smith, What Can We Learn Privately?, SIAM J. Comput. 40(3) (2011) 793–826; cited from arXiv:0803.0924v3. https://arxiv.org/abs/0803.0924 · https://doi.org/10.1137/090756090
• M. Kearns, Efficient noise-tolerant learning from statistical queries, J. ACM 45(6) (1998) 983–1006. https://doi.org/10.1145/293347.293351
• A. Blum, M. Furst, J. Jackson, M. Kearns, Y. Mansour, S. Rudich, Weakly learning DNF and characterizing statistical query learning using Fourier analysis, STOC 1994, 253–262. https://doi.org/10.1145/195058.195147
• C. Dwork, F. McSherry, K. Nissim, A. Smith, Calibrating noise to sensitivity in private data analysis, TCC 2006. https://doi.org/10.1007/11681878_14
• S. L. Warner, Randomized response: a survey technique for eliminating evasive answer bias, J. Amer. Statist. Assoc. 60 (1965) 63–69. https://doi.org/10.1080/01621459.1965.10480775

Motivation

A backward stochastic differential equation (BSDE) prescribes the value of a process at a terminal time and asks for an adapted process that reaches it while following a given drift. In mathematical finance the price of a contingent claim, and its hedging strategy, solve such an equation; when the market has frictions (different borrowing and lending rates, for example) the drift, called the driver, is nonlinear and no closed form exists. Numerical methods for BSDEs are therefore methods for pricing and hedging under nonlinear models, and also for semilinear parabolic PDEs, which BSDEs represent probabilistically.

Gobet, Lemor and Warin (Ann. Appl. Probab. 15 (2005), arXiv:math/0508491) proposed and analysed a simulation scheme in which every conditional expectation of a backward time-stepping recursion is replaced by a least-squares regression on finitely many functions, as in the Longstaff–Schwartz method for American options. Their analysis splits the total error into three parts: time discretization (Theorem 1, from Zhang's results), replacing conditional expectations by $\mathbf L_2$ projections on function bases (Theorem 2), and replacing those projections by empirical regressions on $M$ simulated paths (Theorem 3). This mission formalizes the second part.

Timeline. Zhang (Ann. Appl. Probab. 14 (2004)) and Bouchard and Touzi (Stoch. Proc. Appl. 111 (2004)) established the $\sqrt h$ rate of the time discretization. Bouchard and Touzi's regression error (their reference [6] in the paper, Theorem 4.1 there) was expressed through the residuals of the scheme's own iterates. Gobet, Lemor and Warin (2005) gave the bound in terms of the residuals of the discrete BSDE, together with estimates on $Z$.

Setting

Fix a horizon $T>0$, dimensions $d,q\ge1$, a drift $b(t,x)\in\mathbb R^d$ and a diffusion matrix $\sigma(t,x)\in\mathbb R^{d\times q}$, both Lipschitz in $(t,x)$ ((H1)), and a driver $f(t,x,y,z)\in\mathbb R$ with

$$|f(t_2,x_2,y_2,z_2)-f(t_1,x_1,y_1,z_1)|\le C_f\big(|t_2-t_1|^{1/2}+|x_2-x_1|+|y_2-y_1|+|z_2-z_1|\big)\qquad\textbf{(H2)}.$$

For $N\ge1$ put $h=T/N$ and $t_k=kh$. On a probability space with a filtration $(\mathcal F_k)$, the increments $\Delta W_k\in\mathbb R^q$ are $\mathcal F_{k+1}$-measurable, independent of $\mathcal F_k$ and Gaussian $\mathcal N(0,hI_q)$; $\Delta W_{l,k}$ is the $l$-th component. The Euler scheme is $S^N_{t_0}=S_0$, $S^N_{t_{k+1}}=S^N_{t_k}+b(t_k,S^N_{t_k})h+\sigma(t_k,S^N_{t_k})\Delta W_k$. An $\mathcal F_k$-adapted process $P^N_{t_k}\in\mathbb R^{d'}$ extends $S^N_{t_k}$ by extra state variables, and the terminal value is $\Phi^N(P^N_{t_N})$, square integrable.

Write $\mathbb E_k=\mathbb E(\cdot\mid\mathcal F_k)$. The discrete BSDE is $Y^N_{t_N}=\Phi^N(P^N_{t_N})$ and, for $k<N$,

$$Z^N_{l,t_k}=\tfrac1h\mathbb E_k(Y^N_{t_{k+1}}\Delta W_{l,k}),\qquad Y^N_{t_k}=\mathbb E_k(Y^N_{t_{k+1}})+hf(t_k,S^N_{t_k},Y^N_{t_k},Z^N_{t_k}).$$

A function basis $p_{l,k}(P^N_{t_k})\in\mathbb R^{n_{l,k}}$ ($0\le l\le q$) is square integrable with invertible Gram matrix $\mathbb E(p_{l,k}p_{l,k}^*)$. $\mathcal P_{p}(U)$ is the $\mathbf L_2(\Omega,\mathbb P)$ orthogonal projection of $U$ onto the span of the basis, and $\mathcal R_p(U)=U-\mathcal P_p(U)$.

The projection–Picard scheme with $I$ iterations (Definition 1) produces $Y^{N,i,I}_{t_k}=\alpha^{i,I}_{0,k}\cdot p_{0,k}$ and $Z^{N,i,I}_{l,t_k}=\alpha^{i,I}_{l,k}\cdot p_{l,k}$, starting from $\alpha^{0,I}=0$, where $\alpha^{i,I}_k$ minimizes

$$\mathbb E\Big(Y^{N,I,I}_{t_{k+1}}-\alpha_0\cdot p_{0,k}+hf(t_k,S^N_{t_k},Y^{N,i-1,I}_{t_k},Z^{N,i-1,I}_{t_k})-\sum_{l=1}^q\alpha_l\cdot p_{l,k}\Delta W_{l,k}\Big)^2.\qquad(9)$$

Finally $\mathcal A^N(S_0)=1+|S_0|^2+\mathbb E|\Phi^N(P^N_{t_N})|^2$.

Formalization targets

Goal: Theorem 2 (p. 11)

For $h$ small enough,

$$\max_{0\le k\le N}\mathbb E|Y^{N,I,I}_{t_k}-Y^N_{t_k}|^2+h\sum_{k=0}^{N-1}\mathbb E|Z^{N,I,I}_{t_k}-Z^N_{t_k}|^2\le Ch^{2I-2}\mathcal A^N(S_0)+C\sum_{k=0}^{N-1}\mathbb E|\mathcal R_{p_{0,k}}(Y^N_{t_k})|^2+Ch\sum_{k=0}^{N-1}\sum_{l=1}^q\mathbb E|\mathcal R_{p_{l,k}}(Z^N_{l,t_k})|^2 .$$

Milestones

1. (10)–(11): the minimizer of (9) is given by projections, $Z^{N,i,I}_{l,t_k}=\frac1h\mathcal P_{p_{l,k}}(Y^{N,I,I}_{t_{k+1}}\Delta W_{l,k})$ and $Y^{N,i,I}_{t_k}=\mathcal P_{p_{0,k}}(Y^{N,I,I}_{t_{k+1}}+hf(\dots,Y^{N,i-1,I}_{t_k},Z^{N,i-1,I}_{t_k}))$.
2. (12): $h\,\mathbb E|Z^{N,i,I}_{l,t_k}|^2\le\mathbb E|Y^{N,I,I}_{t_{k+1}}|^2-\mathbb E|\mathbb E_k(Y^{N,I,I}_{t_{k+1}})|^2$.
3. (13): the map $Y\mapsto\mathcal P_{p_{0,k}}(Y^{N,I,I}_{t_{k+1}}+hf(t_k,S^N_{t_k},Y,Z^{N,I,I}_{t_k}))$ is a $(C_fh)$-contraction on $\mathbf L_2(\mathcal F_k)$ with a unique fixed point.
4. The discrete Gronwall lemma with $c$-terms (p. 11, item 3).
5. (19): $\mathbb E|Y^{N,i,I}_{t_k}|^2+h\,\mathbb E|Z^{N,i,I}_{l,t_k}|^2\le C\mathcal A^N(S_0)$, uniformly in $I$, $i$, $k$.

Significance

Theorem 2 shows that the projection errors of a backward regression scheme only add up over the $N$ time steps, with a constant that does not grow with $N$, and that they are measured by the residuals of the discrete BSDE itself. That makes the influence of the basis directly computable (the paper's §6 does so for Voronoi-cell indicators), and shows that $I=2$ Picard iterations already give an error of the order of the time discretization. Combined with Theorem 3 it gives the complete error budget of the algorithm.

The result is proved in the paper; to our knowledge none of it is machine-checked. A complete development provides a formal $\mathbf L_2$-regression calculus for discrete BSDEs (projections on random bases, conditional expectations against Gaussian increments, contraction of Picard maps in $\mathbf L_2(\mathcal F_k)$) and a backward discrete Gronwall lemma, all reusable for other regression schemes. One printed step, (14), fails at $i=1$ (see below), so a formal proof also certifies that the theorem survives the repair.

Difficulty

The obvious argument compares the scheme with the discrete BSDE one step at a time and applies Gronwall. It fails for $Z$: $Z^{N,i,I}_{t_k}$ carries a factor $1/h$, and a naive bound $\mathbb E|Z|^2\le h^{-1}\mathbb E|Y|^2$ summed over $N=T/h$ steps explodes. A usable bound has to account for the conditional variance of $Y^{N,I,I}_{t_{k+1}}$ given $\mathcal F_k$, not only its second moment. The second difficulty is that the projection does not commute with the driver: projection errors enter at every step through the nonlinear $f$, and must be bounded by residuals of $Y^N$ and $Z^N$, not of the scheme's iterates. A third is the Picard step: at $i=1$ the iterate is computed with $Z^{N,0,I}=0$, so it is not an iterate of the contraction of milestone 3, and the printed inequality (14) $\mathbb E|Y^{N,\infty,I}_{t_k}-Y^{N,i,I}_{t_k}|^2\le(C_fh)^{2i}\mathbb E|Y^{N,\infty,I}_{t_k}|^2$ fails there; an extra term in $\mathbb E|Z^{N,I,I}_{t_k}|^2$ is needed.

Formalization scope

The Lean development lives in the namespace RegMCBSDE.Projection. Points are in EuclideanSpace ℝ (Fin d), the matrix norm in (H1) is the Frobenius norm, and all expectations of squares are lower Lebesgue integrals in $[0,\infty]$, so no junk value of a Bochner integral can make an inequality vacuous. Component $m$ (from $0$) of $\Delta W_k$ is the paper's $\Delta W_{m+1,k}$, and the bases are indexed by Fin (q+1) with $l=0$ for $Y$.

Committed readings:

• (H3) is dropped. It constrains the continuous terminal functional, which no statement involves.
• The filtration is abstract. Any filtration with $\mathcal F_{k+1}$-measurable increments independent of $\mathcal F_k$ and of law $\mathcal N(0,hI_q)$; the Brownian filtration is one. The Markov representation of $P^N$ is not used and is dropped.
• Schemes are relations. $(Y^N,Z^N)$ is any solution of (5)–(6), and $\alpha$ is any family satisfying the arg-min rule (9) for every $i\ge1$ (the paper runs $i\le I$; (19) refers to all $i\ge0$). (10)–(11) are a milestone, not the definition.
• The projection is Mathlib's orthogonal projection in $\mathbf L_2(\Omega,\mathbb P)$ onto the span of the basis coordinates. It is not defined by the normal equations.
• Constants. In Theorem 2 and (19), $C$ and the threshold $h_0$ of "$h$ small enough" are chosen after $(T,d,q,b,\sigma,f,C_f,L)$ and before $N$, $I$, $S_0$, $d'$, the probability space, $P^N$, $\Phi^N$ and the bases. A constant chosen after the scheme data would make the theorem trivially true and is ruled out by this quantifier order.
• Pinned readings. $\max_k$ is "for every $k\le N$". In (13) the argument $Z^{N,i-1,I}$ is read as $Z^{N,I,I}$, as the displayed (13) shows, and "$h$ small enough" is $C_fh<1$. (12) is multiplied by $h$ to avoid subtraction. (19) is stated for $k\le N-1$ and every $1\le l\le q$. (10) is stated where $Y^{N,I,I}_{t_{k+1}}\Delta W_{l,k}$ is square integrable, since $\mathcal P$ acts on $\mathbf L_2$.
• Not stated. (14), which is false at $i=1$, and the steps whose printed derivation passes through it ((15)–(18), (20)–(26)); Theorem 1 and Propositions 1 and 3.

Contributions welcome: proofs of the milestones, and the Mathlib-level lemmas they need (conditional expectation of a product with an independent centered Gaussian, $\mathbf L_2$ moments of the Euler scheme, the projection identity $\mathcal P_p(U)=\mathcal P_p(\mathbb E_kU)$ for $\mathcal F_k$-measurable bases).

Selected references

• E. Gobet, J.-P. Lemor, X. Warin, A regression-based Monte Carlo method to solve backward stochastic differential equations, Ann. Appl. Probab. 15(3), 2172–2202, 2005. arXiv:math/0508491, doi:10.1214/105051605000000412
• J. Zhang, A numerical scheme for BSDEs, Ann. Appl. Probab. 14(1), 459–488, 2004. doi:10.1214/aoap/1075828058
• B. Bouchard, N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stoch. Proc. Appl. 111(2), 175–206, 2004. doi:10.1016/j.spa.2004.01.001
• F. A. Longstaff, E. S. Schwartz, Valuing American options by simulation: a simple least-squares approach, Rev. Financ. Stud. 14(1), 113–147, 2001. doi:10.1093/rfs/14.1.113

Motivation

Backward stochastic differential equations (BSDEs) describe the price and the hedge of a contingent claim in models with nonlinear pricing rules (differential interest rates, funding costs, reflected or constrained claims), and give probabilistic representations of semilinear parabolic PDEs (El Karoui, Peng and Quenez, 1997). Their numerical solution in moderate dimension is done by simulation: a backward recursion over a time grid in which each conditional expectation is replaced by a least-squares regression on simulated paths, the same device as the regression method for Bermudan options of Longstaff and Schwartz (2001).

Gobet, Lemor and Warin (2005) split the error of such a scheme into three parts: time discretization (Theorem 1), projection on finite function bases (Theorem 2), and the replacement of $\mathbf L_2$ projections by empirical regressions on $M$ simulated paths (Theorem 3). This mission formalizes the third part, which the authors describe as the major contribution of the paper. Its point is that the error from the simulations is controlled nonasymptotically, step by step, without blowing up as the time step $h$ shrinks, even though every regression of the backward recursion reuses the same simulated paths.

Setting

A model consists of a horizon $T>0$, a drift $b$, a diffusion $\sigma$ satisfying the Lipschitz condition (H1), and a driver $f(t,x,y,z)$ satisfying (H2): $|f(t_2,x_2,y_2,z_2)-f(t_1,x_1,y_1,z_1)|\le C_f(|t_2-t_1|^{1/2}+|x_2-x_1|+|y_2-y_1|+|z_2-z_1|)$. For $N\ge1$, $h=T/N$, $t_k=kh$, the Euler scheme is $S^N_{t_{k+1}}=S^N_{t_k}+b(t_k,S^N_{t_k})h+\sigma(t_k,S^N_{t_k})\Delta W_k$, where the increments $\Delta W_k\sim\mathcal N(0,hI_q)$ are independent of the past. A chain $P^N_{t_k}\in\mathbb R^{d'}$ extends $S^N_{t_k}$, and $\Phi^N(P^N_{t_N})$ is the terminal value.

At each time $t_k$, function bases $p_{0,k}$ (for $Y$) and $p_{l,k}$, $1\le l\le q$ (for the components of $Z$) are fixed, orthonormal in the sense $\mathbb E[p_{l,k}(P^N_{t_k})p_{l,k}(P^N_{t_k})^*]=\mathrm{Id}$. The projection–Picard scheme (Definition 1) computes coefficients $\alpha^{i,I}_k$ by $I$ Picard iterations of an $\mathbf L_2$ least-squares problem, and sets $Y^{N,I,I}_{t_k}=\alpha^{I,I}_{0,k}\cdot p_{0,k}$, $Z^{N,I,I}_{l,t_k}=\alpha^{I,I}_{l,k}\cdot p_{l,k}$.

The empirical scheme (4) replaces the expectation by an average over $M$ independent simulations $(P^{N,m},\Delta W^m)$ of the path: $\alpha^{i,I,M}_k$ minimizes

$$\frac1M\sum_{m=1}^M\Big(Y^{N,I,I,M,m}_{t_{k+1}}-\alpha_0\cdot p^m_{0,k}+hf^m_k(\alpha^{i-1,I,M}_k)-\sum_{l=1}^q\alpha_l\cdot p^m_{l,k}\Delta W^m_{l,k}\Big)^2 .$$

The outputs are truncated: with $\rho^N_{l,k}(x)=\max(1,C_0|p_{l,k}(x)|)$ and a smooth profile $\xi$ equal to the identity on $[-3/2,3/2]$, $Y^{N,I,I,M}_{t_k}=\hat\rho^N_{0,k}(\alpha^{I,I,M}_{0,k}\cdot p_{0,k})$ with $\hat\rho^N_{l,k}(x)=\rho^N_{l,k}(P^N_{t_k})\xi(x/\rho^N_{l,k}(P^N_{t_k}))$. The regression vector is $[v_k]^*=(p_{0,k}^*,p_{1,k}^*\Delta W_{1,k}/\sqrt h,\dots,p_{q,k}^*\Delta W_{q,k}/\sqrt h)$, and the good event $\mathbf A^M_k$ (27) asks that the empirical matrices $V^M_j=\frac1M\sum_mv^m_j[v^m_j]^*$ and $P^M_{l,j}=\frac1M\sum_m p^m_{l,j}[p^m_{l,j}]^*$ be close to the identity for all $j\ge k$.

Formalization targets

Goal: Theorem 3

For $I\ge3$, orthonormal bases with $\mathbb E|p_{l,k}|^4<\infty$, $C_0$ such that the bounds of Proposition 2 hold, and $h$ small enough, for $0\le k\le N-1$,

$$\mathbb E|Y^{N,I,I}_{t_k}-Y^{N,I,I,M}_{t_k}|^2+h\sum_{j=k}^{N-1}\mathbb E|Z^{N,I,I}_{t_j}-Z^{N,I,I,M}_{t_j}|^2\le 9\sum_{j=k}^{N-1}\mathbb E(|\rho^N_j|^2\mathbf 1_{[\mathbf A^M_k]^c})+Ch^{I-1}\sum_{j=k}^{N-1}[1+|S_0|^2+\mathbb E|\rho^N_j|^2]+\frac{C}{hM}\sum_{j=k}^{N-1}\epsilon_j,$$

where $\epsilon_j$ collects second moments of $v_jv_j^*-\mathrm{Id}$, $|v_j|^2|p_{0,j+1}|^2$ and $|v_j|^2(1+|S^N_{t_j}|^2+\dots)$, written out in full in the goal statement. The constant $C$ and the threshold on $h$ depend only on the model and on $\xi$.

Milestones

1. Proposition 2: the a priori bounds $|Y^{N,i,I}_{t_k}|\le\rho^N_{0,k}$, $\sqrt h|Z^{N,i,I}_{l,t_k}|\le\rho^N_{l,k}$ that fix the truncation levels.
2. (28)–(29): the empirical least-squares solution and its contraction inequality $\lambda_{\min}(V^M)|\theta_x|^2\le|\theta_x\cdot v|^2_M\le|x|^2_M$.
3. Lemma 1: on $\mathbf A^M_k$ the empirical Picard iterations contract at rate $Ch$ to a unique fixed point, with error $[Ch]^I$ after $I$ steps.
4. (32): the pathwise bound $|\theta^{i,I,M}_k|^2\le C(\mathcal A^{N,M}_{k+1}+h\mathcal B^{N,M}_k)$ on $\mathbf A^M_k$.
5. (34): the expectation formula $\theta^{\infty,I}_k=\mathbb E(v_k[Y^{N,I,I}_{t_{k+1}}+hf_k(\alpha^{\infty,I}_k)])$.

Significance

Theorem 3 is nonasymptotic: together with Theorems 1 and 2 it lets one compare the three error sources and choose $h$, the bases and $M$ jointly for a target accuracy. The $1/(hM)$ rate shows how many paths a finer time grid requires, and the $h^{I-1}$ term shows that $I=3$ Picard iterations suffice. The term involving $[\mathbf A^M_k]^c$ isolates the event on which the empirical regression matrices are badly conditioned, which the truncation keeps under control.

The result is proved in the paper; nothing in this mission is open mathematically. No machine-checked proof of any part of it is known. A formal proof would check a long chain of estimates whose constants the paper tracks only as a generic $C$, and would settle the two misprints in the printed statement (see below). The definitions layer (Euler scheme, regression schemes, empirical regression matrices) is reusable for other regression Monte Carlo schemes for BSDEs and for optimal stopping.

Difficulty

The obvious approach is to treat each regression as an independent statistical estimation problem and apply a variance bound per time step. This fails because all regressions of the backward recursion use the same $M$ paths: the response at time $t_{k+1}$ is itself a function of the simulations, so the regression at $t_k$ is not a regression of a fixed variable on independent samples. Moreover, the empirical matrix $V^M_k$ may be singular, and on that event the empirical coefficients are unbounded unless truncated. A naive per-step bound also produces a factor $1/h$ at each of the $N=T/h$ steps, which explodes; the proof must keep the accumulated constants of order one.

Formalization scope

• All objects are defined in the namespace RegMCBSDE.Simulation. Continuous time is not used: the increments $\Delta W_k$ are any family that is $\mathcal F_{k+1}$-measurable, independent of $\mathcal F_k$ and $\mathcal N(0,hI_q)$-distributed for some filtration. That covers the Brownian case. (H3), a condition on the terminal functional of the continuous path, is dropped. So is the Markov representation of $P^N$, which no statement uses.
• Expectations of squared quantities are taken in $[0,\infty]$, on both sides of every inequality.
• Euclidean norms are written as sums of squares. $\|A\|\le c$ for symmetric $A$ is written as $|x^*Ax|\le c|x|^2$ for all $x$.
• The schemes are defined as predicates (any minimizer of each least-squares problem), never through a matrix inverse. The empirical coefficients are required to be measurable functions of the simulations, and any such minimizer is allowed, as the paper says the choice is arbitrary.
• The reference path at which the fitted coefficients are evaluated is assumed independent of the $M$ simulations.
• Picard iterations are imposed for every $i\ge1$.
• "For $h$ small enough" is formalized as $T/N<h_0$. Every constant $C$ and every $h_0$ is chosen before $N$, $I$, $M$, $k$, $C_0$, the bases, $S_0$ and the probability space. A formalization in which $C$ is chosen after the scheme data would make every inequality with a positive right-hand side trivially true, and is ruled out.
• "$C_0$ large enough" in Theorem 3 is the hypothesis that the bounds of Proposition 2 hold for $C_0$.
• Proposition 2 is stated for $h$ small enough. The printed statement omits this condition, which its proof needs through the uniform bound (19).
• In (32) and Lemma 1, $\rho^N_{0,N}$, which the paper leaves undefined, is read through the terminal response $\Phi^N(P^{N,m}_{t_N})$.
• Corrections to Theorem 3 (both from the paper's proof, p. 21):
  1. The printed $j=N-1$ summand contains the undefined $p_{0,N}$. It is replaced by $\mathbb E(|v_{N-1}|^2|\Phi^N(P^N_{t_N})|^2)$.
  2. The printed factor $\mathbb E|\rho^N_{0,j}(P^N_{t_j})|^2$ next to $\mathbb E(|v_j|^2|p_{0,j+1}|^2)$ becomes $\mathbb E|\rho^N_{0,j+1}(P^N_{t_{j+1}})|^2$.

Contributions welcome: proofs of the milestones, in particular (28)–(29) (finite-dimensional linear algebra) and (34) (independence and orthonormality), and a construction showing that measurable minimizers of (4) exist.

Selected references

• E. Gobet, J.-P. Lemor, X. Warin, A regression-based Monte Carlo method to solve backward stochastic differential equations, Ann. Appl. Probab. 15(3), 2172–2202, 2005. https://doi.org/10.1214/105051605000000412, preprint https://arxiv.org/abs/math/0508491v1
• N. El Karoui, S. Peng, M. C. Quenez, Backward stochastic differential equations in finance, Math. Finance 7(1), 1–71, 1997. https://doi.org/10.1111/1467-9965.00022
• F. A. Longstaff, E. S. Schwartz, Valuing American options by simulation: a simple least-squares approach, Rev. Financ. Stud. 14(1), 113–147, 2001. https://doi.org/10.1093/rfs/14.1.113
• J. Zhang, A numerical scheme for BSDEs, Ann. Appl. Probab. 14(1), 459–488, 2004. https://doi.org/10.1214/aoap/1075828058
• B. Bouchard, N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stochastic Process. Appl. 111(2), 175–206, 2004. https://doi.org/10.1016/j.spa.2004.01.001

Motivation

Sample covariance matrices are the basic object of multivariate statistics: principal component analysis, signal detection and factor models all start from the eigenvalues of $S = \frac1M\sum_{k=1}^M \vec y_k \vec y_k^{\,*}$ computed from $M$ samples of an $N$-dimensional vector. When $N$ is comparable to $M$, the eigenvalues of $S$ no longer approximate those of the population covariance $\Sigma$. The question is then whether a few large population eigenvalues ("spikes") are visible in the sample spectrum at all, and how the top sample eigenvalue fluctuates when they are.

Baik, Ben Arous and Péché (Ann. Probab. 33 (2005) 1643–1697) answered this for complex Gaussian samples. They found a sharp threshold $1+\gamma^{-1}$, where $\gamma^2 = M/N$, now called the BBP phase transition. This mission formalizes the supercritical half of the answer, Theorem 1.1(b).

Timeline.

• 2000–2001: Johansson (Comm. Math. Phys. 209, 2000) and Johnstone (Ann. Statist. 29) proved that for $\Sigma = I$ the largest eigenvalue, centred at $(1+\gamma^{-1})^2$ and scaled by $M^{2/3}$, has the Tracy–Widom law. Johansson treated the complex case, Johnstone the real case.
• 2005: Baik, Ben Arous and Péché treated $\Sigma$ with finitely many eigenvalues different from $1$. Spikes at or below $1+\gamma^{-1}$ give $M^{2/3}$ fluctuations with limit laws $F_k$ (part (a), a separate mission of this series). Spikes above $1+\gamma^{-1}$ give $\sqrt M$ fluctuations with limit $G_k$ (part (b)).
• 2006: Baik and Silverstein (J. Multivariate Anal. 97) located the outlying eigenvalue for general, non-Gaussian samples. Bai and Yao (Ann. IHP 44, 2008) obtained Gaussian fluctuations of the outliers for general samples.

Setting

Let $g_{kj}$, $1 \le k \le M$, $1 \le j \le N$, be independent standard complex Gaussians: $g = a + ib$ with $a, b$ independent real normal variables of mean $0$ and variance $1/2$. Fix a unitary matrix $U$ and positive population eigenvalues $\ell_1, \dots, \ell_N$. The samples are

$$\vec y_k = U\,\mathrm{diag}(\sqrt{\ell_1},\dots,\sqrt{\ell_N})\,\vec g_k,$$

which are independent mean-zero complex Gaussian vectors with covariance $\Sigma = U\,\mathrm{diag}(\ell)\,U^*$. The sample covariance matrix is $S = \frac1M\sum_k \vec y_k\vec y_k^{\,*}$, and $\lambda_1$ is its largest eigenvalue.

The regime has $M, N \to \infty$ with $M/N = \gamma^2$ and $\gamma$ in a compact subset of $[1,\infty)$. A fixed number $r$ of the $\ell_j$ differ from $1$. For some $1 \le k \le r$ the top $k$ coincide, $\ell_1 = \cdots = \ell_k$, with common value in a compact subset of $(1+\gamma^{-1},\infty)$. The others, $\ell_{k+1},\dots,\ell_r$, lie in a compact subset of $(0,\ell_1)$.

The limit law is the finite GUE distribution. Let $Z_k = \int_{\mathbb R^k}\prod_{i<j}|\xi_i-\xi_j|^2\prod_i e^{-\xi_i^2/2}\,d\xi$. Then

$$G_k(x) = \frac1{Z_k}\int_{(-\infty,x]^k}\prod_{1\le i<j\le k}|\xi_i-\xi_j|^2\prod_{i=1}^k e^{-\xi_i^2/2}\,d\xi_1\cdots d\xi_k,$$

which is the law of the largest eigenvalue of a $k\times k$ Gaussian unitary ensemble matrix (Definition 1.2).

Formalization targets

Goal: Theorem 1.1(b)

$$\mathbb P\left(\Big(\lambda_1 - \Big(\ell_1 + \frac{\ell_1\gamma^{-2}}{\ell_1-1}\Big)\Big)\cdot\frac{\sqrt M}{\sqrt{\ell_1^2 - \ell_1^2\gamma^{-2}/(\ell_1-1)^2}} \le x\right) \longrightarrow G_k(x)\qquad\text{for every } x \in \mathbb R.$$

Companion: Corollary 1.1(b)

$$\lambda_1 - \ell_1\Big(1 + \frac{\gamma^{-2}}{\ell_1-1}\Big) \longrightarrow 0 \quad\text{in probability.}$$

Milestones

The milestones follow the proof in §4, in attack order:

• the critical points of the phase $f$, (218)–(219);
• the contour estimates, Lemmas 4.1 and 4.2;
• the uniform kernel asymptotics, Proposition 4.1;
• the Hermite-polynomial identities (287), (295), (288) and (292);
• the closed form (299) of the limiting kernel $K_2$;
• the Fredholm-determinant representation of $G_k$, Lemma 1.1:

$$G_k(x) = \det\big(1 - \mathbf H^{(k)}_x\big),\qquad H^{(k)}(u,v) = \frac{c_{k-1}}{c_k}\,\frac{p_k(u)p_{k-1}(v) - p_{k-1}(u)p_k(v)}{u-v}\,e^{-(u^2+v^2)/4}.$$

Here $p_n$ are the orthonormal polynomials for the weight $e^{-x^2/2}$ and $c_n$ their leading coefficients.

Significance

The result. Below the threshold the top eigenvalue sticks to the bulk edge $(1+\gamma^{-1})^2$. Above it, Theorem 1.1(b) and Corollary 1.1(b) show that $\lambda_1$ separates from the bulk to the explicit location $\ell_1(1+\gamma^{-2}/(\ell_1-1))$. Its fluctuations shrink from order $M^{-2/3}$ to order $M^{-1/2}$, and their law is the top eigenvalue of a $k\times k$ GUE, with $k$ the multiplicity of the spike. This gives:

• a detection threshold for spiked signals in high dimension;
• the centring and scaling of tests based on the top sample eigenvalue;
• the first instance of a $k$-dependent family of finite-GUE limits in a spiked model.

Formalizing it. The result is proved in the paper and is not open. Nothing in this mission has a machine-checked proof yet. A complete development would include:

• the first formal statements of the spiked complex Wishart model;
• the GUE law $G_k$ and its Fredholm representation;
• a uniform steepest-descent analysis with explicit contours.

The two contour lemmas (4.1, 4.2) and the Hermite identities are elementary and are footholds. Proposition 4.1 and Lemma 1.1 are substantial. The paper does not prove Lemma 1.1 but cites it as a standard result of random matrix theory (its references [28, 41]).

Difficulty

The obvious route is to diagonalise the sample matrix, write down the eigenvalue density and take a limit. That density involves the Harish-Chandra–Itzykson–Zuber integral. Its large-$N$ limit is not accessible directly, because the spike enters through a determinant with $N$ nearly coincident columns.

The paper instead starts from an exact Fredholm-determinant formula (Proposition 2.1, posed in mission 1 of this series). In it the kernel factors into two contour integrals $\mathcal H$ and $\mathcal J$ with phase $f(z) = -\mu(z-q) + \log z - \gamma^{-2}\log(1-z)$. Above the threshold the two factors are governed by different points. $\mathcal J$ has a nondegenerate saddle at $\pi_1 = \ell_1^{-1}$. For $\mathcal H$ the natural saddle $1/(\mu\pi_1)$ lies beyond the pole at $\pi_1$, so the contour must be deformed through a pole of order $k$, and the leading term is a residue rather than a saddle contribution. The estimates must hold uniformly in $\gamma$ and in the remaining spikes, and the convergence must be strong enough (Hilbert–Schmidt) to pass to Fredholm determinants.

Formalization scope

Model. Complex Gaussian samples with mean zero and no centring. $S = \frac1M\sum_k\vec y_k\vec y_k^{\,*}$, with the factor $1/M$. $\Sigma = U\,\mathrm{diag}(\ell)\,U^*$ for every unitary $U$. This is the model of the paper's (59), (61) and Proposition 2.1; the introduction's mentions of $1/N$, of centring and of the real density (1) are inconsistent with those formulas. $\lambda_1$ is the supremum of the eigenvalues of the Hermitian matrix $S$, and probabilities are measures of sets of sample arrays.

Regime. The regime is in sequence form:

• $N_n \to \infty$ and $\gamma_n = \sqrt{M_n/N_n} \in [1,\gamma_0]$;
• $1+\gamma_n^{-1}+c \le \ell_1 = \cdots = \ell_k \le C$;
• $c \le \ell_j \le \ell_1 - c$ for $k < j \le r$;
• $\ell_j = 1$ for $j > r$.

The fixed margins $c, C$ encode the compact subsets whose open ends move with $\gamma$. Convergence is pointwise in $x$.

Analytic conventions.

• The Fredholm determinant is the Fredholm series $\sum_n\frac{(-1)^n}{n!}\int_{(x,\infty)^n}\det[K(u_i,u_j)]$.
• $H^{(k)}$ takes its continuous diagonal value.
• $p_n = \mathrm{He}_n/((2\pi)^{1/4}\sqrt{n!})$ with Mathlib's probabilists' Hermite polynomials, equal to the paper's (31).
• The closed contours $\Gamma$, $\Sigma$ in $\mathcal H$, $\mathcal J$ are explicit circles satisfying the paper's constraints.
• $\Sigma_\infty$ is the imaginary axis.
• Residues of $a^{-k}\varphi(a)$ are Taylor coefficients.

Ruled-out trivialisations. Five shortcuts would make the statements trivial, and none is available:

• division by zero on the kernel diagonal is replaced by the diagonal value;
• the Fredholm series is shown to be $1$ on the zero kernel, so it is not identically junk;
• the conditionally convergent real form of contour integrals is not used;
• $\Sigma$ is not specialised to a diagonal matrix;
• the regime is shown satisfiable in a sorry-free check.

Infrastructure. The needed infrastructure is the complex Wishart model, eigenvalues of random Hermitian matrices, Fredholm determinants of integral operators, and contour integrals with uniform estimates. The Hermite and Fredholm layers can be reused for any orthogonal-polynomial ensemble. Contributions welcome:

• proofs of the elementary milestones ((218)–(219), (287), (295), (288), (292), (299), Lemmas 4.1, 4.2);
• Lemma 1.1 via Christoffel–Darboux and Andréief;
• the operator-theoretic step from Hilbert–Schmidt convergence of kernels to convergence of Fredholm series.

Selected references

• J. Baik, G. Ben Arous, S. Péché, Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices, Ann. Probab. 33(5) (2005) 1643–1697. https://doi.org/10.1214/009117905000000233
• K. Johansson, Shape fluctuations and random matrices, Comm. Math. Phys. 209 (2000) 437–476. https://doi.org/10.1007/s002200050027
• I. M. Johnstone, On the distribution of the largest eigenvalue in principal components analysis, Ann. Statist. 29 (2001) 295–327. https://doi.org/10.1214/aos/1009210544
• J. Baik, J. W. Silverstein, Eigenvalues of large sample covariance matrices of spiked population models, J. Multivariate Anal. 97 (2006) 1382–1408. https://doi.org/10.1016/j.jmva.2005.08.003
• Z. Bai, J. Yao, Central limit theorems for eigenvalues in a spiked population model, Ann. Inst. H. Poincaré Probab. Statist. 44 (2008) 447–474. https://doi.org/10.1214/07-AIHP118

Motivation

A standard way to build a random network with a prescribed degree sequence is to choose a graph uniformly at random from the set of all simple graphs with that degree sequence. The model arises in testing whether the exponential family with the degree sequence as sufficient statistic fits network data, in simulating networks with given degrees, and, for constant degrees, as the random regular graph (Blitzstein–Diaconis, Chatterjee–Diaconis–Sly). Quantities such as the expected number of triangles of such a graph are usually estimated by simulation.

For dense graphs, those whose number of edges is of order $n^2$, the Lovász–Szegedy theory of graph limits (Lovász–Szegedy 2006) describes a large graph by a symmetric function $W:[0,1]^2\to[0,1]$. Chatterjee, Diaconis and Sly identify this limit for uniformly random graphs with a given degree sequence. The answer is the limit of the β-model, a random graph with independent edges whose edge probabilities are $e^{\beta_i+\beta_j}/(1+e^{\beta_i+\beta_j})$. The limit gives exact asymptotic formulas for subgraph counts without simulation.

Setting

Graphs and densities. For a finite simple graph $H$ on $[k]=\{1,\dots,k\}$ and a simple graph $G$ on $n$ vertices, the homomorphism density is

$$t(H,G)=\frac{|\hom(H,G)|}{n^k},$$

where $\hom(H,G)$ is the set of edge-preserving maps $V(H)\to V(G)$. For $W:[0,1]^2\to\mathbb R$,

$$t(H,W)=\int_{[0,1]^k}\prod_{\{i,j\}\in E(H)}W(x_i,x_j)\,dx_1\cdots dx_k,$$

with one factor per edge. Graphs $G_n$ on $n$ vertices converge to the limit represented by $W$ if $t(H,G_n)\to t(H,W)$ for every finite simple graph $H$.

Degree sequences and scaling limits. For each $n$ let $d^n=(d^n_1\ge\dots\ge d^n_n)$ be the degree sequence of some simple graph on $n$ vertices. The sequence $\{d^n\}$ has scaling limit $f$, a nonincreasing function on $[0,1]$, if

$$\left|\frac{d^n_1}{n}-f(0)\right|+\left|\frac{d^n_n}{n}-f(1)\right|+\frac1n\sum_{i=1}^n\left|\frac{d^n_i}{n}-f\!\left(\frac in\right)\right|\longrightarrow0 .$$

$D'[0,1]$ is the set of nonincreasing functions on $[0,1]$ that are left continuous on $(0,1)$. It carries the modified $L^1$ norm $\|f\|_{1'}=|f(0)|+|f(1)|+\int_0^1|f|$. $\mathcal F\subseteq D'[0,1]$ is the set of all scaling limits of degree sequences.

The random graphs. $G_n$ is chosen uniformly from the simple graphs on $\{1,\dots,n\}$ with degree sequence $d^n$. In a random graph with independent edges each pair $\{i,j\}$ is an edge with probability $p_{ij}$, independently of the other pairs. The β-model $P_\beta$ is the case $p_{ij}=e^{\beta_i+\beta_j}/(1+e^{\beta_i+\beta_j})$.

Formalization targets

Goal: Theorem 1.1 (pp. 4–5)

If $f$ lies in the interior of $\mathcal F$ for $\|\cdot\|_{1'}$, there is a function $g\in D'[0,1]$, unique on $[0,1]$, such that

$$W(x,y)=\frac{e^{g(x)+g(y)}}{1+e^{g(x)+g(y)}}\quad\text{satisfies}\quad f(x)=\int_0^1W(x,y)\,dy\ \ \text{for all }x\in[0,1],$$

and almost surely $t(H,G_n)\to t(H,W)$ for every finite simple graph $H$.

Milestones (in the order the proof uses them)

• Lemma 2.2 (p. 14): the fixed-point map $\varphi$ of the β-model likelihood equations satisfies $|\varphi(x)-\varphi(y)|_1\le2e^{2K}|x-y|_1$ when $|x|_\infty,|y|_\infty\le K$.
• §6.2, p. 30: for nonincreasing $d$, the Erdős–Gallai slack $\mathcal E(B)$ over sets of size $k$ is minimized at $\{1,\dots,k\}$.
• Lemma 6.1 (p. 26): for independent-edge random graphs, $\mathbb P(|t(H,G)-\mathbb E t(H,G)|>\varepsilon)\le2e^{-C\varepsilon^2n^2}$ with $C=C(H)$.
• Claim 6.3 (p. 27): integer margins close to the row and column sums of a matrix with entries in $[\delta,1-\delta]$ are realized by a 0–1 contingency table.
• Lemma 6.2 (p. 26): if $\delta\le p_{ij}\le1-\delta$ and $d_i=\sum_{j\ne i}p_{ij}$, the independent-edge graph has degree sequence $d$ with probability at least $\tfrac12\delta^{\,n^{3/2+\varepsilon}}$ for large $n$.
• p. 34: conditioned on its degree sequence being $d$, the β-model graph is uniform on graphs with degree sequence $d$.

The proof also uses three results posed in sibling missions of this series: Proposition 1.2 (mission 4, the interior of $\mathcal F$), Lemma 4.1 (mission 3, existence and boundedness of the MLE) and Theorem 1.5 (mission 1, geometric convergence of the iteration $x\mapsto\varphi(x)$). They are not restated here.

Significance

The theorem turns questions about a uniformly random graph with given degrees, a law with no independence, into questions about an explicit graphon. Every subgraph density, the limiting degree distribution, and every graph parameter continuous for the graph-limit topology of the uniform model is computed from $W$. The function $g$ solves a continuum version of the β-model maximum-likelihood equations, which connects the combinatorial model with exponential-family statistics.

The result is proved in the paper, which is published in the Annals of Applied Probability (2011). To our knowledge it has no machine-checked proof. Formalizing it requires a Lean development of dense graph limits (homomorphism densities, graphon densities, graph-limit convergence), none of which exists in Mathlib. It also requires a bounded-differences concentration inequality, and a transfer argument from independent-edge models to models conditioned on degrees. Each of these is reusable beyond this mission.

Difficulty

The uniform measure on graphs with a fixed degree sequence has no independence, so concentration inequalities do not apply to it directly. The obvious route is to condition an independent-edge model on its degree sequence. This works only if the probability of the conditioning event is larger than the concentration bound, which is $e^{-cn^2}$. A crude lower bound of $\delta^{\binom n2}$ is too small. Lemma 6.2 supplies $\delta^{n^{3/2+\varepsilon}}$, and it needs a combinatorial completion step (Claim 6.3).

The second difficulty is analytic. The maximum-likelihood parameters $\beta^n$ for different $n$ live in different dimensions. One must show that their step-function profiles converge in $L^1$ to a limit $g$ and that this $g$ solves the continuum equation at every point of $[0,1]$, not only almost everywhere. Several steps of the paper's argument are only sketched. Among them are the almost-sure convergence of the β-model graphs to $W$ and the estimate $f_n(x)=d^n_{\lceil nx\rceil}/n+O(1/n)$ (pp. 33–34).

Formalization scope

• Vertices are Fin n (the paper's vertex $i$ is $i-1$), graphs are SimpleGraph (Fin n), and a test graph $H$ is a SimpleGraph (Fin k), which loses nothing because densities are invariant under relabelling.
• Functions on $[0,1]$ are ℝ → ℝ, and only their values on $[0,1]$ matter. Uniqueness of $g$ is equality on $[0,1]$.
• $t(H,W)$ takes one factor per unordered edge, and the integral is the Lebesgue integral over $[0,1]^k$.
• Laws on graphs are finite sums of Dirac masses. The σ-algebra on SimpleGraph (Fin n) is discrete.
• The interior of $\mathcal F$ is taken in $D'[0,1]$ for $\|\cdot\|_{1'}$, as the ball form "every $h\in D'[0,1]$ with $\|h-f\|_{1'}<\varepsilon$ lies in $\mathcal F$".
• The paper does not say how the $G_n$ for different $n$ are coupled. The goal quantifies over every probability space and every family of measurable $G_n$ with the uniform marginals, so the almost-sure convergence holds for every joint law. The almost-sure event contains "for every $H$".
• Lemma 6.2's printed bound $\tfrac12\exp(-\log(\delta)n^{3/2+\varepsilon})$ exceeds $1$. The milestone states the intended bound $\tfrac12\exp(\log(\delta)n^{3/2+\varepsilon})$, with $N$ depending only on $(\delta,\varepsilon)$.
• Constants are explicit quantifiers. In Lemma 6.1, $C>0$ is chosen from $H$ before $n$, $p$ and $\varepsilon$.

A formalization that replaces almost-sure convergence by convergence in expectation or in probability, that fixes a single coupling, or that asserts the existence of $g$ only almost everywhere states a weaker theorem and does not close this mission.

Welcome contributions include a general graph-limit library (densities, convergence, the counting lemma), McDiarmid's bounded-differences inequality for finitely many independent Bernoulli variables, the Gale–Ryser theorem in the form needed by Claim 6.3, and proofs of the milestones.

Selected references

• S. Chatterjee, P. Diaconis, A. Sly, Random graphs with a given degree sequence, Ann. Appl. Probab. 21(4), 1400–1435, 2011. https://arxiv.org/abs/1005.1136v5
• L. Lovász, B. Szegedy, Limits of dense graph sequences, J. Combin. Theory Ser. B 96, 933–957, 2006. https://doi.org/10.1016/j.jctb.2006.05.002
• C. Borgs, J. Chayes, L. Lovász, V. Sós, K. Vesztergombi, Convergent sequences of dense graphs I, Adv. Math. 219, 1801–1851, 2008. https://doi.org/10.1016/j.aim.2008.07.008
• J. Blitzstein, P. Diaconis, A sequential importance sampling algorithm for generating random graphs with prescribed degrees, Internet Math. 6(4), 489–522, 2011. https://doi.org/10.1080/15427951.2010.557277