Statistics

Motivation

Many procedures in high-dimensional statistics reduce to computing the probability that a normalized sum of independent random vectors falls into a set: simultaneous confidence intervals for many means, max-type tests, multiple testing with family-wise error control, and inference after model selection. When the dimension $p$ is comparable to or much larger than the sample size $n$, these probabilities are approximated by those of a Gaussian vector with the same covariance, the high-dimensional central limit theorem. That Gaussian law is not available in practice, because its covariance matrix is unknown. The Gaussian multiplier bootstrap replaces it by a Gaussian vector built from the data, whose covariance is the sample covariance, and simulates from it.

Chernozhukov, Chetverikov and Kato, Central limit theorems and bootstrap in high dimensions (Ann. Probab. 45 (2017)), prove that this bootstrap is valid uniformly over hyperrectangles and over simple convex sets (sets sandwiched between a polytope with polynomially many facets and its small enlargement), with error depending on $p$ only through $\log p$. This mission formalizes their abstract multiplier bootstrap theorem, Theorem 4.1.

Timeline:

• 2013: the same authors prove Gaussian approximation and multiplier bootstrap validity for maxima of sums of high-dimensional vectors (Ann. Statist. 41 (2013)).
• 2015: comparison and anti-concentration bounds for maxima of Gaussian vectors (Probab. Theory Related Fields 162 (2015)), reference [20] of the paper; its Theorem 1 is the Gaussian-to-Gaussian comparison used here.
• 2017: the present paper extends both the CLT and the bootstrap from maxima to hyperrectangles and simple convex sets.

Setting

Let $n\ge4$ and $p\ge3$. Let $X_1,\dots,X_n$ be independent random vectors in $\mathbb R^p$ with coordinates $X_{ij}$, centred ($\mathrm E[X_{ij}]=0$) with $\mathrm E[X_{ij}^2]<\infty$. Let $Y_1,\dots,Y_n$ be independent Gaussian vectors with $Y_i\sim N(0,\mathrm E[X_iX_i'])$ and set

$$S_n^Y=\frac1{\sqrt n}\sum_{i=1}^nY_i .$$

The multiplier bootstrap: let $e_1,\dots,e_n$ be i.i.d. $N(0,1)$, independent of the data $X_1^n=\{X_1,\dots,X_n\}$, let $\bar X=n^{-1}\sum_iX_i$, and

$$S_n^{eX}=\frac1{\sqrt n}\sum_{i=1}^ne_i(X_i-\bar X).$$

Given the data, $S_n^{eX}$ is Gaussian with covariance $\widehat\Sigma=n^{-1}\sum_i(X_i-\bar X)(X_i-\bar X)'$, while $S_n^Y$ has covariance $\Sigma=n^{-1}\sum_i\mathrm E[X_iX_i']$.

A hyperrectangle is a set $\{w:a_j\le w_j\le b_j\ \forall j\}$ with $-\infty\le a_j\le b_j\le\infty$. For a finite set $\mathcal V$ of unit vectors and thresholds $s(v)$, the polyhedron $A^m=\bigcap_{v\in\mathcal V}\{w:w'v\le s(v)\}$ has the enlargement $A^{m,\epsilon}=\bigcap_{v\in\mathcal V}\{w:w'v\le s(v)+\epsilon\}$. For constants $a,d>0$, a Borel set $A$ satisfies condition (C) if $A^m\subseteq A\subseteq A^{m,a/n}$ for some such polyhedron with $m\le(pn)^d$ facets; $A^m(A)$ denotes this polyhedron and $\mathcal V(A^m)$ its normals. Condition (M.1′) with constant $b>0$ asks $n^{-1}\sum_i\mathrm E[(v'X_i)^2]\ge b$ for every $v\in\mathcal V(A^m)$. Finally,

$$\Delta_n(\mathcal A)=\sup_{A\in\mathcal A}\max_{v_1,v_2\in\mathcal V(A^m(A))}\big|v_1'(\widehat\Sigma-\Sigma)v_2\big|,\qquad \Delta_{n,r}=\max_{j,k}|\widehat\Sigma_{jk}-\Sigma_{jk}| .$$

Formalization targets

Goal: Theorem 4.1

Let $\mathcal A$ be a class of sets satisfying (C) and (M.1′). There is $C$ depending only on $a,b,d$ such that for every $\bar\Delta_n>0$, on the event $\Delta_n(\mathcal A)\le\bar\Delta_n$,

$$\rho_n^{MB}(\mathcal A)=\sup_{A\in\mathcal A}\big|P(S_n^{eX}\in A\mid X_1^n)-P(S_n^Y\in A)\big|\le C\big\{\bar\Delta_n^{1/3}\log^{2/3}(pn)+n^{-1}\log^{1/2}(pn)\big\}.$$

The theorem is deterministic in the data: it holds at every realization on the event, and turning it into a rate requires only a bound on $\Delta_n(\mathcal A)$.

Milestones

In the order of the proof (App. E.2, p. 2342):

1. Conditional Gaussianity: given $X_1^n$, $S_n^{eX}\sim N(0,\widehat\Sigma)$.
2. Display (21): $0\le F_\beta(w)-\max_j(w_j-y_j)\le\beta^{-1}\log p$ for the smooth max $F_\beta(w)=\beta^{-1}\log\sum_je^{\beta(w_j-y_j)}$.
3. The comparison display $|\mathrm E[g(F_\beta(S^{eX}_n))\mid X_1^n]-\mathrm E[g(F_\beta(S^Y_n))]|\le(\|g''\|_\infty/2+\beta\|g'\|_\infty)\Delta_{n,r}$.
4. Lemma A.1 (Nazarov's inequality): $P(Y\le y+a)-P(Y\le y)\le Ca\sqrt{\log p}$ for a centred Gaussian $Y$ with variances at least $b$.
5. The smoothed step: $|P(S^{eX}_n\le y-\phi^{-1}\mid X_1^n)-P(S^Y_n\le y-\phi^{-1})|\le C\{\phi^{-1}\log^{1/2}p+(\phi^2+\beta\phi)\Delta_{n,r}\}$ with $\beta=\phi\log p$.
6. Display (39): $\sup_y|P(S^{eX}_n\le y\mid X_1^n)-P(S^Y_n\le y)|\le C\Delta_{n,r}^{1/3}\log^{2/3}p$ under (M.1): $n^{-1}\sum_i\mathrm E[X_{ij}^2]\ge b$.
7. Remark 4.1 = (38): $\rho_n^{MB}(\mathcal A^{re})\le C\bar\Delta_n^{1/3}\log^{2/3}p$ on $\Delta_{n,r}\le\bar\Delta_n$, $C$ depending only on $b$.
8. The reduction display: for one simple convex set, the error at $A$ is at most $C\epsilon\log^{1/2}(pn)+\bar\rho$, where $\bar\rho$ is the larger error at $A^m$ and $A^{m,\epsilon}$.

Significance

Theorem 4.1 is the abstract step behind the paper's explicit multiplier bootstrap rates: Proposition 4.1 combines it with maximal inequalities that bound $\Delta_n(\mathcal A)$ with high probability, and obtains, with probability at least $1-\alpha$, the rate $(B_n^2\log^5(pn)\log^2(1/\alpha)/n)^{1/6}$ for simple convex sets with sparse facet normals. In applications it justifies bootstrap critical values for max-type statistics and simultaneous confidence regions when $p\gg n$. Remark 4.1, the hyperrectangle case, is the version most used in practice.

The results are proved on paper; none is formalized on the platform. A formal development produces machine-checked versions of the smooth-max sandwich, a Gaussian anti-concentration bound (Nazarov's inequality), a quantitative Gaussian comparison inequality and the bootstrap theorem itself. The comparison display and Lemma A.1 are proved in the paper only by citation ([20] and Klivans–O'Donnell–Servedio, Thm. 20), so formal proofs of them are new work rather than transcriptions.

Difficulty

The obvious route compares $S^{eX}_n$ and $S^Y_n$ through the Kolmogorov distance of their laws or through a total-variation bound between two Gaussians. Both depend polynomially on $p$ and fail when $p\gg n$. The proof instead compares smooth functionals of the coordinatewise maximum, where the error depends on $p$ only through $\beta\sim\log p$, and then converts back to probabilities of orthants with an anti-concentration inequality whose constant grows like $\sqrt{\log p}$. Both conversions must hold uniformly in the threshold $y$, and the final bound for simple convex sets requires passing to the $m$-dimensional vectors $(v'X_i)_{v\in\mathcal V}$ with $m\le(pn)^d$, which is why $\log(pn)$ appears. The comparison inequality for $g\circ F_\beta$ is the step with the least existing infrastructure: it is a Slepian–Stein interpolation between two Gaussian laws with explicit control of second derivatives of $F_\beta$.

Formalization scope

• Vectors live in EuclideanSpace ℝ (Fin p); $N(0,\Sigma)$ is Mathlib's multivariateGaussian 0 Σ. Independence is iIndepFun; $Y_i$'s law is fixed by P.map (Y i) = multivariateGaussian 0 (covMat P (X i)), and square integrability of each $X_{ij}$ makes covMat positive semidefinite, so this is the genuine Gaussian law.
• $n\ge4$ and $p\ge3$ (the paper's standing assumption) are hypotheses of every theorem.
• $P(\cdot\mid X_1^n)$ is read per realization: the bootstrap law mbLaw x is the push-forward of $N(0,1)^{\otimes n}$ under $e\mapsto n^{-1/2}\sum_ie_i(x_i-\bar x)$, and the bound is asserted at every $\omega$ on the event, for the data $X(\omega)$. Since $e$ is independent of the data, this is a version of the conditional law; no conditional expectation appears.
• Suprema over the class, over $v_1,v_2$ and over hyperrectangles are encoded as bounds for every member.
• The class $\mathcal A$ is an indexed family of measurable sets, each with its approximating polyhedron given by a finite set of unit vectors and real thresholds. The page uses facet normals and support values; allowing any finite family of unit vectors is a mild generalization that implies the page's theorem.
• Constants: $C$ is quantified after $a,b,d$ (after $b$ alone for (38), (39), Lemma A.1 and the smoothed step) and before $n$, $p$, the probability space, the data, the class, $\bar\Delta_n$ and the realization. A constant chosen after $n$ or $p$ would trivialize the goal.
• Moments enter only through lower bounds ((M.1), (M.1′)) on integrals of square-integrable functions, so no integral takes a default value; the comparison display integrates continuous functions of linear growth against Gaussian laws. The goal does not mention $F_\beta$, $g$, $\varrho_n^{MB}$, (38) or (39).

Infrastructure needed: multivariate Gaussian laws and their linear images, Gaussian integration by parts or Slepian interpolation, Nazarov's inequality, and properties of log-sum-exp. Lemma A.1, display (21) and the comparison display are reusable beyond this mission, for the high-dimensional CLT of the companion mission and for any max-type Gaussian approximation. Proofs of any milestone are welcome, as are proofs of the conditional Gaussianity identity from Mathlib's Gaussian API.

Selected references

• V. Chernozhukov, D. Chetverikov, K. Kato, Central limit theorems and bootstrap in high dimensions, Ann. Probab. 45(4):2309–2352, 2017. https://doi.org/10.1214/16-AOP1113
• V. Chernozhukov, D. Chetverikov, K. Kato, Comparison and anti-concentration bounds for maxima of Gaussian random vectors, Probab. Theory Related Fields 162:47–70, 2015. https://doi.org/10.1007/s00440-014-0565-9
• V. Chernozhukov, D. Chetverikov, K. Kato, Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors, Ann. Statist. 41(6):2786–2819, 2013. https://doi.org/10.1214/13-AOS1161
• A. Klivans, R. O'Donnell, R. Servedio, Learning geometric concepts via Gaussian surface area, FOCS 2008. https://doi.org/10.1109/FOCS.2008.64
• F. Nazarov, On the maximal perimeter of a convex set in ℝⁿ with respect to a Gaussian measure, Geometric Aspects of Functional Analysis, Lecture Notes in Math. 1807, 169–187, 2003. https://doi.org/10.1007/978-3-540-36428-3_14

Motivation

Many statistical procedures look at a large number of averages at once: simultaneous confidence intervals for $p$ means, multiple testing of $p$ hypotheses, the maximum of $p$ $t$-statistics. Their calibration needs the joint law of a $p$-dimensional normalized sum, and in modern applications $p$ is comparable to, or much larger than, the sample size $n$. Classical multivariate central limit theorems with explicit error, such as Bentkus's bound over convex sets (Bentkus 2003), require $p$ to grow slowly with $n$ (roughly $\sqrt p = o(n^{1/7})$), which rules out these applications.

Chernozhukov, Chetverikov and Kato showed that when the class of sets is restricted to hyperrectangles, the Gaussian approximation error depends on the dimension only through $\log p$ (Ann. Probab. 45 (2017)).

Timeline. In 2013 the same authors proved a Gaussian approximation for the maximum of a sum (sets of the form $\{w : w_j \le a \ \forall j\}$) with rate $n^{-1/8}$ (Ann. Statist. 41 (2013)). The 2017 paper extends this to all hyperrectangles, improves the rate to $n^{-1/6}$, and adds bootstrap results.

Setting

Let $n \ge 4$ and $p \ge 3$. Let $X_1, \dots, X_n$ be independent random vectors in $\mathbb R^p$ with coordinates $X_{ij}$, each centred ($\mathrm E[X_{ij}] = 0$) with $\mathrm E[X_{ij}^2] < \infty$. Let $Y_1, \dots, Y_n$ be independent centred Gaussian vectors with $Y_i \sim N(0, \mathrm E[X_iX_i'])$, the same covariance as $X_i$. The normalized sums are

$$S^X_n := \frac{1}{\sqrt n}\sum_{i=1}^n X_i, \qquad S^Y_n := \frac{1}{\sqrt n}\sum_{i=1}^n Y_i.$$

A hyperrectangle is a set $A = \{w \in \mathbb R^p : a_j \le w_j \le b_j \text{ for all } j\}$ with $-\infty \le a_j \le b_j \le \infty$; $\mathcal A^{\mathrm{re}}$ is the class of all of them. The quantity to bound is

$$\rho_n(\mathcal A^{\mathrm{re}}) := \sup_{A \in \mathcal A^{\mathrm{re}}}\big|\mathrm P(S^X_n \in A) - \mathrm P(S^Y_n \in A)\big|.$$

The error is measured through the third-moment parameter $L_n := \max_j n^{-1}\sum_i \mathrm E|X_{ij}|^3$ and the truncated maximal moments

$$M_{n,X}(\phi) := \frac1n\sum_{i=1}^n \mathrm E\Big[\max_j|X_{ij}|^3\,1\Big\{\max_j|X_{ij}| > \frac{\sqrt n}{4\phi\log p}\Big\}\Big],$$

$M_{n,Y}(\phi)$ the same with $Y$, and $M_n(\phi) := M_{n,X}(\phi) + M_{n,Y}(\phi)$. Here $\log$ is the natural logarithm.

Formalization targets

Goal: Theorem 2.1

If $n^{-1}\sum_i \mathrm E[X_{ij}^2] \ge b > 0$ for all $j$, there are constants $K_1, K_2 > 0$ depending only on $b$ such that for every $\bar L_n \ge L_n$,

$$\rho_n(\mathcal A^{\mathrm{re}}) \le K_1\left[\left(\frac{\bar L_n^2\log^7 p}{n}\right)^{1/6} + \frac{M_n(\phi_n)}{\bar L_n}\right], \qquad \phi_n := K_2\left(\frac{\bar L_n^2\log^4 p}{n}\right)^{-1/6}.$$

The constants are left existential, as on the page: the goal asserts the shape of the bound, not particular values.

Milestones

1. Display (21): the smooth maximum $F_\beta(w) = \beta^{-1}\log\sum_j e^{\beta(w_j - y_j)}$ is within $\beta^{-1}\log p$ of $\max_j(w_j - y_j)$.
2. Lemma B.1: association inequalities for nondecreasing nonnegative functions, including (20) without independence.
3. Lemma A.1 (Nazarov's inequality): $\mathrm P(Y \le y + a) - \mathrm P(Y \le y) \le Ca\sqrt{\log p}$ for a centred Gaussian $Y$ with $\mathrm E[Y_j^2] \ge b$.
4. Lemma 5.1, the key lemma: a recursive bound on $\varrho_n := \sup_{y, v\in[0,1]}|\mathrm P(\sqrt v S^X_n + \sqrt{1-v}S^Y_n \le y) - \mathrm P(S^Y_n \le y)|$.
5. Corollary 5.1: the same bound for $\varrho_n'$, the supremum over hyperrectangles.
6. Lemma C.1: an exponential tail bound gives $\mathrm E[\xi^3 1\{\xi > t\}] \le 6A(t+B)^3e^{-t/B}$.

Companion: Proposition 2.1

Under the moment conditions (M.1)–(M.2) and either an exponential (E.1) or a polynomial (E.2) moment bound with parameter $B_n \ge 1$,

$$\rho_n(\mathcal A^{\mathrm{re}}) \le C\Big(\frac{B_n^2\log^7(pn)}{n}\Big)^{1/6} \quad\text{or}\quad \rho_n(\mathcal A^{\mathrm{re}}) \le C\Big\{\Big(\frac{B_n^2\log^7(pn)}{n}\Big)^{1/6} + \Big(\frac{B_n^2\log^3(pn)}{n^{1-2/q}}\Big)^{1/3}\Big\}.$$

Significance

The result. Under the exponential moment condition (E.1), Proposition 2.1 makes the Gaussian approximation error tend to zero when $B_n^2\log^7(pn) = o(n)$, allowing $p$ much larger than $n$. Under the polynomial moment condition (E.2), its second error term must also tend to zero. Theorem 2.1 is the probabilistic basis of the paper's multiplier and empirical bootstrap results for simultaneous inference (Chernozhukov, Chetverikov and Kato 2017).

Formalizing it. The result is proved on paper; no machine-checked version of this theorem was found in the platform catalog during this mission's prior-art search. A formal proof needs interpolation between $S^X_n$ and $S^Y_n$, Gaussian anti-concentration (Nazarov's inequality), and third-order Taylor estimates of smooth approximations of indicator functions. Nazarov's inequality is used in the paper by citation, so its formal proof is a separate milestone.

Difficulty

The obvious route, a Lindeberg replacement of each $X_i$ by $Y_i$ after smoothing the indicator of a hyperrectangle coordinate by coordinate, loses a polynomial factor in $p$: the derivatives of a product of $p$ smoothed one-dimensional indicators add up over coordinates, so the error grows like a power of $p$ rather than of $\log p$. Any smoothing must also be paid for by the probability that the Gaussian vector falls in a thin neighbourhood of the boundary of a hyperrectangle, and a dimension-free bound on that probability is not available for general hyperrectangles; the $\sqrt{\log p}$ anti-concentration bound is the best one can use. Finally, a direct comparison yields an inequality in which the quantity to be bounded appears on both sides, and the rate $n^{-1/6}$ is obtained only after that inequality is solved; a one-shot estimate gives a worse power of $n$. The truncation in $M_n(\phi)$ handles coordinates too large for a third-order Taylor expansion, since no moment condition beyond third moments is assumed.

Formalization scope

Vectors live in EuclideanSpace ℝ (Fin p) and the coordinate $X_{ij}$ is X i ω j. The Gaussian comparison vectors satisfy P.map (Y i) = multivariateGaussian 0 (covMat P (X i)) with covMat P Z = (E[Z_j Z_k])_{j,k}, which is positive semidefinite, so this is the genuine $N(0, \mathrm E[X_iX_i'])$; the $Y_i$ are independent. Hyperrectangle endpoints are extended reals used only in comparisons. Every theorem carries the standing assumptions $n \ge 4$, $p \ge 3$, independence, centring and finite second moments. Items that involve $L_n$ or $M_n$ add $\mathrm E|X_{ij}|^3 < \infty$; when a third moment is infinite the paper's bounds are vacuous, so this loses nothing. $M_n(\phi)$ is defined by its formula for every real $\phi$, because Theorem 2.1 evaluates it at $\phi_n$, which may be below $1$. Lemma 5.1 and Corollary 5.1 assume the $Y$'s independent of the $X$'s, as the definition of $\varrho_n$ does; Theorem 2.1 and Proposition 2.1 do not. A bound on a supremum is stated set by set.

Constants that the paper says depend only on $b$ (or on $b$ and $q$) are quantified after $b$ and before $n$, $p$, the probability space and every random vector; a constant chosen after $n$, $p$ or the law would trivialize every statement. Upper bounds on expectations in hypotheses ((M.2), (E.1), (E.2)) and in Lemma C.1's conclusion are lower Lebesgue integrals, so a non-integrable function cannot make a hypothesis hold with a junk value of $0$. Dropping the independence of the $Y_i$ or using a wrong covariance would make the goal false.

A complete development needs: a multivariate Stein or Slepian interpolation, third-order Taylor expansion with remainder for functions on $\mathbb R^p$, Gaussian integration by parts, and Nazarov's inequality. The anti-concentration and smooth-max lemmas are reusable for any max-type Gaussian approximation, including the bootstrap results of the second mission of this series. Proofs of individual milestones, of Nazarov's inequality in particular, are welcome on their own.

Selected references

• V. Chernozhukov, D. Chetverikov, K. Kato, Central limit theorems and bootstrap in high dimensions, Ann. Probab. 45(4):2309–2352, 2017. https://doi.org/10.1214/16-AOP1113
• V. Chernozhukov, D. Chetverikov, K. Kato, Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors, Ann. Statist. 41(6):2786–2819, 2013. https://doi.org/10.1214/13-AOS1161
• V. Bentkus, On the dependence of the Berry–Esseen bound on dimension, J. Statist. Plann. Inference 113:385–402, 2003. MR1965117, https://mathscinet.ams.org/mathscinet-getitem?mr=1965117
• F. Nazarov, On the maximal perimeter of a convex set in ℝⁿ with respect to a Gaussian measure, Geometric Aspects of Functional Analysis, Lecture Notes in Math. 1807, 169–187, Springer, 2003. MR2083397, https://mathscinet.ams.org/mathscinet-getitem?mr=2083397

Motivation

Approximate message passing (AMP) is an iterative method for large systems with random matrices. At each step, every coordinate receives a weighted combination of nonlinear messages from all other coordinates, together with a correction for dependence created by earlier steps. For symmetric random matrices, the correction makes a simple description of the high-dimensional iterates possible: a finite collection of covariance matrices predicts their limiting distributions. Bayati, Lelarge and Montanari prove this state evolution statement for polynomial AMP on a broad class of independent sub-Gaussian matrix entries, rather than only Gaussian entries (paper, §§1.2–1.3). It gives a precise distributional account of an algorithm whose coordinates are coupled at every iteration.

The paper has three distinct headline results. Its Theorem 3 compares AMP moments under two matrix laws with equal entry variances, and its Theorem 2 concerns projected cross-polytopes. This mission targets Theorem 4, which identifies the limit of empirical averages of AMP iterates with a Gaussian expectation. The companion missions address the other results; the statements here use their own copy of the paper's AMP model because draft definitions cannot yet be imported across proposals (paper, pp. 5, 8–9).

Setting

For each size $N$, let $A(N)$ be a random symmetric $N\times N$ matrix with zero diagonal. Its entries above the diagonal are independent, centered, and sub-Gaussian at scale $C/N$. Each coordinate $i$ carries a polynomial map $f^i(\cdot;t):\mathbb R^q\to\mathbb R^q$ at iteration $t$ and an initial vector $x_i^0$. Their coefficients have degree at most $d$ and a uniform bound $C$. The AMP orbit is

$$x_i^{t+1}=\sum_j A_{ij}f^j(x_j^t;t) -\sum_j A_{ij}^{2}Jf^j(x_j^t;t)f^i(x_i^{t-1};t-1).$$

The second sum is the Onsager memory term; it is absent at $t=0$. Here $Jf^j$ is the Jacobian, so its product with the preceding $q$-vector is another $q$-vector. The initial vectors obey an exponential bound on their squared Euclidean norms. These assumptions form the paper's $(C,d)$-regularity condition (Definitions 2–4, pp. 6–7).

A converging sequence has a fixed finite set of classes $a\in[k]$. Coordinate $i$ belongs to $C_a^N$ for one class $a$, whose proportion approaches $c_a\in(0,1)$. A symmetric nonnegative matrix $W$ specifies entry variances: if $i\in C_a^N$ and $j\in C_b^N$, then $\mathbb E A_{ij}^2=W_{ab}/N$. Independent labels $Y(i)$ have class-dependent laws $P_a$, each a finite mixture of Gaussian laws that may be degenerate. The polynomial maps have the common form $f^i(x;t)=g(x,Y(i),a,t)$ for $i\in C_a^N$ (Definition 5, pp. 8–9).

The empirical second moments of $g(x_i^0,Y(i),a,0)$ converge in probability to a prescribed positive semidefinite matrix $\widehat\Sigma_a^0$. For $t\ge1$, state evolution defines

$$\Sigma_a^t=\sum_b c_bW_{ab}\widehat\Sigma_b^{t-1},\qquad \widehat\Sigma_a^t=\mathbb E\big[g(Z_a^t,Y_a,a,t)g(Z_a^t,Y_a,a,t)^{\mathsf T}\big],$$

where $Z_a^t\sim\mathcal N(0,\Sigma_a^t)$ is independent of $Y_a\sim P_a$. The index $t-1$ in the first formula and $t$ in the second are essential (equations (1.9)–(1.11), p. 9).

Formalization targets

Coordinate moments and concentration

For a nonnegative multi-index $m$ and a coordinate sequence eventually in class $a$, Proposition 4 identifies the Gaussian limit of each moment. Proposition 5 says the variance of the class-average monomial tends to zero:

$$\mathbb E[(x_{i(N)}^t)^m]\longrightarrow\mathbb E[(Z_a^t)^m],\qquad \operatorname{Var}\!\left(\frac1{|C_a^N|}\sum_{i\in C_a^N}(x_i^t)^m\right)\longrightarrow0.$$

The diagonal two-time identity (4.4) is another milestone: $\Sigma_a^{t,t}$ has four blocks, each equal to $\Sigma_a^t$. These are the paper's stated intermediate targets (pp. 15, 30, 32).

Main goal: Theorem 4

For $t\ge1$, class $a$, and every locally Lipschitz function $\psi$ with $|\psi(x,y)|\le K(1+\|x\|_2^2+\|y\|_2^2)^K$ for some $K$, prove convergence in probability:

$$\frac1{|C_a^N|}\sum_{j\in C_a^N}\psi(x_j^t,Y(j)) \longrightarrow \mathbb E[\psi(Z_a^t,Y_a)].$$

The target covers polynomial-growth test functions, including functions beyond monomials. The expectation on the right is under the product law of the Gaussian state and the independent class label (Theorem 4, p. 9).

Significance

The theorem assigns a computable asymptotic distribution to each class of AMP coordinates at every fixed positive iteration time. It can be used to evaluate empirical observables without tracking the full coupled $N$-dimensional random trajectory. The class weights and variance profile permit heterogeneous coordinates, while the sub-Gaussian assumption covers more matrix laws than a Gaussian-only statement (paper, Introduction and Theorem 4).

The paper proves this result. The mission's Lean statements are compilation-checked targets, with proofs still to be formalized. A completed development would establish the measure-theoretic meaning of the covariance recursion and its Gaussian expectations, as well as the moment and concentration claims needed for the empirical limit. The polynomial representation, random-matrix regularity predicate, and Gaussian transport are reusable for other AMP results from the same paper.

Difficulty

The matrix $A(N)$ is reused at every AMP iteration, so a new coordinate is correlated with earlier coordinates and with the matrix entries that produced them. A direct independent-summand argument for its Gaussian limit is unavailable. The memory term is tuned to this dependence, and any formal proof must account for its exact Jacobian and squared-entry form. Passing from individual moments to empirical averages also requires control of correlations between different coordinates. The target allows locally Lipschitz functions with polynomial growth, so convergence for a fixed monomial alone does not finish the job (paper, §4.7).

Formalization scope

Vectors in $\mathbb R^q$ are functions Fin q → ℝ; squared Euclidean norms are explicit sums of coordinate squares. Polynomials use bounded real coefficient arrays as in (4.10), with an explicit Jacobian. Fin N uses zero-based indices. The orbit stores the current and preceding iterates and omits the memory term at $t=0$. A single probability space carries all sizes; matrix entries, labels, coefficients, and initial vectors have explicit measurability and independence conditions.

The paper's displayed sub-Gaussian formula conflicts with its stated $C/N$ scale, so the model uses the scale required by its moment estimates and imposes the moment-generating bound for all real arguments. Its initial exponential bound is read almost surely at every size, which ensures the expectations used by the paper are well-defined. Label-dependent polynomial coefficients are measurable. The full label vector is independent of the matrix and initial vectors; this makes explicit the independence needed for Theorem 4's Gaussian-label product law. Independence of the induced polynomial coefficients alone would not imply it when $g$ ignores part of a label. The initial covariance limit is entrywise convergence in probability; finite-dimensional equivalence makes this the same matrix convergence. Class averages divide by the class cardinality, even though that expression has Lean's default value at an empty class for a small $N$. Positive limiting class fractions make those finite exceptions irrelevant.

Gaussian laws live on Mathlib's Euclidean space and are transported to coordinate vectors. Theorems explicitly assert positive semidefiniteness or integrability where a non-PSD covariance, non-integrable expectation, or non-square-integrable variance could otherwise give a default value. The source's $Y(i)$ in (1.12) is read as $Y(j)$, its summation index, and its final covariance as $\Sigma_a^t$. The Gaussian covariance is computed by the recursion above; it is not defined from the proposed limit. Contributions establishing moment bounds, covariance positivity, Gaussian integrability, the three milestones, or the full Theorem 4 are within scope.

Selected references

• M. Bayati, M. Lelarge and A. Montanari, Universality in polytope phase transitions and message passing algorithms, Annals of Applied Probability 25(2), 2015, 753–822. arXiv:1207.7321v2; DOI:10.1214/14-AAP1010.

Motivation

A set-valued predictor returns, for each input $x$, a set of plausible responses rather than a single guess. Such sets are used wherever a model's output feeds a decision with asymmetric costs: a medical classifier that lists every diagnosis it cannot rule out, a segmentation model that marks every pixel that may belong to a tumour, a protein-structure predictor that reports a range of distances. Bates, Angelopoulos, Lei, Malik and Jordan, Distribution-Free, Risk-Controlling Prediction Sets (J. ACM 68(6), 2021), calibrate a nested family of such predictors $\{\mathcal T_\lambda\}$ so that a chosen loss is controlled with high probability, without distributional assumptions.

Calibration says nothing about which family to calibrate. Any nested family can be made risk-controlling; families differ in how large their sets are, and large sets are uninformative. §4 of the paper asks which family is best in the population: which predictor achieves a given risk with the smallest expected set size. For a weighted miscoverage loss (Theorem 7, §4.2) the answer thresholds a weighted conditional density, in the spirit of the optimality results of Sadinle, Lei and Wasserman (2019) for classification. Theorem 8 (§4.3) extends the answer to every loss of integral form. This mission formalizes Theorem 8.

Setting

Let $(X, Y)$ take values in $\mathcal X \times \mathcal Y$ with law $P$, marginal $P_X$, and conditional law $\kappa(x)$ of $Y$ given $X = x$. Prediction sets are subsets of a measurable space $\mathcal Z$ (in the paper's examples $\mathcal Z = \mathcal Y$). A set-valued predictor is a map $\mathcal T : \mathcal X \to 2^{\mathcal Z}$.

Fix a finite measure $\mu$ on $\mathcal Z$ and a nonnegative cost $\ell : \mathcal Y \times \mathcal Z \to [0, \infty)$; $\ell(y, z)$ is the cost of leaving $z$ out of the prediction set when the true response is $y$. The loss of a set $\mathcal S$ is the total cost of what it leaves out,

$$L(y; \mathcal S) = \int_{z \in \mathcal S^c} \ell(y, z)\, d\mu(z),$$

and the risk of a predictor is $R(\mathcal T) = \mathbb E[L(Y; \mathcal T(X))]$ (§2.1). Larger sets have smaller loss. The conditional expected cost of omitting $z$ at $x$ is

$$\mathbb E[\ell(Y; z) \mid X = x] = \int \ell(y, z)\, d\kappa(x)(y),$$

and for $\lambda < 0$ the threshold predictor of display (11) is

$$\mathcal T_\lambda(x) = \{z \in \mathcal Z : \mathbb E[\ell(Y; z) \mid X = x] \ge -\lambda\}.$$

It includes exactly the points whose omission is expected to cost at least $-\lambda$. The size of a set is its $\mu$-measure, $|\mathcal S| = \mu(\mathcal S)$, and the expected size of a predictor is $\mathbb E[|\mathcal T(X)|] = \int \mu(\mathcal T(x))\, dP_X(x)$.

In Lean: setLoss μ ℓ y S, risk P μ ℓ T, condLoss κ ℓ x z, thresholdSet κ ℓ lam x and expSize P μ T in the namespace RiskControl.Optimal.

Formalization targets

Goal: Theorem 8

For every $\lambda < 0$ and every set-valued predictor $\mathcal T'$ (with measurable graph),

$$R(\mathcal T') \le R(\mathcal T_\lambda) \implies \mathbb E\big[|\mathcal T_\lambda(X)|\big] \le \mathbb E\big[|\mathcal T'(X)|\big].$$

No constant is involved; the statement is uniform over $\ell$, $\mu$, $P$ and $\mathcal T'$.

Milestones: the steps of the proof (p. 27)

1. Moving the conditional expectation inside the $\mu$-integral: for a predictor $\mathcal A$ with measurable graph, $\mathbb E[\int_{\mathcal A(X)} \ell(Y; z)\, d\mu(z)] = \mathbb E[\int_{\mathcal A(X)} \mathbb E[\ell(Y; z) \mid X]\, d\mu(z)]$.
2. From the risk comparison to the set differences:

$$\mathbb E\Big[\int_{\mathcal T'(X) \setminus \mathcal T_\lambda(X)} \mathbb E[\ell(Y; z) \mid X]\, d\mu\Big] \ge \mathbb E\Big[\int_{\mathcal T_\lambda(X) \setminus \mathcal T'(X)} \mathbb E[\ell(Y; z) \mid X]\, d\mu\Big].$$

3. From costs to sizes: $\mathbb E|\mathcal T'(X) \setminus \mathcal T_\lambda(X)| \ge \mathbb E|\mathcal T_\lambda(X) \setminus \mathcal T'(X)|$.
4. From set differences to sets: $\mathbb E|\mathcal A(X) \setminus \mathcal B(X)| \le \mathbb E|\mathcal B(X) \setminus \mathcal A(X)|$ implies $\mathbb E|\mathcal A(X)| \le \mathbb E|\mathcal B(X)|$.

Significance

The result. Theorem 8 identifies the population-optimal shape of a set-valued predictor for every loss that charges for omitted points through a nonnegative cost and a finite measure. It justifies building nested families by thresholding an estimate of the conditional expected cost: if the estimate were exact, no other predictor with the same or smaller risk could have smaller sets on average. It contains the weighted-miscoverage case (Theorem 7) when $\mathcal Y$ is discrete and $\mu$ is counting measure, and it covers losses for which no such special case exists, such as costs that depend on the distance between $z$ and $y$.

Formalizing it. The proof on p. 27 is a chain of eight displays. A machine-checked version settles the measure-theoretic points the chain passes over: the conditional expectation inside an integral over a random set, the cancellation of possibly infinite terms, and the range of $\lambda$. To our knowledge no formal proof of Theorem 8 or of the Neyman–Pearson-type optimality results it generalizes for set predictors exists in Lean or Mathlib.

Difficulty

The natural reading of the proof subtracts the quantity $\mathbb E\int_{\mathcal Z} \mathbb E[\ell(Y; z) \mid X]\, d\mu(z)$ from both sides of the risk comparison to pass from complements to sets. Nothing in the setting makes this quantity finite, so that step is not valid as written, and a proof that follows the displays literally breaks there. The theorem itself does not need finiteness of $\ell$, and no such hypothesis is in scope.

The other obstacle is measure-theoretic bookkeeping. The conditional expectation $\mathbb E[\ell(Y; z) \mid X]$ is indexed by a point $z$ of a second space, and must be a jointly measurable function of $(x, z)$ before it can be integrated over the random set $\mathcal T(X)$; here it is given by a regular conditional distribution, and every integral over $\mathcal T(x)$ needs the measurability of the graph of $\mathcal T$.

Formalization scope

• Spaces and measures. $\mathcal X, \mathcal Y, \mathcal Z$ are arbitrary measurable spaces; $P$ is a probability measure on $\mathcal X \times \mathcal Y$; $\mu$ is a finite measure on $\mathcal Z$; $\ell$ takes values in $[0, \infty)$ and is jointly measurable.
• Conditional law. $\mathbb E[\ell(Y; z) \mid X = x]$ is $\int \ell(y, z)\, d\kappa(x)(y)$ for a Markov kernel $\kappa$ with $P = P_X \otimes \kappa$, which is exactly the statement that $\kappa$ is a regular conditional distribution of $Y$ given $X$. Such a kernel exists, for instance, when $\mathcal Y$ is standard Borel; the formalization does not assume this. The threshold predictor is built from the same $\kappa$ and the same $\ell$ as the risk; neither is a free object.
• Integrals. Every integral is a lower Lebesgue integral with values in $[0, \infty]$; no integrability hypothesis is made, and an infinite risk is $+\infty$.
• Range of $\lambda$ (correction of the print). The paper states Theorem 8 for $\lambda \in \Lambda \subset (-\infty, 0]$. At $\lambda = 0$ it is false: with $\ell \equiv 0$ and $\mu(\mathcal Z) > 0$, $\mathcal T_0(x) = \mathcal Z$ and $\mathcal T' \equiv \emptyset$ both have risk $0$, but $\mathbb E|\mathcal T'(X)| = 0 < \mu(\mathcal Z)$. The proof divides by $-\lambda$. Theorem 8 and its milestones are posed for $\lambda < 0$.
• Set size (correction of the print). p. 13 describes $|\cdot|$ as Lebesgue or counting measure, a sentence written for Theorem 7. In Theorem 8 the proof measures sets with the measure $\mu$ of the loss, and $|\cdot| = \mu$ throughout.
• Measurable graphs. The competing predictor $\mathcal T'$ is required to have a measurable graph $\{(x, z) : z \in \mathcal T'(x)\}$; otherwise its risk and expected size are not integrals of measurable functions. The page leaves this implicit. The graph of $\mathcal T_\lambda$ is measurable by construction.
• Not a trivialization. The goal is stated about the threshold sets (11) built from the given $\ell$, $\mu$ and $\kappa$, for every competing predictor; it does not assume any of the intermediate inequalities, and its hypotheses are satisfiable (for instance $\ell \equiv 1$, $\lambda = -1/2$, where $\mathcal T_\lambda(x) = \mathcal Z$).

The development needs Mathlib's composition of a measure with a kernel, Tonelli's theorem for kernels and lower integrals over sections of a measurable set. Lemmas about expected $\mu$-sizes of random sets (milestone 4) and about moving a kernel integral inside a set integral (milestone 1) are reusable beyond this mission. Proofs of the milestones, and of the goal by other routes, are welcome.

Selected references

• S. Bates, A. Angelopoulos, L. Lei, J. Malik, M. I. Jordan, Distribution-Free, Risk-Controlling Prediction Sets, Journal of the ACM 68(6), 2021. arXiv:2101.02703v3. https://arxiv.org/abs/2101.02703 — Theorem 8, p. 13; proof p. 27.
• M. Sadinle, J. Lei, L. Wasserman, Least Ambiguous Set-Valued Classifiers With Bounded Error Levels, Journal of the American Statistical Association 114(525), 2019. https://doi.org/10.1080/01621459.2017.1395341

Why ask whether a Wasserstein distance is Hilbertian

Kernel methods, multidimensional scaling and low-distortion embeddings all work best on data whose distance comes from a Hilbert space. A distance $d$ on a set $\mathcal Z$ is Hilbertian when $\mathcal Z$ can be mapped into a Hilbert space so that $d$ becomes the norm distance; for such a $d$, the kernel $e^{-d^p/t}$ is positive definite for $0\le p\le 2$ and $t>0$ (Berg, Christensen and Ressel, 1984), and the classical Euclidean toolbox applies. Optimal transport distances are popular for comparing probability measures, so it is natural to ask whether they are Hilbertian. In dimension one the answer is yes for $\mathcal W_2$: the map sending a measure to its quantile function is an isometry into $L^2([0,1])$, and between univariate Gaussians $\mathcal W_2$ is the Euclidean distance between (mean, standard deviation). Peyré and Cuturi, Computational Optimal Transport (FnT ML 2019), §8.3, show that this does not extend to the plane: for $p=1,2$ and $d\ge2$ the $p$-Wasserstein distance on $\mathbb R^d$ is not Hilbertian (Proposition 8.2, p. 507).

The tool is a classical characterization. Schoenberg (1938) proved that a distance is Hilbertian exactly when its square is conditionally negative definite; Berg, Christensen and Ressel (1984, Prop. 3.2) give the modern form. The easy direction turns non-embeddability into a finite check, and the book's proof of Proposition 8.2 performs that check numerically on 35 measures supported on the four corners of the unit square (Figure 8.6, p. 508). Stronger, quantitative statements are known: planar $\mathcal W_1$ does not even embed into $L^1$ with bounded distortion (Naor and Schechtman, 2007), and Andoni, Naor and Neiman (2018) study which powers of Wasserstein distances embed.

Setting

Let $\mathcal Z$ be a set. A function $\varphi:\mathcal Z\times\mathcal Z\to\mathbb R$ is negative definite (Definition 8.3, p. 501, in the conditional form used in §8.3) if it is symmetric and, for every $n\ge0$, every $x_1,\dots,x_n\in\mathcal Z$ and every $r\in\mathbb R^n$ with $\sum_ir_i=0$,

$$\sum_{i,j=1}^n r_ir_j\,\varphi(x_i,x_j)\le0 .$$

A function $d:\mathcal Z\times\mathcal Z\to\mathbb R$ is Hilbertian (Definition 8.4, p. 506) if there are a real Hilbert space $\mathcal H$ and a map $\phi:\mathcal Z\to\mathcal H$ with $d(z,z')=\|\phi(z)-\phi(z')\|_{\mathcal H}$ for all $z,z'$.

On $\mathbb R^d$ with the Euclidean norm $\|\cdot\|_2$, let $\mathcal P_p(\mathbb R^d)$ be the Borel probability measures $\mu$ with $\int\|x\|_2^p\,d\mu<\infty$. A coupling of $\mu$ and $\nu$ is a measure $\pi$ on $\mathbb R^d\times\mathbb R^d$ whose two marginals are $\mu$ and $\nu$, and the $p$-Wasserstein distance (2.18) is

$$\mathcal W_p(\mu,\nu)=\Big(\inf_{\pi}\int\|x-y\|_2^p\,d\pi(x,y)\Big)^{1/p},$$

finite on $\mathcal P_p(\mathbb R^d)$. For points $x_1,\dots,x_n$ and a histogram $a\in\Sigma_n$ (nonnegative, summing to one) the discrete measure is $\sum_ia_i\delta_{x_i}$; for two histograms $a,b$, $U(a,b)$ is the set of nonnegative matrices with row sums $a$ and column sums $b$, and $L_C(a,b)=\min_{P\in U(a,b)}\sum_{i,j}C_{i,j}P_{i,j}$ (2.11).

The configuration of the book's proof consists of the corners $x^1=[0,0]$, $x^2=[1,0]$, $x^3=[0,1]$, $x^4=[1,1]$ of the unit square, the 35 histograms of $\Sigma_4$ with entries in $\{0,\tfrac14,\tfrac12,\tfrac34,1\}$, and the centering matrix $J=I_n-\tfrac1n\mathbb 1_{n,n}$.

Formalization targets

Goal: Proposition 8.2

For $d\ge2$ and $p\in\{1,2\}$,

$$\mathcal W_p \text{ on } \mathcal P_p(\mathbb R^d) \text{ is not Hilbertian.}$$

Milestones

1. Proposition 8.1, proved direction (pp. 506–507): if $d$ is Hilbertian then $d^2$ is negative definite.
2. Squared distance (p. 507): $\mathcal W_p^2$ on $\mathcal P_p(\mathbb R^d)$ is not negative definite, $d\ge2$, $p=1,2$.
3. Reduction to the plane (proof of Prop. 8.2): if $\mathcal W_p^2$ on $\mathcal P_p(\mathbb R^d)$, $d\ge2$, were negative definite, so would be $\mathcal W_p^2$ on $\mathcal P_p(\mathbb R^2)$.
4. Discrete measures (Remark 2.13, p. 375): $\mathcal W_p\big(\sum_ia_i\delta_{x_i},\sum_jb_j\delta_{x_j}\big)^p=L_C(a,b)$ with $C_{i,j}=\|x_i-x_j\|_2^p$.
5. Centering criterion (proof of Prop. 8.2): a matrix $M$ admits a zero-sum $r$ with $r^\top Mr>0$ if and only if $JMJ$ has a positive eigenvalue.
6. The grid counterexample (proof of Prop. 8.2, Figure 8.6): for $p=1,2$ there are grid histograms $a^1,\dots,a^n$ and a zero-sum $r$ with

$$\sum_{i,j}r_ir_j\,\mathcal W_p^2\Big(\sum_k a^i_k\delta_{x^k},\sum_k a^j_k\delta_{x^k}\Big)>0.$$

Significance

A Hilbertian distance gives positive definite Gaussian and Laplace kernels, Euclidean multidimensional scaling and Johnson–Lindenstrauss dimension reduction for free. Proposition 8.2 says none of these can be obtained for $\mathcal W_1$ or $\mathcal W_2$ on $\mathbb R^d$, $d\ge2$, by an isometric embedding; this is why the literature turned to approximate embeddings, sliced distances and entropic surrogates (pp. 507–508). The negative-definiteness test of milestone 1 is a reusable certificate for non-embeddability of any finite metric.

The result is known, and the book's proof is a floating-point eigenvalue computation. The formal content of this mission is a machine-checked proof with an exact certificate: rational or algebraic Wasserstein distances between explicit discrete measures and an explicit zero-sum vector. As far as this mission is aware, no formal proof of the statement exists in Mathlib or on the platform. Schoenberg's converse is cited by the book and is not part of the mission.

Difficulty

The abstract part is short: expanding $\|\phi(z_i)-\phi(z_j)\|^2$ and using $\sum_ir_i=0$ gives milestone 1. The obstacle is the certificate. The book observes that $J\mathbf D_p^2J$ has a positive top eigenvalue (about $1.2$ for $p=1$ and $0.7$ for $p=2$), but a proof needs exact transport costs between the chosen measures, each the value of a small linear program whose optimality has to be established (for instance by a dual certificate), together with a vector $r$ for which the quadratic form is provably positive. For $p=1$ the costs involve $\sqrt2$, the diagonal of the square. A second, separate difficulty is connecting the measure-theoretic $\mathcal W_p$ to the discrete program (milestone 4) and transporting a counterexample from $\mathbb R^2$ into $\mathbb R^d$ (milestone 3), both of which require working with couplings as measures on a product space.

Formalization scope

• $\mathbb R^d$ is EuclideanSpace ℝ (Fin d); $\mathcal W_p$ is the published definition WassersteinDRO.Duality.wassersteinDistance (an $[0,\infty]$-valued infimum over all measures with the two marginals, raised to $1/p$), converted to a real number. The distance is considered on $\mathcal P_p(\mathbb R^d)$, the set where it is finite; outside it the real conversion is a junk $0$. The book names no domain; since the counterexample consists of finitely supported measures, the finite computation gives non-Hilbertianity on $\mathcal P_p$ and on every smaller set containing those measures.
• Negative definiteness is conditional. Definition 8.3 as printed quantifies over all $r$; under that reading no nonzero squared distance is negative definite, so Proposition 8.1 would be false and "not negative definite" trivially true. The zero-sum condition, used by the book in the proof of Proposition 8.1, is part of the definition; this rules out the trivial formalization.
• Hilbert spaces are real, complete inner-product spaces in the universe of $\mathcal Z$. Indices are 0-based (Fin n). Histograms are real vectors; the reduction milestone is stated for every $p\ge1$.
• Discrete transport objects ($U(a,b)$, $L_C$) are redefined in this mission's namespace CompOT.NotHilbertian, duplicating chunk II of the series.

Contributions welcome: the forward direction of Schoenberg's characterization, the discrete-measure bridge (reusable wherever discrete optimal transport meets measure-theoretic couplings), the isometric-embedding invariance of $\mathcal W_p$, and exact certificates for the grid computation.

Selected references

• G. Peyré and M. Cuturi, Computational Optimal Transport, Foundations and Trends in Machine Learning 11(5–6):355–607, 2019. https://doi.org/10.1561/2200000073
• I. J. Schoenberg, Metric spaces and positive definite functions, Transactions of the AMS 44(3):522–536, 1938. https://doi.org/10.1090/S0002-9947-1938-1501980-0
• C. Berg, J. P. R. Christensen and P. Ressel, Harmonic Analysis on Semigroups, Graduate Texts in Mathematics 100, Springer, 1984. https://doi.org/10.1007/978-1-4612-1128-0
• A. Naor and G. Schechtman, Planar earthmover is not in $L_1$, SIAM Journal on Computing 37(3):804–826, 2007. https://doi.org/10.1137/05064206X
• A. Andoni, A. Naor and O. Neiman, Snowflake universality of Wasserstein spaces, Annales scientifiques de l'École normale supérieure 51(3):657–700, 2018. https://arxiv.org/abs/1509.08677

Motivation

Best-subset selection with ridge shrinkage, $\min_\beta \tfrac12\|y-X\beta\|_2^2+\lambda_0\|\beta\|_0+\lambda_2\|\beta\|_2^2$, is a mixed-integer program that statisticians want to solve to certified optimality at the scale of modern data ($p$ up to $10^7$ features). Hazimeh, Mazumder and Saab (arXiv:2004.06152, Mathematical Programming 2022) built a branch-and-bound solver, L0BnB, whose node relaxations are solved in the primal space by active-set coordinate descent instead of by an interior-point method. Branch-and-bound prunes a node only with a valid dual bound, a certified lower bound on the node's relaxation value. A primal method returns an approximate minimizer, not such a bound, so the solver must turn an inexact primal point into a dual feasible point and needs to know how much is lost in doing so. This mission formalizes the paper's answer, its main theorem (Theorem 3).

Setting

Data are $X\in\mathbb R^{n\times p}$ with columns $X_1,\dots,X_p$, $y\in\mathbb R^n$, and parameters $\lambda_0,\lambda_2,M>0$. The reverse Huber penalty is $\mathcal B(t)=|t|$ for $|t|\le1$ and $(t^2+1)/2$ for $|t|\ge1$. The penalty $\psi(b)$ equals $\psi_1(b)=2\lambda_0\mathcal B(b\sqrt{\lambda_2/\lambda_0})$ if $\sqrt{\lambda_0/\lambda_2}\le M$ and $\psi_2(b)=(\lambda_0/M+\lambda_2M)|b|$ otherwise. The reduced relaxation (5) is

$$\min_{\beta\in\mathbb R^p}\ F(\beta)=\tfrac12\|y-X\beta\|_2^2+\sum_{i}\psi(\beta_i)\quad\text{s.t.}\quad \|\beta\|_\infty\le M .$$

Throughout, the columns of $X$ and $y$ have unit $\ell_2$ norm (the standing assumption of Section 3 of the paper).

Algorithm 2 (active-set coordinate descent) returns a point $\hat\beta$ in the box such that the set

$$V=\{i\notin\operatorname{Supp}(\hat\beta):\ 0\ne\operatorname*{arg\,min}_{|t|\le M}F(\hat\beta_1,\dots,t,\dots,\hat\beta_p)\}$$

is empty: no coordinate outside the support wants to move. Write $\hat r=y-X\hat\beta$, $k=\|\hat\beta\|_0$, and let $\beta^*$ be an optimal solution of (5) with $r^*=y-X\beta^*$. The primal gap is $\epsilon=\|X(\beta^*-\hat\beta)\|_2$.

The two duals of (5) (Theorem 2 of the paper) are, with $v(\alpha,\gamma_i)=[(\alpha^\top X_i-\gamma_i)^2/(4\lambda_2)-\lambda_0]_++M|\gamma_i|$,

$$h_1(\alpha,\gamma)=-\tfrac12\|\alpha\|_2^2-\alpha^\top y-\textstyle\sum_i v(\alpha,\gamma_i),\qquad h_2(\rho,\mu)=-\tfrac12\|\rho\|_2^2-\rho^\top y-M\|\mu\|_1,$$

the latter under $|\rho^\top X_i|-\mu_i\le\lambda_0/M+\lambda_2M$. The dual variables built from $\beta^*$ are $\alpha^*=\rho^*=-r^*$, $\gamma^*_i=\mathbb 1_{[|\beta^*_i|=M]}(\alpha^{*\top}X_i-2M\lambda_2\operatorname{sign}(\alpha^{*\top}X_i))$ and $\mu^*_i=\mathbb 1_{[|\beta^*_i|=M]}(|\rho^{*\top}X_i|-\lambda_0/M-\lambda_2M)$. The dual points built from $\hat\beta$ are $\hat\alpha=\hat\rho=-\hat r$, with $\hat\gamma$ a maximizer of $h_1(\hat\alpha,\cdot)$ (25) and $\hat\mu$ a maximizer of $h_2(\hat\rho,\cdot)$ under the constraints (27).

Formalization targets

Goal: Theorem 3 with the proof's constants

If $\sqrt{\lambda_0/\lambda_2}\le M$, with $c_i=(2\lambda_2)^{-1}$ when $|\beta^*_i|<M$ and $c_i=M$ when $|\beta^*_i|=M$,

$$h_1(\hat\alpha,\hat\gamma)\ \ge\ h_1(\alpha^*,\gamma^*)-2\epsilon-\tfrac12\epsilon^2-\sum_{i\in\operatorname{Supp}(\hat\beta)}\big(c_i\epsilon+(4\lambda_2)^{-1}\epsilon^2\big).\qquad(55)$$

If $\sqrt{\lambda_0/\lambda_2}>M$,

$$h_2(\hat\rho,\hat\mu)\ \ge\ h_2(\rho^*,\mu^*)-\epsilon(2+Mk)-\tfrac12\epsilon^2.\qquad(59)$$

The paper states these as $-kO(\epsilon)-kO(\epsilon^2)$ (29) and $-kO(\epsilon)-O(\epsilon^2)$ (30); the goal states the expressions its proof establishes.

Milestones

In the order the proof uses them: the closed-form coordinate updates (14) and (15); Proposition 3, $V=\{i\notin\operatorname{Supp}(\hat\beta): |\langle\hat r,X_i\rangle|>c(\lambda_0,\lambda_2,M)\}$; the bound $\|\alpha^*\|_2\le1$; Lemma 2, which bounds $v(\hat\alpha,\hat\gamma_i)$ coordinate by coordinate and makes it vanish off the support; inequality (53); the closed form (28) of $\hat\mu$; the optimality conditions (58); the bound (57) $|\hat\mu_i|\le\epsilon+|\mu^*_i|$ on the support; and inequality (56).

Significance

The result. In both regimes the loss of the dual bound is controlled by the primal gap times the sparsity $k$ of the iterate, with constants depending only on $M$ and $\lambda_2$; the number of features $p$ does not enter. Since L0BnB seeks solutions with $k\ll p$, the cheap dual bound obtained from an inexact primal solution is nearly as good as the exact one, which is what allows pruning without an interior-point solve at every node. The paper notes that with plain coordinate descent in place of Algorithm 2 the same argument gives $p$ in place of $k$.

Formalizing it. The result is proved in the paper; it has no machine-checked proof. The mission produces a checked version of the main theorem with explicit constants, a formal model of what "output of an active-set method" means for the analysis, and reusable closed forms for the boxed soft-thresholding updates of an $\ell_1$- or reverse-Huber-penalized box-constrained least squares problem.

Difficulty

The bound is not a consequence of weak duality alone: weak duality says only that each dual value is below the primal optimum, and says nothing about how far the constructed dual point is from the dual optimum. The obvious estimate, a Lipschitz bound on $h_1$ or $h_2$ summed over all coordinates, gives a loss proportional to $p$. Getting $k$ instead requires showing that the dual contribution of every coordinate outside $\operatorname{Supp}(\hat\beta)$ vanishes exactly, which uses the emptiness of $V$ through Proposition 3 and the closed forms (14)–(15) of one-dimensional nonsmooth box-constrained problems. On the support, the loss must be bounded using the structure of $\gamma^*$ and $\mu^*$ at coordinates where $\beta^*$ hits the box, which is where $c_i=M$ and (58) enter.

Formalization scope

Indices are Fin p; norms are explicit sums; $\operatorname{sign}$ is Real.sign with $\operatorname{sign}(0)=0$; $[a]_+$ is max a 0. The unit norms of the columns of $X$ and of $y$ are hypotheses, never built into the definitions, and each theorem assumes only the ones it uses. $\beta^*$ is any point of the box minimizing $F$ over the box. The output of Algorithm 2 is modelled by its two properties used in the paper: $\|\hat\beta\|_\infty\le M$ and $V=\emptyset$ ("$0\ne\arg\min$" encoded as "$0$ is not a minimizer"); the iterations are not formalized. $\hat\gamma$ and $\hat\mu$ are arbitrary maximizers satisfying (25) and (27). The dual variables (23)–(24) are defined by their formulas; their optimality (Theorem 2, a separate mission) is neither assumed nor needed, and the goal is an inequality between explicit numbers.

Instantiations of $O(\cdot)$. (29) is stated as (55): $kO(\epsilon)+kO(\epsilon^2)$ becomes $2\epsilon+\tfrac12\epsilon^2+\sum_{i\in\operatorname{Supp}(\hat\beta)}(c_i\epsilon+(4\lambda_2)^{-1}\epsilon^2)$ with $c_i\in\{(2\lambda_2)^{-1},M\}$ decided by $\beta^*$. (30) is stated as (59): $kO(\epsilon)+O(\epsilon^2)$ becomes $\epsilon(2+Mk)+\tfrac12\epsilon^2$. The proof's final rearrangement of (29), which counts $|\hat\beta_i|$ instead of $|\beta^*_i|$, is not used.

A formalization that takes $\hat\beta=\beta^*$, $\hat\gamma=\gamma^*$ or $\hat\mu=\mu^*$, that drops the hypothesis $V=\emptyset$, or that replaces the constants by an existential "$\exists C$" would be trivial or a different statement; the goal quantifies over every optimal $\beta^*$, every box-feasible $\hat\beta$ with $V=\emptyset$ and every maximizer $\hat\gamma$, $\hat\mu$, with the explicit constants above.

Needed infrastructure: one-dimensional convex minimization on an interval with piecewise penalties, Cauchy–Schwarz for finite sums, and first-order optimality conditions for box-constrained composite problems. Proofs of any milestone are welcome independently.

Selected references

• H. Hazimeh, R. Mazumder, A. Saab, Sparse Regression at Scale: Branch-and-Bound rooted in First-Order Optimization, arXiv:2004.06152v2, 2021; Mathematical Programming 196 (2022). https://arxiv.org/abs/2004.06152
• P. Tseng, Convergence of a block coordinate descent method for nondifferentiable minimization, J. Optim. Theory Appl. 109 (2001). https://doi.org/10.1023/A:1017501703105
• J. Friedman, T. Hastie, R. Tibshirani, Regularization paths for generalized linear models via coordinate descent, J. Stat. Softw. 33 (2010). https://doi.org/10.18637/jss.v033.i01