Variational Inference in high-dimensional linear regression

Introduction

In this age of big-data, statisticians routinely analyze large, high-dimensional datasets arising from applications in genomics, finance, public policy etc., with the goal of discovering relationships between the response variable, and the observed features. The linear regression model

is arguably the most common framework for this task when the response is continuous. Under a Bayesian formalism, the statistician posits a prior distribution $\pi$ for the coefficient vector $\bm{\beta}$, and constructs the corresponding posterior. Subsequent inference is based solely on this posterior distribution. Two questions are naturally relevant in this setting:

What are the statistical properties of Bayesian procedures, particularly in high-dimensions?

Are these procedures computationally tractable for large datasets with high-dimensional features?

In general, computationally tractable Bayesian inference for high-dimensional models presents significant challenges. MCMC based strategies have been explored extensively for this purpose. Despite rapid progress in MCMC methodology and supporting theory, these methods are still slower than competing frequentist methodology e.g. those based on convex optimization. Variational methods [WJ08] provide an attractive general option to the Bayesian statistician. The simplest form of variational inference approximates the true posterior distribution using a product distribution—this version is often referred to as naive mean-field Variational Bayes (nVB). Computing the best approximating product distribution is computationally fast, and thus these methods provide a practical option for modern datasets. We refer the interested reader to [BKM17] for an introduction to variational Bayes methods in statistics and machine learning.

Although variational methods provide a computationally feasible strategy, supporting theoretical evidence has been relatively scarce. In [WT06, WB19b], the authors established the correctness of this approach in parametric models in the classical fixed $p$ setting. Subsequently [WB19a] study variational Bayes in misspecified models. In the context of the linear model, early work by [NOW14, OYM17] focused on variational inference in the low dimensional linear model, while [CS12] provides an early variational approximation for variable selection.

The past two years have witnessed rapid progress in the analysis of variational methods for high-dimensional models [ARC16, AR20, HY19, RS21, RSC20, YPB20]. These results focus on high-dimensional models with a sparse underlying truth $\bm{\beta}_{0}$, and study the contraction properties of the variational posterior. In particular, they derive sufficient conditions for the variational posterior to contract at the minimax optimal rate. Correctness of variational methods have also been established for community detection [BCCZ13, ZZ20], a poisson mixed model [HPWW11], frequentist models [WM19] and mixture models [CAA18].

There is thus an immediate need to understand general properties of statistical problems which ensure the correctness of the nVB approximation. In this paper, we study the accuracy of this approximation in Bayesian linear regression with product priors. We provide easily verifiable conditions which ensure the asymptotic accuracy of the nVB approximation. Further, we illustrate that the approximation can yield detailed information regarding the statistical properties of this model, which might be unavailable using other techniques. We elaborate on our specific contributions below.

Our main contributions in this paper are as follows:

In Theorem 1, we provide sufficient conditions for the asymptotic tightness of the nVB approximation. The conditions are easily verifiable, and can be explicitly checked in specific applications. This provides rigorous theoretical support for widely used mean-field approximation based methodology.

Assuming a true frequentist linear model for the data, we derive a limiting variational formula for the log normalizing constant in Theorem 2. We emphasize that in contrast to existing results, we do not assume any sparsity on the true regression coefficients. We refer the interested reader to [LPM20] for a motivating discussion on the importance of such scenarios in scientific applications.

Under an additional “separation" condition (11), we establish that the limiting variational problem has a unique optimizer (see Theorem 3 for a precise statement). Further, this optimizer governs the probabilistic properties of the posterior distribution (Corollary 4). We also provide interpretable sufficient conditions which enforce this separation condition (Lemma 5).

We emphasize that in existing analyses of high-dimensional Bayesian linear regression, properties of the posterior are directly established [CSHVdV15], and the variational posterior is analyzed independently [RS21]. Our approach is inherently different—we first establish the correctness of the mean-field approximation, and then study the posterior through the lens of the mean-field approximation formula.

We further illustrate our general results by applying them to three specific examples—a two factor ANOVA model, a gaussian design setting with spiked covariance, and a sparse bernoulli design. In each case, we identify the specific limiting functional which determines the limiting log normalizing constant.

From a theoretical perspective, our results crucially utilize recent advances in the theory of non-linear large deviations. Initiated in the seminal paper [CD16], the theory of non-linear large deviations was originally conceived to answer some deep questions concerning large deviations of sub-graph counts in sparse random graphs. In [BM17], one of the authors successfully utilized this framework to establish the tightness of the naive mean field approximation for the log-partition function in a family of Potts models. To the best of our knowledge, these developments have not been utilized previously for other statistical models. In this paper, we demonstrate the usefulness of these tools in the context of high dimensional statistics; we hope that this spurs an in-depth study of their applicability for other high dimensional problems. We consider this to be a key conceptual contribution of this paper.

2 Non-linear large deviations and related results

In a breakthrough paper, Chatterjee and Dembo [CD16] introduced the theory of non-linear large deviations with the goal of studying large deviations for non-linear functions of bernoulli random variables. As an application of this general machinery, they characterized sharp deviation probabilities for sub-graph counts in sparse Erdős-Rényi random graphs. Subsequent extensions by Eldan [Eld18] and Augeri [Aug19] allow one to track a wider regime of sparsity for the binary variables. In a different direction, Yan [Yan20] extended the Chatterjee-Dembo framework to general bounded Banach-space valued variables. See also [Aus19] for related decompositions of general Gibbs measures using information theoretic ideas. These results have galvanized the study of large deviations for sub-graph counts on sparse random graphs, and the past three years have witnessed rapid progress in this direction. We refer the interested reader to [CDP21] and references therein for a survey of recent progress in this area.

At the heart of the Chatterjee-Dembo framework lies a tight approximation bound for the log-normalizing constant for general Gibbs type distributions in terms of the naive mean-field approximation formula. This framework was utilized by one of the authors [BM17] to derive asymptotic limits for the log normalizing constant of Potts models on several sequences of graphs. Similar results were independently derived by [JKM18, JKR19] using different techniques.

In this paper, we study Bayesian linear regression through the lens of non-linear large deviations, and uncover precise statistical properties of these models by analyzing the mean-field variational problem.

3 Setup

In the display above, $\pi$ is a probability distribution supported on $$, and we consider an iid prior on the regression coefficients. For our subsequent discussion, the precise interval$$ for the prior support is not crucial— our arguments go through unchanged, as long as the prior has bounded support. Note that in the context of Bayesian inference for linear regression, the regression parameters are often drawn from a parametric family, and the parameters specifying the prior are, in turn, sampled from a hyper-prior. In our discussions, we will restrict ourselves to the simpler setting where the prior $\pi$ is fixed—however, we do not make any assumptions on $\pi$ apart from the bounded support assumption. Throughout, we assume that the noise variance $\sigma^{2}>0$ is fixed and known to the statistician.

Given the Bayesian model (1), one naturally constructs the posterior distribution

Using this definition we can write the posterior distribution $\mu$ as

A central object in the theoretical study of these posterior distributions is the normalizing constant, also referred to as the “partition function" in statistical physics parlance. Formally, we define

The partition function is intractable for most priors, unless special conjugacy properties are satisfied between the prior and the likelihood. Henceforth, we suppress the dependence of $\mu$, $Z_{p}$ on $\mathbf{y}$, $\mathbf{X}$ whenever it is clear from the context. The classical Gibbs variational principle characterizes the partition function as

where the supremum ranges over probability distributions $Q$ on $^{p}$ (see e.g. [WJ08]). In fact, the supremum in the above variational problem is attained by $Q=\mu$. The naive mean-field approximation to $Z_{p}$ restricts the supremum to product distributions, and thus obtains a universal lower bound. Formally, we have,

The optimizer in the display above provides the best approximation to $\mu$ among product distributions under KL divergence. This paper focuses on the tightness of the naive mean-field lower bound in the context of linear regression. For an in-depth survey of variational inference, we refer the interested reader to [Bis06, BKM17, WJ08].

Outline: The rest of the paper is structured as follows. We collect our results in Section 2. We discuss some directions for future enquiry in Section 3. The main results are established in Sections 4. We defer some proofs to the Appendix.

Acknowledgments: The authors thank Pragya Sur for discussions on high-dimensional regression. SM gratefully thanks NSF (DMS 1712037) for support during this research.

Results

We collect our main results in this section. To this end, we first discuss some elementary facts regarding exponential families. The next result collects some analytic properties of the cumulant generating function $c(\cdot)$, which will be relevant for our subsequent discussion. For the sake of completeness, we provide a proof in Section 4.

Armed with these basic facts, we can formally state the naive-mean field approximation to the log-normalizing constant (2).

where $\pi_{\infty}$ and $\pi_{-\infty}$ are degenerate distributions which puts mass 1 at $1$ and $-1$ respectively.

We will need the following facts about the derivatives of $G$. We defer the proof of Lemma 2 to the Appendix.

This result follows immediately from elementary facts about exponential families. We refer the interested reader to [WJ08, Section 5.3].

Throughout the paper, we use the usual Landau notation $O(\cdot)$ and $o(\cdot)$ for deterministic sequences dependent on $p$.

Assume that the matrix $A_{p}$ satisfies the two conditions

We have, setting $R_{p}:=\sup_{\mathbf{u}\in^{p}}M_{p}({\bf{u}})$, as $p\to\infty$,

Suppose there exists $\hat{\bf u}\in^{p}$ which satisfies that for every $\eta>0$ we have

The careful reader would have already noticed that Theorem 1 is stated for deterministic data $(\mathbf{y},\mathbf{X})$. In practical applications, it is often more natural to assume that the data $(\mathbf{y},\mathbf{X})$ is, in turn, sampled from some underlying distribution $\mathcal{P}$. The conclusions of Theorem 1 continue to hold as long as the sufficient conditions hold asymptotically with high probability under $\mathcal{P}$.

The third part of our theorem provides further insights into the posterior distribution, assuming that the naive mean field approximation is “dominated" by a unique factorized distribution. In the language of Statistical Physics, these distributions are referred to be in a “pure phase".

We now demonstrate a few applications of Theorem 1 to concrete examples, which cover both deterministic and random design matrices. We defer the proofs of these corollaries to Appendix C.

Suppose that the $i^{th}$ row of the design matrix ${\mathbf{X}}$ equals $n^{-1/2}{\mathbf{x}}_{i}$, where $\{\mathbf{x}_{i}\}_{1\leq i\leq n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,\Gamma_{p})$. Assume that $p=o(n)$, and the following conditions hold:

The off diagonal part $\Gamma_{p,{\rm off}}$ of the covariance matrix $\Gamma_{p}$ satisfies ${\rm tr}(\Gamma_{p,{\rm off}}^{2})=o(p)$.

Sparse design matrices arise routinely in coding theory [Gal62, Mac99] and genomics [WLC+11, TJL+14]. Our next corollary discusses the accuracy of the nVB approximation in the context of a linear regression problem with a sparse bernoulli design. We assume that the entries of the design are independent, but not necessarily identical.

2 Scaling limit for the log-normalizing constant

Under the asymptotic validity of the naive mean-field approximation, we derive an asymptotic scaling limit for the log-normalizing constant. We will subsequently establish that this limiting description captures crucial information regarding the behavior of the optimizers at finite $n,p$. To this end, we will require some notation.

The cut norm of a graphon $W$ is given by

The cut norm is equivalent to the $L^{\infty}\mapsto L^{1}$ operator norm defined by

More precisely, we have $\lVert{W}\rVert_{\square}\leq\lVert{W}\rVert_{\infty\mapsto 1}\leq 4\lVert{W}\rVert_{\square}.$ It also follows from (9) that the cut norm is weaker than the $L^{1}$ norm, i.e. convergence in $L^{1}$ implies convergence in cut norm.

Given a symmetric $p\times p$ matrix $B$ with real entries, define a piecewise constant function on $^{2}$ by dividing $^{2}$ into $p^{2}$ smaller squares each of length $1/p$, and set

We will also need the following notion for embedding vectors into functions.

Suppose we are in the setting of Theorem 1. Assume further that

Then there exists $(\xi_{1},\cdots\xi_{p})\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,1)$ such that

We will derive a limiting formula for the log-normalizing constant (2) in terms of a variational problem on a space of probability distributions. This requires the following definition.

where $(X,X^{\prime})\stackrel{{\scriptstyle i.i.d.}}{{\sim}}U$ and $(Z,Z^{\prime})\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,1)$ are mutually independent.

$d_{\square}(W_{pA_{p}},W)\to 0$ for some $W\in\mathcal{W}$.

$w_{\bm{\beta}_{0}}\stackrel{{\scriptstyle L^{1}}}{{\rightarrow}}\phi(\cdot)$, and $w_{{\bf D}}\stackrel{{\scriptstyle L^{1}}}{{\rightarrow}}\psi(\cdot)$.

Lemma 4 and Theorem 2 are established in Section 4.2.

3 Uniqueness of the optimizer

Theorem 1 identifies general conditions for the asymptotic tightness of the naive mean-field lower bound to the log-normalizing constant. The Gibbs Variational Principle [WJ08] establishes that under these settings, the distribution $\mu$ can be approximated, to the leading order, by a product distribution. However, this does not specify whether the “best" approximation is unique. Indeed, the ferromagnetic Ising model on the complete graph, henceforth referred to as the Curie-Weiss model, provides a classical example where the mean-field lower bound is tight, but without an external magnetization, the model has two distinct optimizers at “low temperature".

In this section, we identify some conditions which guarantee the uniqueness of the minimizers. From a statistical perspective, these conditions allow us to conclude that the posterior distribution roughly behaves like a product distribution. To this end, our next lemma identifies a set of sufficient conditions.

Suppose all assumptions of Theorem 2 hold. Assume further that there exists ${\bf u}_{p}^{*}\in^{p}$ such that for all $\varepsilon>0$ there exists $\delta^{\prime}>0$ such that

The limiting variational problem (10) has a unique optimizer $F^{*}\in\mathscr{F}$.

$L_{p}^{(\mathbf{u}_{p}^{*})}\stackrel{{\scriptstyle P}}{{\to}}\mu^{*},$ where $\mu^{*}$ is the law of $(X,Z,F^{*}(X,Z))$ with $X\sim U$ and $Z\sim N(0,1)$ mutually independent.

$F^{*}$ satisfies the fixed point equation:

where $X\sim U$ and $Z\sim N(0,1)$ are mutually independent.

Combining Theorem 1 and Theorem 3, we get the following corollary, which deduces a Law of Large numbers under the posterior distribution.

Suppose (6), (11), and the three conditions (a), (b), (c) of Theorem 2 hold. Then the optimization in Theorem 2 has a unique solution $F^{*}\in\mathscr{F}$ that satisfies

To use this corollary in concrete examples, one needs to compute the functions $(W,g,\psi)$, and verify the conditions (5), (6) and (11). Of these, the separation condition (11) is somewhat implicit, while the other two conditions are relatively easy to verify directly. The following lemma provides two sufficient conditions to this end. In particular, it shows that (11) holds in the so called “high temperate regime”, or if $\pi$ has a density (with respect to Lebesgue measure) which is log concave.

Suppose there exists $\lambda>0$ and ${\bf u}_{p}^{*}\in^{p}$ such that for all ${\bf u}\in^{p}$ we have

for some $\lambda>0$, free of ${\bf u}$ and $p$. Then the condition (11) in Theorem 3 holds.

In particular (13) holds under either of the following conditions:

Let the prior $\pi$ have a density with respect to Lebesgue measure on $$,

We collect the proofs of Theorem 3 and Corollary 5 in Section 4.3.

4 Applications

To illustrate the utility of Theorem 2 and Corollary 4, we apply our results to specific examples in this section. The proofs are deferred to Appendix C.

Spiked Covariance Matrix: We consider linear regression with mean-zero gaussian features. We assume a spike covariance structure [JL09] on the features.

Consider the following two cases: either (i) $\lambda<\sigma^{2}$, or (ii) the conditions of Lemma 5 Part (b) (ii) hold. Then the conclusions of Corollary 4 hold.

The next example re-visits the sparse bernoulli design setting introduced in Corollary 3.

Suppose $\sigma^{2}>{\rm ess}\,\sup S(X)$, where $X\sim U$;

the conditions of Lemma 5 Part (b) (ii) hold.

Then the conclusions of Corollary 4 hold.

As our last example, we consider a sequence of design matrices arising in the study of two-way ANOVA designs.

Suppose that we have a two factor ANOVA model of the form

Assume that the true regression coefficient $\bm{\beta}_{0}\in^{p}$ satisfies $w_{\bm{\beta}_{0}}\stackrel{{\scriptstyle L^{1}}}{{\to}}\phi$, for some function $\phi\in L^{1}$.

and $g(x)=\int_{}W(x,y)\phi(y)dy+\frac{1}{2}\phi(x)$.

the conditions of Lemma 5 Part (b) (ii) hold.

Then the conclusions of Corollary 4 hold.

Discussions

We discuss some limitations of our current results, and collect some questions for future enquiry.

The bounded support assumption on the prior—A vital technical assumption in our analysis concerns the bounded support assumption on the prior. We note that for a general prior with unbounded support, the posterior $\mu_{\mathbf{y},\mathbf{X}}$ might not even be a proper probability distribution. One intuitively expects that under appropriate “tail-decay" conditions on the prior, the results in this paper should generalize. Going beyond the bounded support assumption requires extending the theory of non-linear large deviations to probability measures on unbounded spaces, and is thus beyond the scope of this paper.

Extensions to Gibbs posteriors and fractional posteriors—A careful study of our proof reveals that our main results do not depend strongly on the correct model specification. As a result, we expect similar techniques to be broadly useful in the study of Gibbs and fractional posteriors [ARC16, YPB20].

Extensions to other GLMs—Another natural question of interest concerns the applicability of these ideas to more general models, e.g. logistic regression. We consider this to be an extremely interesting question, and plan to explore this in the future.

Extensions to models with latent characteristics—Variational methods are ubiquitous in applications with latent characteristics e.g. topic modeling, community detection etc. In contrast, the relevant variables are all observed in the linear model framework. It will be interesting to explore the applicability of our techniques to models with latent features.

Proofs

We prove the main results in this section. Theorem 1 is established in Section 4.1, Theorem 2 is established in Section 4.2, while Theorem 3 is proved in Section 4.3.

Our first lemma collects some basic facts about the exponential family $\psi_{\gamma}$. The proof is deferred to the Appendix.

In the setting of Lemma 1, we have the following conclusions.

We now turn to the proof of Theorem 1. To this end, set

We will prove that certain statistics under the posterior distribution $\mu$ can be “well-approximated" by the vector of conditional means. To this end, we establish the following results.

We defer the proofs of these lemmas to the end of this section, and establish Theorem 1, given these lemmas.

Proof of Part (i). Since $\text{tr}(A^{2})=o(p)$, [Yan20, Theorem 4] implies

Note that it suffices to establish that for all $\varepsilon>0$, $\mu({\bf b}\in C_{p}(\varepsilon))\to 0$ as $p\to\infty$. To this end, define the event

Note that (20), in combination with Chebychev inequality, implies that $\mu(E_{p}(\varepsilon/2)^{c})\to 0$ as $p\to\infty$. Therefore,

where the first inequality follows using the definition of $E_{p}(\varepsilon/2)$, while the last inequality follows from the definition of $C_{p}(\varepsilon)$. Using the display above in combination with (7), it suffices to establish that

This argument would be relatively straight-forward if the $\theta_{i}$ were fixed constants independent of $\bm{\beta}$, as

We caution the reader that this is not the case, and the $\theta_{i}$ are themselves functions of $\bm{\beta}$, as specified in (16).

For $\bm{\beta}\in\mathcal{P}(\mathbf{p})$, setting $\theta_{i}^{\mathbf{p}}=\frac{z_{i}-p_{i}}{\sigma^{2}}$, we have,

where the last inequality uses the definition of $\mathcal{P}(\mathbf{p})$. This gives,

where the first inequality follows from Cauchy-Schwarz, and the second inequality follows from the smoothness of $c(\cdot,\cdot)$. Thus we have,

This concludes the proof, as $\delta>0$ is arbitrary.

where $H$ is defined as in (14).To complete the proof, it thus suffices to show that

Fixing $K<\infty$ and setting $d_{K}(i):=d_{i}1\{|d_{i}|\leq K\}$, using Lemma 6 Part (vi) we have

which goes to as $p\to\infty$ followed by $K\to\infty$, using the uniform integrability assumption on $\frac{1}{p}\sum_{i=1}^{p}\delta_{d_{i}}$. Also, with $\delta>0$ and setting

Fixing $K>0$, we now claim that the RHS of (4.1) converges to as $\delta\to 0$. Given this claim, it follows from (23) and (4.1) that

A similar argument with $\hat{u}_{i}$ replaced by $b_{i}$ gives

where $b_{\delta}(i):=b_{i}1\{|b_{i}|\leq 1-\delta\}$.

It suffices to show that for every $\delta,K$ fixed we have

which converges to zero in probability under the posterior distribution $\mu_{\mathbf{y},\mathbf{X}}(\cdot)$. Here the last estimate uses part (ii) of this Theorem.

To complete the argument, it thus remains to verify the claim involving (4.1), for which it suffices to verify continuity of the function $(x,y)\mapsto H(x,y)$ at $x=\pm 1$, uniformly for $y\in[-K,K]$. By symmetry, it suffices to show that

Suppose this is not true. Then there exists sequences $\{x_{k}\}_{k\geq 1},\{y_{k}\}_{k\geq 1}$ with $\lim_{k\to\infty}x_{k}=1$, and $|y_{k}|\leq K$, such that $|H(x_{k},y_{k})-1|>\varepsilon$ for all $k$, for some $\varepsilon>0$. Without loss of generality, by passing to a subsequence, we can assume that $y_{k}$ converges to $y\in[-K,K]$. Using Lemma 6 Part (iii) we have that $\pi_{(h(x_{k},y_{k}),y_{k})}$ converges weakly to $\delta_{1}$, the point mass at $1$. Consequently, using DCT we have

a contradiction. This verifies the claim, and hence completes the proof of part (iii).

Next we turn to the proof of Lemmas 7 and 8.

Since (18) and (19) follow immediately from [Yan20, (4.40)] and [Yan20, (4.41)] respectively, it suffices to prove (20). But this follows upon observing that

For any $i,j\in[p]$ with $i\neq j$, setting

Since $t_{j}^{(i)}$ does not depend on $\beta_{i}$, we have

For estimating the second term in the RHS of (25), setting

Using the last two displays and summing over $i,j$ give

2 Proof of Theorem 2

Without loss of generality assume that the samples $(Z_{1},\cdots,Z_{p})$ are arranged in increasing order of variance, i.e. increasing order of $\{D_{p}(i,i)\}_{i=1}^{p}$. Let $\{\delta_{p}:p\geq 1\}$ be a positive sequence converging to , such that

The existence of such a sequence follows from the assumption (5). Let

where we use $\sum_{i=1}^{p}D_{p}(i,i)=O(p)$. For $q\leq i\leq r$, setting $Y_{i}:=\frac{Z_{i}}{\sqrt{D_{p}(i,i)}}$, let $\Gamma_{p}$ denote the $(r-q+1\times r-q+1)$ dimensional covariance matrix of the vector ${\bf Y}:=(Y_{q},\cdots,Y_{r})^{\top}$. Define $\mathcal{E}=(\xi_{q},\cdots,\xi_{r})^{\top}:=\Gamma_{p}^{-1/2}{\bf Y}$, and note that $(\xi_{q},\cdots,\xi_{r})\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,1)$. Setting $v_{i}:=u_{i}\sqrt{D_{p}(i,i)}$ and ${\bf v}:=(v_{q},\cdots,v_{r})$, we have

where $B:=\Gamma_{p}^{1/2}-{\bf I}_{r-q+1}$. We now claim that

Finally with $\{\xi_{i}\}_{[p]\backslash[q,r]}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,1)$ independent of $\{\xi_{i}\}_{i=q}^{r}$, using arguments similar to the derivation of (26) we get

Combining (26), (29) and (30) the desired conclusion follows.

It thus remains to verify (28). To this effect, note that $\Gamma_{p}(i,i)=1$, and $\Gamma_{p}(i,j)=\frac{A_{p}(i,j)}{\sqrt{D_{p}(i,i)D_{p}(j,j)}}$. Setting $C_{\delta}(i,j):=\Gamma_{p}(i,j)1\{i\neq j\}$, and using Spectral Theorem, we have,

where C:=\sup_{x\geq-1}\Big{|}\frac{\sqrt{1+x}-1}{x}\Big{|}<\infty, and the last bound uses (31). The last term above is $o(p)$ by the choice of $\delta_{p}$, and so we have verified (28). This completes the proof of the Lemma. ∎

To establish Theorem 2, we need the following definitions.

The following stability estimates will be crucial in our proof of Theorem 2. The proof is deferred to the Appendix.

We have, for $W,W^{\prime}\in\mathcal{W}$, $\phi,\psi,\phi^{\prime},\psi^{\prime}\in L^{1}()$,

$W_{k},W\in\mathcal{W}$ is such that $d_{\square}(W_{k},W)\to 0$.

$\phi_{k},\phi\in\mathcal{L}$ is such that $\int_{}|\phi_{k}(x)-\phi(x)|dx\to 0$.

$\psi_{k},\psi\in L_{1}$ is such that $\int_{}|\psi_{k}(x)-\psi(x)|dx\to 0$.

We have, for any $W\in\mathcal{W}$, $g,\psi\in L^{1}()$,

Then there exists a sequence $\delta_{p}\to 0$ as $p\to\infty$ such that setting

The proof of Lemma 11 is straightforward, and is thus omitted.

Suppose we are in the setting of Theorem 2. Fix $F\in\mathscr{F}$ and $\delta,K>0$ and $p\geq 1$. Let $V_{i}\sim U([(i-1)/p,i/p])$ be independent of $\xi_{i}$ arising in Lemma 4. Set $d_{i}^{K}=d_{i}\mathbf{1}(|d_{i}|\leq K)$ and define

We prove Theorem 2 assuming Lemma 10, 12 and 13. The corresponding proofs are deferred to the end of the section.

where $\sup_{\mathbf{u}\in^{p}}|\mathcal{E}_{p}^{(2)}({\bf u})|\stackrel{{\scriptstyle P|\mathbf{X}}}{{\to}}0$ as $p\to\infty$. Here, $(b)$ uses the definition (32), and $(c)$ uses Lemma 9 part (ii). To invoke Lemma 9, we use the fact that

To complete the upper bound, invoking Lemma 10, it suffices to establish that for all $\varepsilon>0$

We first turn to the proof of (37) by contradiction. Suppose there exists $\varepsilon>0$, such that setting

Lower Bound: It suffices to show that for all $\varepsilon>0$,

To this end, let $F\in\mathscr{F}$ satisfy

where $\psi_{K}(x):=\psi(x)\mathbf{1}(\psi(x)\leq K)$, and we have used Lemma 9 Part (i). Further, setting $d_{i}^{K}=\min\{d_{i},K\}$, a similar argument gives

For any $p\geq 1$, let $\xi_{1},\cdots,\xi_{p}$ denote the iid $\mathcal{N}(0,1)$ random variables arising in the optimization problem $M_{p}$. For each $i\in[p]$, generate $V_{i}\sim U([(i-1)/p,i/p])$ independent of each other, and of the $\xi_{i}$ variables. Fixing $\delta>0$, we set

Using Lemma 13, for any $\varepsilon>0$, there exists $\delta,K>0$ (depending on $\varepsilon$), such that with probability at least $1-\varepsilon$,

As $\varepsilon>0$ is arbitrary, this completes the proof of the lower bound. ∎

It remains to prove Lemma 10, 12 and 13. We establish each result in turn.

where the final step follows from Jensen’s inequality and the observation that $\frac{\partial^{2}}{\partial u^{2}}G(u,d)\geq 0$ from Lemma 2. This implies

Finally, noting that $\Phi^{-1}$ is strictly increasing, we have,

where for each $x\in$, $F^{*}(x,\cdot)$ is the monotone rearrangment of $F(x,\cdot)$. This last step follows from Hardy-Littlewood inequality. Observe that as the Uniform distribution is invariant under measure preserving transformations, the other terms in the functional above are unchanged by this rearrangement. Thus

Since $F^{*}\in\mathscr{F}$, the proof is complete.

The desired conclusion follows once we establish that

But this follows on using Fatou’s lemma and the lower semicontinuity of G\Big{(}\cdot,\frac{\psi(X)}{\sigma^{2}}\Big{)} (Lemma 6 part (iv)). ∎

The definition of $\mathbf{u}$ in (35) implies

Combining (42), (43), (44) and (45), for any $\varepsilon>0$, we have,

We now claim that for any $\varepsilon>0$,

We first turn to the quadratic form. We have,

by the weak law of large numbers, and the observation $|\mathcal{M}_{p}(i,j)|\leq K/p^{2}$. Here, the $o_{P}(1)$ term converges to zero in probability under the joint distribution of $\{(V_{i},\xi_{i}):1\leq i\leq p\}$. Further, we have,

where $X\sim U()$ and $Z\sim\mathcal{N}(0,1)$ are independent. Finally,

where (a) uses $\|W_{\mathcal{M}_{p}}-W^{(K)}\|_{1}\to 0$ as $p\to\infty$ and (b) uses $\|W-W^{(K)}\|_{1}\to 0$ as $K\to\infty$. Combining (48), (49), (50), (51), the first conclusion of (47) follows upon sending $p\to\infty$, $\delta\to 0$ and $K\to\infty$.

The other conclusions of (47) follow using analogous arguments, and are thus omitted. ∎

3 Proof of Theorem 3 and related results

The following continuity statement will be critical for the proof of Theorem 3. The proof is deferred to the Appendix.

For any $W\in\mathcal{W}$, if $\nu_{p}\stackrel{{\scriptstyle w}}{{\to}}\nu$, then

We break the proof into several steps. Note that using Lemma 10, both variational problems attain their suprema.

We establish uniqueness of the maximizer assuming (11). If possible, let $F_{1}$ and $F_{2}$ be two distinct optimizers of $\mathcal{G}_{W,g,\psi}$ in $\mathscr{F}$. Therefore, using Theorem 2,

Let $\xi_{1},\cdots,\xi_{p}$ be iid $\mathcal{N}(0,1)$ random variables arising in the functional $\widetilde{M}_{p}(\cdot)$. Further, let $V_{i}\sim U([(i-1)/p,i/p])$ be independent of the $\xi_{i}$ variables. Fixing $\delta>0$, following the recipe in (35), we define two sequences $\mathbf{u}^{(1)}_{p,\delta,K}$, $\mathbf{u}^{(2)}_{p,\delta,K}$ in $^{p}$ corresponding to $F_{1},F_{2}$ respectively. Using Lemma 13, for $i=1,2$, for $\delta,K>0$ we have

where $\limsup_{K\to\infty}\limsup_{\delta\to 0}\mathcal{E}(\delta,K)=0$.

Moreover, using Theorem 2 and the assumption that $F_{1}$ and $F_{2}$ are optimizers of the limiting variational problem, we have, for $i=1,2$,

On the event $\{\|\mathbf{u}_{p,\delta,K}^{(i)}-\mathbf{u}_{p}^{*}\|_{2}^{2}>p\varepsilon,i=1,2\}\cap\{\frac{1}{p}M_{p}(\mathbf{u}^{*}_{p})<\mathcal{G}_{W,g,\psi}(F_{i})+\frac{\delta^{\prime}}{4}\}$, for $i=1,2$ we have,

Thus, with $X\sim U()$, $Z\sim\mathcal{N}(0,1)$ independent, we have,

Finally, note that $\varepsilon>0$ is arbitrary, and thus $F_{1}=F_{2}$ a.s.. This establishes the desired uniqueness.

Recalling the 2-Wasserstein distance $d(\cdot,\cdot)$ from Lemma 11, we have,

where the last equality follows from the law of large numbers and the boundedness of $\eta$. Here, the $o_{P}(1)$ term converges to zero in probability under the joint distribution of $\{(V_{i},\xi_{i}):1\leq i\leq p\}$. Finally,

For any $p\geq 1$, the function $\widetilde{M}_{p}(\cdot)$ is upper semicontinuous on $^{p}$, and thus attains its maximum. Let $\{\hat{\mathbf{u}}_{p}:p\geq 1\}$ denote any sequence of global maximizers of $\widetilde{M}_{p}(\cdot)$. Lemma 4 and the separation condition (11) together imply that $\frac{1}{p}\|\mathbf{u}_{p}^{*}-\hat{\mathbf{u}}_{p}\|_{2}^{2}\stackrel{{\scriptstyle P|\mathbf{X}}}{{\to}}0$. In turn, this implies

Next, differentiating $\widetilde{M}_{p}(\cdot)$ at $\mathbf{u}\in(-1,1)^{p}$, we obtain

where we use $\ddot{c}(\theta,d)\leq 1$, and observe that

Then, upon observing that $\sup_{\gamma_{1},\gamma_{2}}|\frac{\partial^{2}}{\partial\gamma_{1}\partial\gamma_{2}}c(\gamma_{1},\gamma_{2})|\lesssim 1$ and $\sup_{\gamma_{1},\gamma_{2}}|\frac{\partial^{2}}{\partial\gamma_{1}^{2}}c(\gamma_{1},\gamma_{2})|\lesssim 1$, we have, on the event $\mathscr{E}(\varepsilon)$,

where the last inequality uses Lemma 14 part (i).

Combining (54), (55), (56), (57), we have,

Sending $p\to\infty$, and then $\varepsilon\to 0$, we obtain

The desired conclusion follows upon recalling that $(X,Z,F^{*}(X,Z))\sim\nu^{*}$.

We next turn to the proof of Lemma 5. To this end, we require the following auxiliary lemma. The proof is deferred to the Appendix.

where $H_{p}(\cdot)$ is the Hessian of $r_{p}(\cdot)$. Then there exists a unique global maximizer ${\bf x}_{0}\in^{p}$, and further for any ${\bf x}\in^{p}$ we have

Proof of (i) The equation (13) implies that the event

Proof of (ii)(a) Differentiating $M_{p}(\cdot)$, twice, the Hessian is given by

Note that $\ddot{c}(h(u,d),d)$ is the variance of a random variable supported on $$, and thus$\ddot{c}(h(u,d),d)\leq 1$. Under the assumption of (i), there exists$\delta>0$such that for all$p$large enough and all${\bf x}\in(-1,1)^{p}$,

It then follows from Lemma 15 that there exists ${\bf u}^{*}\in^{p}$ which is a unique optimizer of $M_{p}(\cdot)$, and further, for any ${\bf u}\in^{p}$ we have

Consider now a scale family of parametric distributions $\{\mathcal{P}_{\theta}:\theta\geq 0\}$ with

Since $V(\cdot)$ is even, it follows that the first moment under $\mathcal{P}_{\theta}$ is for all $\theta$. Also, since $V(\cdot)$ is an increasing, this family has monotone likelihood ratio in $T(x)=|x|$, and thus the second moment under the law $\mathcal{P}_{\theta}$ is decreasing in $\theta$ for $\theta>0$, and so

Proceeding to bound the RHS of the above display, for $d>0$ we have

This is exactly the variance of a truncated standard Gaussian distribution, truncated to the interval $[-\sqrt{d},\sqrt{d}]$. Indeed, the truncated variance is $1-\frac{2\sqrt{d}\phi(\sqrt{d})}{2\Phi(\sqrt{d})-1}<1$ [JK71], where $\phi(\cdot)$ and $\Phi(\cdot)$ represent the pdf and cdf of the standard Gaussian distribution. ∎

References

Appendix A Some technical lemmas

We collect some basic results about exponential families in the first subsection. We prove Lemma 15 in the next subsection.

We prove Lemmas 1, 2 and 6 in this section.

This follows by direct computation (see e.g. [LC06]).

Without loss of generality, we consider the case $\gamma_{1,k}\to\infty$. The case $\gamma_{1,k}\to-\infty$ follows using the same argument, with obvious modifications. Observe that as $\gamma_{1,k}\to\infty$ and $\limsup_{k}|\gamma_{2,k}|<\infty$, for $k$ sufficiently large, the function $\gamma_{1,k}z-\frac{\gamma_{2,k}}{2}z^{2}$ is increasing on $$. This implies, for any$\varepsilon>0$,

On taking ratios of the above two displays we get

Also note that the ratio of probabilities under $\pi$ in the RHS is positive, as $\pm 1$ is in the support of $\pi$. The desired conclusion now follows upon letting $k\to\infty$, upon noting that $\gamma_{1,k}\to\infty$, and $\gamma_{2,k}$ stays bounded.

The proof follows directly from the lower semicontinuity of KL divergence under weak convergence [Pos75, Theorem 1].

For $x\in\{-1,1\}$, $H(x,\cdot)$ is independent of $y$, and thus $\frac{\partial H(x,y)}{\partial y}=0$. We assume henceforth that $x\in(-1,1)$. Differentiating, we obtain,

where the last equality uses $\dot{c}(h(x,y),y))=x$. Therefore,

We complete the proof by contradiction. Suppose, if possible, that $M:=\lim_{t\to 1}h(t,\gamma_{2})<\infty$. In this case, as $-1$ is in the support of $\pi$, $\pi()>0$. This implies that

To complete the argument, it suffices to show that $\liminf_{t\to 1}\pi_{(h(t,\gamma_{2}),\gamma_{2})}()>0$. But this follows on noting that the function $t\to c(h(t,\gamma_{2}),\gamma_{2})$ is non-decreasing on $$, and thus

The lemma follows by direct computation. First, note that

where the last equality follows upon noting that $\dot{c}(h(u,d),d)=u$. For the second derivative, note that

where the last equality follows from differentiating the equation $\dot{c}(h(u,d),d)=u$. ∎

A.2 Proof of Lemma 15

Since $r_{p}(\cdot)$ is upper semi-continuous there exists a global maximizer in $^{p}$, say $\widetilde{\bf x}_{0}$. Fixing any ${\bf x}$ in $(-1,1)^{p}$, consider the function $f(t):=r_{p}((1-t)\widetilde{\bf x}_{0}+t{\bf x})$. Then $f$ is twice differentiable on $(0,1)$, and

Consequently, for any $\varepsilon\in(0,1)$ Taylor’s expansion implies

Taking limits as $\varepsilon\to 0$ we get

where the last inequality uses the fact that $f^{\prime}(0+)\leq 0$, as $t=0$ is a global maximizer of $f$. The last display is equivalent to

This inequality then extends to $x\in^{p}$ by upper semi-continuity of $r_{p}$. ∎

Appendix B Stability Estimates for functionals

Next, we turn to the proof of Part (ii). Using the definition of cut norm (as in Definition 4) and Part (i), we have,

Using Part (i) again, it suffices to show that

To this end, note that $|W.\phi_{k}+\phi_{k}(x)\psi(x)-W.\phi-\phi(x)\psi(x)|$ converges to in measure, and

which is an integrable function. This completes the argument using DCT.

B.2 Proof of Lemma 14

For any $\varepsilon>0$, let $S^{+}(\varepsilon)=\{x:m(\mu,W,x)-m(\mu,W^{\prime},x)>\varepsilon\}$. For $X\sim U()$, we have,

Setting $S^{-}(\varepsilon)=\{x:m(\mu,W,x)-m(\mu,W^{\prime},x)<-\varepsilon\}$, the same argument now yields that

using the uniform continuity of $W^{\prime}$. The desired conclusion follows upon combining the two displays above.

Appendix C Proofs of Examples

We establish Corollaries 1-3 and 5-7 in this section. Throughout this section $o_{P}(1)$ terms converge to zero in probability under the marginal distribution of the design matrix $\mathbf{X}$.

Since $\{{\bf x}_{i}\}_{1\leq i\leq n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,\Gamma_{p})$, and

where the last equality uses $p=o(n)$. Noting that

Consequently, $A_{p}$ satisfies (5) with high probability, as ${\rm tr}(\Gamma_{p,{\rm off}}^{2})=o(p)$. Further, since $\lVert\Gamma_{p}\rVert_{2}=O(1)$, it follows from (65) gives that $\lVert A_{p}\rVert_{2}=O_{P}(1)$. Finally, noting that

we have $\frac{1}{p}\sum_{i=1}^{p}\delta_{D_{p}(i,i)}$ is uniformly integrable with high probability. This completes the proof of the corollary.

It suffices to verify the same conditions on the matrices $(A_{p},D_{p})$ as in Corollary 2.

To this end, note that for any $i\neq j$ we have

which is $o(p)$ as $p=o(n)$. This verifies (5).

and so (6) holds with high probability. Finally we have

and so $\frac{1}{p}\sum_{i=1}^{p}D_{p}(i,i)$ is uniformly integrable. ∎

C.2 Limiting variational formula

The desired conclusion follows from Theorem 2, once we can verify

We begin by verifying condition (11). To this effect, set

and use the fact that $\max_{i\in[p]}\sum_{j\neq i}\Gamma_{p}(i,j)\leq\lambda<\sigma^{2}$ along with part (a) of Lemma 5 to get the existence of ${\mathbf{u}}_{p}^{*}$ such that

This in turn shows that for any $\varepsilon>0$ we have

Thus we have verified (11). The desired conclusion then follows from Corollary 4.

Part (b)(ii) follows from part (b)(ii) of Lemma 5 and Corollary 4.

The desired conclusion follows from Theorem 2, once we can verify

Proceeding to verify the above display, note that

Thus, setting $\widetilde{D}_{i}:=\frac{1}{n}\sum_{k=1}^{n}G(k/n,i/p)$, Using Chernoff bounds it follows that

where the last equality uses the almost sure continuity of $G(.,.)$.

To this end, note that for any $i\neq j$ we have

From this, using the condition $p=o(\sqrt{n})$ gives

where the last equality again uses the almost sure continuity of $G(\dot{,}\dot{)}$.

(i) Once again, the desired conclusion of part (b) follows from Corollary 4, once we verify condition (11).

To this effect, fixing $\delta>0$ and setting $R_{i}:=\sum_{j\neq i}^{p}A_{p}(i,j)$, define a $p\times p$ matrix $B_{p,\delta}$ by setting

We now claim that there exists $\delta>0$ such that

But this is immediate from Lemma 5 part (a), on noting that

It thus remains to verify (70). To this effect, recall from part (a) that $W_{pA_{p}}$ converges to $W$ in the cut metric, which in turn implies $\frac{1}{p}\sum_{i=1}^{p}\delta_{R_{i}}$ converges weakly in probability to the law of $S(X)$, where $X\sim U$ [BCCG15, Theorem 2.16]. This gives

It thus suffices to show that there exists $\delta>0$ such that the RHS above is . But this follows on noting that ${\rm esssup}(S(X))<\sigma^{2}$. This completes the proof of part (i).

Part (b)(ii) follows from part (b)(ii) of Lemma 5 and Corollary 4, as before.

It then follows that $d_{\square}(W_{pA_{p}},W)\to W$, and $d_{L_{1}}(w_{\bf D},\psi)\to 0$. The desired conclusion then follows from Theorem 2 as before.

Since $\limsup\limits_{p\to\infty}\max_{i\in[p]}\sum_{j\neq i}|A_{p}(i,j)|=\frac{1}{2}$, the result is immediate from Corollary 4 and Lemma 5.

Appendix D Relevant concentration inequalities

To prove Lemma 16, we first obsere that if $X$ is sub-exponential, then $X^{2}$ is sub-Weibull [KC18].

Let $Y_{1},\cdots,Y_{n}$ be independent sub-exponential random variables with common sub-exponential parameter $c>0$. For any $\delta>0$, there exists $C:=C(\delta)>0$ such that

We first claim that $Y_{i}^{2}$ are sub-Weibull, i.e., there exists a constant $C^{\prime}>0$ (depending only on $c$) such that

Using (72), along with the deviation bound above, we have,

for some constant $C^{\prime}>0$. The last inequality uses that $Y_{i}$ is sub-exponential. (72) can be verified by direct computation.

Using (71) and [KC18, Proposition A.3], we have, $\max_{i\leq n}\|Y_{i}^{2}\|_{\psi_{1/2}}$ is uniformly bounded in $n$. Consequently, using [KC18, Theorem 3.1], for any $t>0$, we have,

where $C_{1}$ is a constant free of $n$, and $L_{n}=C_{2}\|b\|_{\infty}/\|b\|_{2}$ for some constant $C_{2}>0$ independent of $n$. Here, $b=(\|Y_{1}^{2}\|_{\psi_{1/2}},\cdots,\|Y_{n}^{2}\|_{\psi_{1/2}})$. We set $C_{1}\|b\|_{2}(\sqrt{t}+L_{n}t^{2})=n\delta$ for some $\varepsilon>0$. Direct calculation yields that $t(\delta)=\sqrt{\frac{n\delta}{C_{1}C_{2}\|b\|_{\infty}}}(1+o(1))$, so that

where $C_{3}>0$ depends on $\delta$. This concludes the proof. ∎

By Chernoff bounds, for any fixed $S\subseteq[p]$, $Y_{i}=\sum_{k\in S}B_{ik}$ is a sub-exponential random variable. Using Lemma 17,

The required conclusion follows from (74) and the deviation bound (73). It thus remains to prove (74). To this effect, note that

which on using triangle inequality gives (74).