1 Introduction↩︎

Laplace-type integrals of the form \[\label{lap} I(\lambda):=\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d} g(x)e^{-\lambda f(x)}\mathrm{d}x\tag{1}\] are a fundamental object in asymptotic analysis and a ubiquitous tool for deriving tractable approximations to otherwise intractable quantities, such as normalizing constants and expectations. When \(\lambda\) is large, these integrals are dominated by neighborhoods of minimizers of \(f\), and Laplace asymptotics yields explicit leading-order formulas and systematic higher-order corrections. Such approximations underpin both rigorous analysis (e.g. sharp tail probabilities and free-energy expansions) and practical computation in settings where direct numerical integration is infeasible.

In the classical regime where the dimension \(d\) is fixed and \(\lambda\to\infty\), Laplace’s method and its higher-order refinements are well developed [1], [2]. At the opposite extreme, there is also an infinite-dimensional theory (e.g. on Wiener space) in which \(\lambda^{-1}\) plays the role of a small-noise parameter [3]. What remains comparatively less understood is the intermediate regime in which the dimension grows with the large parameter. This “growing-\(d\)” regime is now very common in modern statistics, where Laplace-type integrals (or ratios thereof) appear as marginal likelihoods and posterior expectations in models whose parameter dimension increases with sample size. See Section [ssec:statistics95motivation] for more details. Likewise, the absence of rigorously justified “growing-\(d\)” expansions created a significant void in areas such as statistical physics, Euclidean quantum field theories (QFT), and chemistry, where formal expansions have been in use for a long time without adequate justification, see Section 1.1.

Recent work by the second author establishes a high-dimensional Laplace expansion (LE) of \(I(\lambda)\), with explicit remainder control, under the scaling \(d^2/\lambda\to 0\) [4]. The proof centers around an approximation of \(f\) in the exponent by its second-order Taylor expansion. Thus the expansion of \(I(\lambda)\) is closely tied to the problem of approximating concentrating densities, of the form \(\pi\propto e^{-\lambda f}\), by Gaussian distributions. In fact, the condition \(d^2\ll\lambda\) arises in numerous contexts as a threshold for Gaussian approximation accuracy. In the context of Gaussian approximation to posteriors in Bayesian inference, \(d^2\ll\lambda\) appears in the works [5]–[8]. The condition \(d^2\ll\lambda\) also arises in proofs of high-dimensional Central Limit Theorems (CLTs); see [9]–[11]. In the CLT, \(\lambda\) (usually denoted \(n\)) is the given number of i.i.d. \(d\)-dimensional random vectors, whose average is of interest. In fact, the aforementioned works [4], [5], [9] establish lower bounds as well, proving \(d^2/\lambda\ll1\) is necessary for the Gaussian approximations to be accurate. Although some works have in fact proved the CLT under much weaker conditions on \(d\) relative to \(\lambda\), the Gaussian approximation holds in a much weaker sense. See Section 3.1 and Remark 2 of [12] for an overview of this line of work.

Thus to summarize, the regime \(d^2/\lambda\to0\) arises as a critical threshold in numerous results on approximating \(d\)-dimensional densities or integrals involving \(e^{-\lambda f}\). Another critical threshold is the regime \(d/\lambda\to0\), which is generically required for distributions \(\propto e^{-\lambda f}\) to concentrate near the minimizer of \(f\); see e.g. [13] and [14]. Clearly, if there is no concentration around the minimizer of \(f\), one should not expect an expansion to hold whose terms are determined solely by the derivatives of \(g\) and \(f\) at this single point.

This leaves an intriguing intermediate region where \(d^2/\lambda\) does not vanish or even diverges to infinity, while \(d/\lambda\) still converges to zero. Due to concentration about the minimizer of \(f\), one expects that \(I(\lambda)\) in 1 can still be characterized by derivatives of \(f\) and \(g\) at the minimizer. But it is unclear how to obtain a closed form approximation to the integral without leaning on the tractable Gaussian integral obtained by replacing \(f\) with its second-order Taylor expansion. To our knowledge, there has been no rigorously justified explicit Laplace-type approximation of \(I(\lambda)\) in this intermediate regime.

In this paper we have achieved an important milestone by completely characterizing this largely unexplored intermediate regime under natural local regularity and global growth conditions. We derive an explicit asymptotic series approximating 1 such that for each expansion order \(L\geq1\), the remainder is negligible as long as \(d^{L+1}/\lambda^L\to0\). The series approximates the logarithm of the integral \(I(\lambda)\), and as we explain below, this is precisely what allows us to push \(d\) above the \(\sqrt\lambda\) barrier. Let \(x_\star\) be the global minimizer of \(f\) and assume without loss of generality that \(\nabla^2f(x_\star)=I_d\). Assuming bounded operator norms of derivatives of \(f\) and \(g\) near \(x_\star\) through orders \(2L+2\) and \(2L\), respectively, and mild global growth conditions on \(f\) and \(g\), we show that I()=_k=1^L-1 b_k(f,g)^-k+(d^L+1/^L) for some coefficients \(b_k(f,g)\) that satisfy \(b_k(f,g)=\mathcal{O}(d^{k+1})\). If \(d\lesssim\log\lambda\), then additional factors of \(\log\lambda\) are present in the remainder. The sum with respect to \(k\) is omitted if \(L=1\), and we assumed without loss of generality that \(g=1\) at \(x_\star\). The coefficients \(b_k\) are already well-known, and can be derived from formal cumulant expansions [15].

In fact, more broadly, our contribution is a powerful technique to tackle Laplace integrals. We use this technique not only to prove [main32res32intro], but also to solve a related problem of approximating a Laplace-type probability density, of the form \(\pi(x)\propto e^{-\lambda f(x)}\). Namely, we construct a transformation \(x_L:\mathbb{R}^d\to\mathbb{R}^d\) such that \(\hat{\pi}_L:=(x_L)_{\#}\mathcal{N}(0,\lambda^{-1}I_d)\) approximates \(\pi\): \[\label{TV-intro} \mathrm{TV}(\pi,\hat{\pi}_L )\lesssim d^{L+1}/\lambda^L.\tag{2}\] Here, \(\mathcal{N}(0,\lambda^{-1}I_d)\) is the Gaussian distribution with zero mean and covariance matrix \(\lambda^{-1}I_d\), while TV stands for the total variation distance. The subscript \(\#\) in the definition of \(\hat{\pi}_L\) denotes the push-forward. This means \(\hat{\pi}_L\) is the law of the random variable \(x_L(Z)\) when \(Z\) is distributed as \(\mathcal{N}(0,\lambda^{-1}I_d)\). The reason \(\hat{\pi}_L\) is so useful is that it is easy to sample from, precisely due to this push-forward construction. Cheaply generating approximate samples from an untractable density \(\pi\propto e^{-\lambda f}\) is an important problem in Bayesian statistics. See Section [ssec:statistics95motivation] for more details.

1.0.0.1 Significance.

A central significance of our result is that it advances the modern asymptotic analysis program of developing expansions of ubiquitous Laplace type integrals that remain accurate in increasingly high-dimensional regimes. The best known results to date work under the assumption \(d^2/\lambda\to0\). In contrast, we extend the range of dimensions for which one can make precise asymptotic statements with explicit formulas until the very limit, because beyond \(d/\lambda\to0\) there is no concentration any longer. In this sense, our work essentially completes the classical Laplace program for concentrating finite-dimensional integrals under natural smoothness and growth assumptions, by proving remainder bounds on the high-dimensional LE up to the concentration threshold. We show that a constructive analytic approximation remains valid whenever \(d^{L+1}/\lambda^L \to 0\) for any \(L\ge 1\), even in the genuinely intermediate region where \(d^2/\lambda\not\to 0\) (indeed, where \(d^2/\lambda\) may diverge), by identifying an explicit exponential correction at the level of the log-integral.

The implications of our results across many areas of science and engineering are numerous. Two particularly noteworthy applications deserve special mention.

The first application is in physics, including statistical physics and QFT, where Laplace-type integrals encode quantities of central importance. In the high-dimensional, many-degrees-of-freedom regimes relevant to these fields, such quantities are often evaluated via formal LEs, typically without appropriate remainder bounds; see Section 1.1 for details. This theoretical gap, which dates back at least to the Darwin–Fowler steepest-descent approach in 1922 [16], has long limited the rigor of many computations. Our results fill the century-old gap by placing these Laplace calculations on firm mathematical footing for a broad class of finite-dimensional large-system models.

Another area where our results are of significant importance is statistics. Here, \(\lambda\) plays the role of sample size, and high-dimensional Laplace-type integrals arise ubiquitously as normalizing constants, marginal likelihoods, and posterior expectations. As mentioned above, a central task is to approximate the posterior density \(\pi \propto e^{-\lambda f}\): one wants to sample from it, compute expectations of observables against it, and evaluate its normalizing constant for model comparison. Our results address all three. We construct an explicit approximation \(\hat{\pi}_L\) to \(\pi\) from which one can easily sample, we provide closed-form approximations to posterior expectations of smooth observables that avoid Monte Carlo error entirely, and we give rigorous asymptotic expansions of the normalizing constant that generalize the Bayesian Information Criterion (BIC) to higher order. Laplace-based surrogates such as the BIC are popular precisely because they provide explicit analytic approximations rather than black-box numerical estimates. Yet until now, theoretical backing for these approximations was limited to the \(d^2/\lambda\to 0\) regime. Our results push these guarantees into substantially higher-dimensional regimes, arbitrarily close to the concentration threshold \(d/\lambda\to0\). See Section [ssec:statistics95motivation] for more details and references.

Besides physics and statistics, high-dimensional Laplace-type integrals arise throughout science and engineering, including in molecular simulation and theoretical chemistry (e.g., partition functions and free-energy calculations) [17], [18], Bayesian inverse problems [19], and others.

Finally, our result has interesting connections to the theory of cumulants and to normalizing flows in machine learning. The connection with the former is that we have solved the problem of bounding the remainder in cumulant expansions. Obtaining a high-dimensional remainder bound using cumulant theory (and the closely related theory of Gaussian chaos [20]–[22]) directly is deeply nontrivial and has not been done before. We have avoided this problem by finding an alternative route. See Section 1.3 for more details.

Regarding normalizing flows, these are sequences of maps \(T_1, \dots, T_L:\mathbb{R}^d\to\mathbb{R}^d\) with the property that if \(Z\sim\mathcal{N}(0, I_d)\) then \(T_1\circ\dots\circ T_L(Z)\) is approximately distributed according to a target distribution \(\pi\) [23], [24]. Typically, \(T_1,\dots,T_L\) are constructed using deep neural networks, the parameters of which are found by minimizing a loss [25]. Here, we approximate \(\pi\) in a similar fashion, pushing forward a Gaussian distribution through a series of transformations. But the key differences with normalizing flows are that 1) our construction is explicit, not through minimizing some loss, and 2) we prove a rigorous error bound. Essentially the only structural assumption we need to obtain this result is that the target \(\pi\) is a concentrating measure.

1.1 Application in physics↩︎

In a wide range of models in statistical physics and Euclidean QFT the partition function, which is defined as \[\label{eq:Z95field} Z=\int \exp\bigl(-\beta\,\mathcal{H}(\phi)\bigr)\,\mathcal{D}\phi,\tag{3}\] plays a central role. Here \(\beta\) is inversely proportional to the temperature, \(\mathcal{H}\) is the Hamiltonian of the model, and \(\int\dots\mathcal{D}\phi\) denotes functional (infinite-dimensional) integration over all allowed configurations \(\phi\) of the model [26], [27], [28]. The normalized logarithm of the partition function is known as free energy: \(\mathcal{F}=-\beta^{-1}\log Z\). The functional integral in 3 is typically interpreted as the limit of the appropriate finite dimensional integrals [27]. After discretization one obtains an ordinary (finite-dimensional) integral over \(d\) degrees of freedom, \[\label{eq:Z95Laplace} Z_{d}(\lambda)= \int_{\mathbb{R}^{d}} \exp\!\bigl(-\lambda f(x)\bigr)\,dx, \quad \mathcal{F}_{d}(\lambda)\;=\; -\lambda^{-1}\log Z_{d}(\lambda),\tag{4}\] where \(\lambda\) is a large prefactor coming from the problem [27]. If the model is discrete to begin with (i.e., it has finitely many degrees of freedom \(d\)), one gets 4 right away [27] and [26].

A basic approximation strategy, known as mean-field approximation, is to replace the integral in 4 by its leading Laplace contribution near a global minimizer \(x_\star\) of \(f\): \[\label{eq:mf} \log Z_{d}(\lambda) \approx -\lambda f(x_\star)-\tfrac12\log\det\bigl(\nabla^2 f(x_\star)\bigr)-\tfrac d2 \log (\lambda/(2\pi)).\tag{5}\] In physics texts, this step is typically justified by a combination of (i) the presence of a large parameter (e.g.volume), and (ii) an a posteriori self-consistency check that the discarded terms are “small” in the regime of interest; see again [27] for an example of this approach.

To go beyond 5 , one expands \(f\) around \(x_\star\), \[\label{eq:Taylor} f(x_\star + u) = f(x_\star)+\tfrac12\,u^\top H u +\sum_{k\ge 3}\tfrac1{k!}\,T_k[u^{\otimes k}], \quad H=\nabla^2 f(x_\star),\tag{6}\] rescales \(u=\lambda^{-1/2}y\), and rewrites 4 as a Gaussian expectation: Z_d()&=e^-f(x_)()^-d/2(H)^-1/2 ,
V_(Y)&:= _k T_k[Y^k], Y~(0,H^-1). Taking logs yields a formal cumulant expansion [15], [29]: \[\label{eq:cumulant} \log \mathbb{E}\bigl[e^{-V_\lambda(Y)}\bigr] = \sum_{m\ge 1}\frac{(-1)^m}{m!}\,\mathrm{cum}\!\bigl(V_\lambda(Y)^{\otimes m}\bigr),\tag{7}\] where \(\mathrm{cum}\) denotes the corresponding cumulant. In field-theoretic language, the terms in 7 are precisely the loop corrections to mean field: Wick expansion of the Gaussian moments [27], [28] produces a sum over Feynman diagrams, and the cumulant picks out the connected diagrams, which correct the free energy \(\mathcal{F}_{d}(\lambda)\) order-by-order in \(\lambda^{-1}\) [26], [28], [29].

Physics literature makes extensive use of truncations of 7 (or equivalent loop expansions). However, derivations are usually done without rigorous error control [27] and [26]. For finite \(d\), the error control follows from classical asymptotic analysis [1], [2]. For growing \(d\), such expansions have never been justified except for some specific cases.

In summary, mean field approximation and loop/cumulant expansions are central computational tools in statistical physics, but rigorous remainder estimates (especially in large systems with many degrees of freedom, \(d\gg1\)) are often absent from the standard presentations.

Our results provide a rigorous version of the loop-correction program for 4 in a joint limit where both the dimension \(d\) (number of effective degrees of freedom retained in the reduced description) and the Laplace prefactor \(\lambda\) grow. Concretely, under mild confining assumptions ensuring integrability and assuming a unique nondegenerate minimum \(x_\star\) together with bounded operator norms of derivatives near \(x_\star\), we obtain for every integer \(L\ge 1\) an expansion of \(\log Z_{d}(\lambda)\) through \(L-1\) loops whose remainder is explicitly controlled (see [main32res32intro]) Z_d() =& + _^L-1+_L(d,),
|_L(d,)| & d^L+1/^L, uniformly over the stated class of \(f\). The scaling condition \(d^{L+1}/\lambda^{L}\to 0\), \(d,\lambda\to\infty\) is exactly the statement that \((L-1)\)-loop-corrected mean field has a provable accuracy guarantee in a growing-system regime. From a physics viewpoint, [our95remainder] can be read as supplying the missing “error bars” for a procedure that is otherwise typically justified heuristically (or by numerics), thereby turning the loop expansion into a rigorously justified approximation scheme in regimes where the effective number of degrees of freedom grows with the large parameter in the exponential.

1.2 Application in statistics↩︎

For the sake of brevity, we focus on Laplace-type integrals and densities in the context of Bayesian statistics; however, these quantities frequently arise in frequentist statistics as well.

1.2.1 Statistics context and state of the art↩︎

Laplace-type integrals \(\int_{\mathbb{R}^d} g(x)e^{-\lambda f(x)}\mathrm{d}x\) and probability distributions \(\pi(x)\propto e^{-\lambda f(x)}\) are omnipresent in statistics. Here, \(\lambda\) (more commonly denoted \(n\) in statistics) plays the role of sample size, and \(f\) depends weakly on \(\lambda\) [30]. In Bayesian inference, \(\pi(x)=P(x\mid \mathrm{data})\) is the posterior probability distribution of an unknown parameter \(x\) given \(\mathrm{data}\), consisting of \(\lambda=n\) independent data points. Once the posterior has been specified, one is typically interested in (1) computing summary statistics, which take the general form \(\mathbb{E}_{X\sim\pi}[g(X)]\), (2) sampling from \(\pi\), and (3) computing the normalizing constant \(\int e^{-\lambda f(x)}\mathrm{d}x\) itself [31]. Summary statistics distill information in the posterior; samples enable exploration and uncertainty quantification. The normalizing constant, known as the model evidence [32], underpins Bayesian model selection. Here, one finds the best model by optimizing the evidence over candidate models [32], [33], [34]. All three tasks become particularly demanding in the growing-\(d\) regime that now arises routinely in applications.

The problem of how to do these computations efficiently has been actively studied for decades. A core approach is to replace \(\pi\) with a tractable approximation \(\hat{\pi}\), and the literature spans both theoretical and numerical methods; see, e.g., [30], [35]–[38] among many others. Within the theoretical analyses, a recent line of work focuses on the dependence of the total variation error \(\mathrm{TV}(\pi,\hat{\pi})\) on both \(d\) and \(\lambda\), either by sharpening error control for existing approximations \(\hat{\pi}\) or by constructing improved approximations. In particular, [14], [19], [39] showed that for the standard Laplace approximation \(\hat{\pi}=\mathcal{N}\!\big(x_\star,(\lambda\nabla^2 f(x_\star))^{-1}\big)\), the TV error scales as \(d\sqrt d/\sqrt\lambda\). Via a tighter analysis, [7], [8], [40] improved this dimension dependence to \(d/\sqrt\lambda\), and [41] showed that the closely related Gaussian variational-inference approximation also achieves TV error \(d/\sqrt\lambda\). Meanwhile, [40] and [37] proposed new approximations incorporating third-order derivative information, improving the \(\lambda\) dependence to \(\mathrm{TV}(\pi,\hat{\pi})\lesssim d^2/\lambda\) and \(\lesssim d^3/\lambda\), respectively.

While approximating \(\pi\) by a tractable \(\hat{\pi}\) is natural for sampling and for evaluating \(\mathbb{E}_{X\sim\pi}[g(X)]\) when \(g\) is nonsmooth, one can often do better when \(g\) is smooth by working directly with the LE. Indeed, \(\mathbb{E}_{X\sim\pi}[g(X)]\) is a ratio of two Laplace-type integrals. Expanding both the numerator and the denominator and then taking the ratio yields an explicit, fully deterministic approximation to \(\mathbb{E}_{X\sim\pi}[g(X)]\), avoiding the Monte Carlo error inherent in estimating expectations via samples from \(\hat{\pi}\). Despite this advantage, LEs have received comparatively little attention as a tool for computing posterior expectations. Notable exceptions are the fixed-\(d\) analyses in [30], [42].

In contrast, LE-based approximations are widely used for model selection problems that maximize the normalizing constant (model evidence) \(\int e^{-\lambda f(x)}\mathrm{d}x\) over a tuning parameter. A key attraction is that the LE provides an explicit analytic surrogate objective as a function of the tuning parameter, enabling efficient optimization. This stands in sharp contrast to methods that typically provide only black-box numerical access to the objective, such as Markov Chain Monte Carlo [43]. However, this convenience of the LE is meaningful only when the approximation error is controlled in a dimension-dependent manner. In particular, concerns about the reliability of the Bayesian Information Criterion (BIC) in high dimensions have been noted explicitly. The BIC introduces an additional approximation on top of the LE to further simplify the optimization objective. In [44], the authors write “for models with a large number of predictors [large \(d\)] it is no longer clear that some version of the BIC, or perhaps rather a Laplace approximation, accurately approximates a marginal likelihood”.

While recent work has begun to provide rigorous, dimension-dependent guarantees for high-dimensional LEs for normalizing constants [4], [45], [46], these require more restrictive scalings (at best, that \(d^2/\lambda\ll 1\)).

1.2.2 Implication of our results in statistics↩︎

Our work essentially completes the program of finding tractable analytic approximations of posterior densities, expectations, and normalizing constants in the high-dimensional but concentrating regime. We propose a combined approach for these related problems which allows dimension to grow relative to \(\lambda\) arbitrarily close to the concentration threshold.

For sampling, we construct a sequence of arbitrarily accurate approximations \(\hat{\pi}_L\) which can be sampled from using an explicit algorithm. For computing expectations of smooth functions, we give an arbitrarily accurate closed-form formula based on the LE of numerator and denominator. For normalizing constants, [main32res32intro] directly applies with \(g\equiv1\). The approximation error for all these tasks is \(d^{L+1}/\lambda^L\), which has the crucial feature that increasing \(L\) not only improves the accuracy but also expands the range of applicable \(d\), since we can take \(d\) as large as \(o(\lambda^{L/L+1})\).

We now give more detail about our methods.

Sampling and nonsmooth expectations.

To approximately sample from \(\pi\), the algorithm consists of drawing \(Z_i\) from \(\mathcal{N}(0,\lambda^{-1}I_d)\) and mapping it through one of the \(x_L\), \(L=1,2,3,\dots\), depending on the desired accuracy. To approximate expectations of nonsmooth functions \(g\), the samples \(X_i=x_L(Z_i) \sim \hat{\pi}_L\) can also be used, via \[\label{eq:monte-carlo-approx} \mathbb{E}_{X \sim \pi}[g(X)] \approx \mathbb{E}_{X \sim \hat{\pi}_L}[g(X)] \approx \frac{1}{N}\sum_{i=1}^{N} g(X_i).\tag{8}\]

Normalizing constants and smooth expectations.

Approximating normalizing constants is immediate using [main32res32intro] with \(g\equiv1\), so we don’t discuss it further. To approximate expectations \(\mathbb{E}_{X \sim \pi}[g(X)] = \int ge^{-\lambda f}/\int e^{-\lambda f}\) of smooth \(g\), we use [main32res32intro] for the numerator and denominator. This gives the approximation \(\mathbb{E}_{X \sim \pi}[g(X)]\approx\exp\bigl(\sum_{k=1}^{L-1}[b_k(f,g) - b_k(f,1)]/\lambda^k\bigr)\), which has accuracy \(\mathcal{O}(d^{L+1}/\lambda^L)\) uniformly over sufficiently smooth \(g\) with bounded derivatives near \(x_\star\).

This estimate improves on 8 in two ways. First, it does not incur any sampling error. Second, it is constructed using fewer derivatives of \(f\): \(2L-1\) versus \(2L+1\) for \(\hat{\pi}_L\). Intuitively, \(\hat{\pi}_L\) does not exploit smoothness of \(g\) and must compensate with more information from \(f\). In practice, this matters because \(f\) encodes the data likelihood and is expensive to differentiate, while the observable \(g\) is typically simple. For example, when \(L = 2\), the closed-form estimate \(\exp([b_1(f,g) - b_1(f,1)]/\lambda)\) has accuracy \(d^3/\lambda^2\) using only three derivatives of \(f\), whereas \(\hat{\pi}_1\) also uses three derivatives but achieves only the lower accuracy \(d^2/\lambda\). Thus using the closed-form asymptotic estimate improves the accuracy by a factor of \(d/\lambda\) over sampling approaches without requiring any more derivatives. We emphasize that most approaches in the literature are of the sampling type and achieve accuracy \(d^2/\lambda\) (our \(\hat{\pi}_1\) and [40]) or \(d^3/\lambda\) [37] at best, even disregarding the Monte Carlo error. See Section 9.4 for a detailed comparison to the two works closest to ours, [37], [40].

1.3 Connection to cumulant theory↩︎

As stated above, the coefficients \(b_k\) in [main32res32intro] can be derived from formal cumulant expansions. Although computing the coefficients using cumulants is easy, bounding the remainder in the expansion of \(\log I(\lambda)\) using cumulant theory seems deeply nontrivial. Indeed, cumulants have been thoroughly studied in the statistics, physics and combinatorics literature, yet rigorous remainder bounds on cumulant expansions in the growing \(d\) regime are virtually non-existent. This suggests that bounding the remainder through the lens of cumulants may be intractable. In contrast, our change-of-variables approach is not only tractable, but also avoids any heavy machinery, as mentioned above.

Our remainder bound makes rigorous the well-known observation that, based on the terms, expanding the cumulant generating function (cgf) is advantageous to expanding the moment generating function (mgf). Specifically, it is well-known that each cumulant is given by a sum of a fewer number of summands than the corresponding moment. For example at the end of Section 3.10.1 of [15], McCullagh writes that the formula for cumulants is a sum only over “connected pairs" of bi-partitions, as opposed to the formula for moments. In [27], the author writes,”When calculating cumulants, only fully connected diagrams (without disjoint pieces) need to be included. This is a tremendous simplification."

The simplification due to this summing over fewer terms in fact leads to a substantial improvement. By analyzing the cumulants and moments directly in our setting, it is possible to show that the moments that contribute to the \(k\)th coefficient in the \(I(\lambda)\) expansion are generically of order \(\mathcal{O}(d^{2k})\), while the cumulants that contribute to the \(k\)th coefficient in the \(\log I(\lambda)\) expansion are generically of order \(\mathcal{O}(d^{k+1})\). Of course, this analysis of the terms does not rigorously prove anything about the remainders. However, it naturally leads to the hypothesis that the remainders of the two expansions behave like the first terms dropped from the truncation, i.e. that the \(L\)th order remainders scale as \(\mathcal{O}(d^{2L}/\lambda^L)\) and \(\mathcal{O}(d^{L+1}/\lambda^L)\), respectively. This is precisely what the previous work [4] and the present paper rigorously prove.

Thus, as mentioned above, expanding the logarithm of \(I(\lambda)\) is what allows us to relax the \(d^2\ll\lambda\) requirement in [4]. Given our result [main32res32intro], it is very easy to see why \(d^2\ll\lambda\) is necessary for an additive expansion of \(I(\lambda)\) to any order. Indeed, we have shown in [main32res32intro] that, to leading order, \(I(\lambda)=\exp(\mathcal{O}(d^2/\lambda))\). Thus if we Taylor expand the exponential to obtain an additive expansion of \(I(\lambda)\), powers of \(d^2/\lambda\) are present throughout. There is no finite number of terms we can subtract from \(I(\lambda)\) to remove the dependence on \(d^2/\lambda\). However, we can remove the dependence on \(d^2/\lambda\) simply by dividing \(I(\lambda)\) by the single term \(\exp(b_1(f,g)/\lambda)\). The multiplicative remainder is then of order \(\exp(\mathcal{O}(d^3/\lambda^2))\).

In Section 9, we illustrate the power of expanding \(\log I(\lambda)\) in two concrete examples: a quartic perturbation of a Gaussian exponent, and a logistic-regression-type likelihood motivated by statistics.

1.4 Related work↩︎

There is a related body of work on Laplace-type integrals arising in statistics (e.g. posterior normalizing constants and marginals) and statistical physics (partition functions) in the “proportional asymptotics" regime. Here, \(d/\lambda\) converges to a constant. The regime lies beyond the concentration threshold, so these integrals can only formally be considered to be Laplace-type. Due to the lack of concentration, entirely different techniques are required to study such asymptotics. See [47]–[51], as well as the referenced cited therein.

Organization↩︎

The paper is organized as follows. Section 2 introduces the notation and conventions used throughout. In Section 3, we describe the problem setting and state the main result, and Section 4 outlines the proof strategy. Section 5 derives the expansion of the main contribution to the integral coming from a neighborhood of zero, while Section 6 computes the expansion coefficients in terms of cumulants. The tail contribution is estimated in Section 7. Section 8 presents our results on the approximation of \(\pi\propto e^{-\lambda f}\) and on computing expectations against \(\pi\). Section 9 compares our results with related work in the literature, and presents two applications of our expansion. Various technical results are collected in the appendices.

2 Notation and conventions↩︎

When we write e.g. \(\sum_{k\geq2}a_k\), the sum is over a finite number of \(k\)’s. Also, any sum of the form \(\sum_{k=1}^0 a_k\) is understood to be omitted.

In some contexts, the term tensor is used to refer to multilinear forms (the case \(j=0\)). Here we adopt the convention of Definition 1 and also call \(A_{k\to 1}\) and \(A_{k\to 2}\) tensors, i.e., vector- and matrix-valued multilinear maps.

The operator norm of a symmetric \(k\)-th order tensor \(A_k\), i.e., \(A_{k\to0}\), is given by [52] \[\label{T32norm} \Vert A_k\Vert_{\mathrm{op}}:=\sup_{\|u\|=1} A_k[u^{\otimes k}].\tag{9}\] For a function \(f\in C^k(\mathbb{R}^d)\), the \(k\)-th derivative tensor is \[(\nabla^kf(x))_{i_1\dots i_k} = \partial_{x_{i_1}}\dots\partial_{x_{i_k}}f(x).\]

We let \(\epsilon=d/\lambda\), and fix an integer \(L\ge 1\).

3 Statement of result↩︎

We study the asymptotics of the integral \[\label{main32int32v0} I(\lambda)=\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d}g(x)e^{-\lambda f(x)}\mathrm{d}x\tag{10}\] as \(d,\lambda\to\infty\). We make the following assumption on \(f\) and \(g\).

In 11 and everywhere below, recall that \(\epsilon=d/\lambda\).

The case of a general minimizer \(x_\star\) of \(f\) and general positive definite \(\nabla^2f(x_\star)\) can be reduced to the one assumed above by an affine change of variables. By linearity, the normalization condition \(g(0)=1\) is not restrictive. A sufficient condition for 11 to be satisfied is \[\begin{align} f(x)&\geq\kappa_f\min\left\{\|x\|^2/2, \sqrt{\epsilon}\left\|x\right\|, d^{-1/(2L)}\log(1+\left\|x\right\|)\right\},\quad x\in\mathbb{R}^d,\tag{12}\\ |g(x)|&\leq \min\left\{e^{\kappa_g\lambda\|x\|^2/2}, e^{\kappa_g\sqrt{d\lambda}\|x\|}, (1+\|x\|)^{\lambda\kappa_g/d^{1/(2L)}}\right\},\quad \|x\|\geq r_0,\tag{13}\\ \kappa&:=\kappa_f-\kappa_g>0. \end{align}\] When \(x\) is small, the minimum in the first line is equal to \(\min\{\|x\|^2/2, \sqrt{\epsilon}\left\|x\right\|\}\). When \(x\) is large, the minimum is given by \(d^{-1/(2L)}\log(1+\left\|x\right\|)\). Thus we only require that \(f\) grows logarithmically at infinity. Meanwhile, \(g\) can grow at most polynomially at infinity under 13 , but the power of the polynomial can become arbitrarily large as \(\lambda\to\infty\).
For positive scalars \(a,b\), the notation \(a\lesssim b\) and \(a=\mathcal{O}(b)\) both mean that \(a\le cb\), where the suppressed constant \(c\) can depend on the constants \(c_{f}\) and \(c_{g}\) appearing in Assumption 1 but not on \(d,\lambda,r_0,\kappa\). Since \(L\ge 1\) is fixed, the dependence of constants on \(L\) is ignored. In a similar vein, when we say “for all \(\delta\) sufficiently small, it holds.…” we mean there exists \(c(c_{f},c_{g})\) such that if \(0<\delta\leq c(c_{f},c_{g})\), then the statement after the ellipses is true. Here, \(\delta\) denotes an arbitrary small parameter that is relevant to the discussion at hand. In most cases, \(\delta=R\sqrt\epsilon\) (see below for the definition of \(R\)).

The following is our first main result.

To illustrate the convention introduced before Theorem 1, we note that the suppressed constant in [main32res] and in the bound on \(|b_k(f,g)|\) may depend on \(c_{f}\), \(c_{g}\) only, but not on \(d,\lambda,r_0,\kappa\). Note also that, according to our convention, the sum with respect to \(k\) in [main32res] is omitted if \(L=1\).

We give a closed-form formula for \(b_1(f,g)\) in Section 6.2. As is seen, in the worst case the first inequality in 14 forces \(R\sim\log\lambda\), up to a factor depending on \(c_{f},c_{g},\kappa\). This occurs when \(d\) remains bounded as \(\lambda\to\infty\). On the other hand, if \(d\geq c(c_{f},c_{g},\kappa)\log\lambda\), then 14 is satisfied with \(R\) being a constant independent of \(d\) and \(\lambda\).

4 Proof outline↩︎

4.1 Key steps↩︎

Here, we outline the main ideas of our approach. Without loss of generality, we suppose \(x_\star=0\) is the global minimizer of \(f\), with \(\nabla^2f(0)=I_d\), so that \(f(x)=\|x\|^2/2+o(\|x\|^2)\) near zero. For simplicity, we only describe the proof in the case \(g\equiv1\). Incorporating a non-constant \(g\) presents no real challenges. On the surface, our proof begins similarly to classical fixed-\(d\) proofs of the expansion based on the Morse Lemma [2]. This lemma states that there is a change of variables \(x=x(t)\) making the exponent \(f(x(t))\) an exact quadratic. This approach, in its pure form, seems to be intractable because the coordinate transformation from the Morse Lemma is difficult to work with.
Instead, we use a variation of the approach. Fix any \(L\ge1\).
Step 1: initial change of variables. We begin by constructing an explicit local polynomial change of variables \(X(t)\), of the form \(X(t)=t+\mathcal{O}(\|t\|^2)\), \(t\to0\), which makes the exponent “more quadratic” but not exactly quadratic. Specifically, it eliminates the third through \((2L+1)\)st order terms in the Taylor expansion of \(f\) around the minimizer, so that \(-\lambda f(X(t))=-\lambda\left\|t\right\|^2/2+\mathcal{O}(\lambda\left\|t\right\|^{2L+2})\). But \(\mathcal{O}(\lambda\left\|t\right\|^{2L+2})\) is of order \(O(d^{L+1}/\lambda^L)\) (and therefore negligible) in the region \(\|t\|\lesssim\sqrt{d/\lambda}\) where the integral concentrates. Thus \(\mathcal{O}(\lambda\left\|t\right\|^{2L+2})\) can be discarded.

The price we pay for this nice change of variables is the appearance of the Jacobian \(\mathrm{det}(X^\prime(t))\) of the coordinate change, which we bring into the exponent. Thus the new exponential function is, upon throwing out the negligible \(\|t\|^{2L+2}\) term, given by \(\exp(-\lambda\|t\|^2/2 + \log\mathrm{det}(X^\prime(t)))\). Crucially, however, \(\log\mathrm{det}(X^\prime(t))\) scales only as \(d\ll \lambda\). Thus the log-Jacobian does not significantly affect the quadratic exponent created by the change of variables. Specifically, we may write \(\log\mathrm{det}(X^\prime(t))=dh(t)\), where \(h(t)\) can be Taylor-expanded as \(h(t)=\sum_{k\geq1}F_k[t^{\otimes k}]+(\text{negligible remainder})\) for some tensors (multi-linear forms) \(F_k\) with bounded operator norms. Here and below, sums in which the upper limit has not been explicitly indicated are understood to mean finite sums. The above arguments lead to I()= & e^(d^L+1/^L)()^d/2_U_1 (-E_1(t))t,
E_1(t):=&t̍̍^2+_kF_k[t^k], where \(\epsilon=d/\lambda\), \(\mathcal{U}_1=\{\|t\|\lesssim\sqrt{\epsilon}\}\), and we have discarded the negligible integral over \(\mathcal{U}_1^c\). Now, we must somehow deal with the terms \(\epsilon F_k[t^{\otimes k}]\), \(k\geq1\). The linear and quadratic terms (\(k=1,2\)) are not a problem, since they can be combined with \(\|t\|^2/2\) to give a new Gaussian measure, and integrating against Gaussians is tractable. For \(k\ge 2L\), the terms \(\lambda\epsilon|F_k[t^{\otimes k}]|\) are sufficiently small uniformly over \(t\in\mathcal{U}_1\), so that \(\exp(-\lambda\epsilon F_k[t^{\otimes k}])\) can be discarded from the integral (meaning, absorbed in the multiplicative remainder \(e^{\mathcal{O}(d^{L+1}/\lambda^L)}\)). But there is an intermediate range of \(k\)’s, namely, \(3\le k\le 2L-1\), for which additional massaging is needed. We do this in the next step. In what follows, all our formulas ignore terms that lead to multiplicative factors \(e^{\mathcal{O}(d^{L+1}/\lambda^L)}\) in \(I(\lambda)\).
Step 2: iterative refinement. We construct a polynomial change of variables of the form \(T_1(s)=s+\epsilon\varphi_1(s)\), with \(\varphi_1(s)=\mathcal{O}(\|s\|^2)\), analogous to \(X(t)\) above. By choosing \(\varphi_1\) appropriately, we can ensure that this change of variables increases the power of \(\epsilon\) in front of \(F_k[t^{\otimes k}]\) for as many \(k\geq3\) as we like. Thus \[\label{E1intro} E_1(T_1(s))=\epsilon s+\epsilon s^2+\tfrac12\|s\|^2+\epsilon^2\sum_{k= 3}^{2L-3}s^k.\tag{15}\] where \(s^k\) is a mnemonic form for \(F_k[s^{\otimes k}]\), and \(F_k\) is some tensor with bounded operator norm. Furthermore, due to the form of \(T_1(s)\), we have \[\label{1laintro} \tfrac1\lambda\log\det(T_1^\prime(s))=\epsilon^2\sum_{k=2}^{2L-3}s^k,\tag{16}\] informally. Here, one of the powers of \(\epsilon\) comes from \(\epsilon\varphi_1(s)\), and the other power of \(\epsilon\) comes from \(\frac{1}{\lambda}\times d\), since \(\log\det\sim d\). Since the powers of \(\epsilon\) in 16 are greater than or equal to the corresponding powers of \(\epsilon\) in 15 , adding \(\frac{1}{\lambda}\log\det(T_1^\prime(s))\) to \(E_1(t_1(s))\) does not change the structure of the latter.

Therefore, implementing the change of variables, the above argument gives I()= & e^(d^L+1/^L)()^d/2_U_2 (-E_2(t))t,
E_2(t):=&t+t^2+t̍̍^2+^2_k=3^2L-3t^k, U_2:=T_1^-1(U_1). We can now repeat the procedure, each time increasing the power of \(\epsilon\) in front of \(\sum_{k\geq3}t^k\). Finally, we arrive at \(E_L(t)=\epsilon t+\epsilon t^2+\frac{1}{2}\|t\|^2\) in the exponent. At this point, all terms in the sum with respect to \(k\) have been absorbed in the multiplicative remainder \(e^{\mathcal{O}(d^{L+1}/\lambda^L)}\). Recalling that \(\epsilon t+\epsilon t^2\) is a mnemonic for \(\epsilon a^\top t+ \epsilon t^\top Bt\), where the norms of \(a\) and \(B\) are bounded, we conclude that I()=& e^(d^L+1/^L)()^d/2_U_L (-),
U_L:=&T_L-1^-1…T_1^-1(U_1).

Step 3: complete the square. We incur a negligible error by extending the integral over \(\mathcal{U}_L\) in [ELintro] to \(\mathbb{R}^d\). But the integral over \(\mathbb{R}^d\) is now exactly computable, since it becomes a Gaussian integral upon completing the square. We obtain precisely that \[I_{\mathrm{in}}(\lambda)=e^{\mathcal{O}(d^{L+1}/\lambda^L)}\exp\left(\tfrac12\lambda\epsilon^2a^\top (I_d+2\epsilon B)^{-1}a-\tfrac12\log\det (I_d+2\epsilon B)\right).\] Finally, the second exponential factor above gives the terms in [main32res32intro]. To show this, we prove that \(a\) can be written as \(a=a_0+\epsilon a_1+\epsilon^2a_2+\dots\), and \(B\) as \(B=B_0+\epsilon B_1+\epsilon^2B_2+\dots\), where the \(a_k\) and \(B_k\) do not depend on \(\epsilon\), and have bounded norms. Substituting this into the exponent in parentheses leads to \(\exp(\lambda\sum_{k\geq2}c_k\epsilon^k +d\sum_{k\geq1}c_k'\epsilon^k)\), where the first sum stems from \(\tfrac12\lambda\epsilon^2a^\top (I_d+2\epsilon B)^{-1}a\) and the second sum stems from \(-\tfrac12\log\det (I_d+2\epsilon B)\). Combining the two sums gives the right orders of magnitude \(d^2/\lambda+d^3/\lambda^2+\dots\) for the terms in [main32res32intro].

It is worth noting that our proof avoids the heavy Gaussian concentration machinery for Lipschitz functions (e.g., via log-Sobolev/Herbst-type arguments), in contrast to the strongest result to date [4].

Step 4: computing the coefficients. In Step 3 above, we have quantified the orders of magnitude of the expansion coefficients. Next, we show that these coefficients can be expressed in terms of cumulants. To do so, it suffices to consider \(f,d\) fixed, view \(\log I(\lambda)\) purely as a function of \(t=\lambda^{-1/2}\), and take the derivatives of this function with respect to \(t\). As a first step towards this goal, we write \(I(\lambda)\) in the following form: \[I(\lambda)=e^{\mathcal{O}(t^{2L})}\mathbb{E}\left[\exp\left(tf_3(Z)+t^2f_4(Z)+\dots+t^{2L-1}f_{2L+1}(Z)\right)\mathbb{1}\{\|Z\|\leq C\log(1/t)\}\right],\]where \(Z\sim\mathcal{N}(0, I_d)\) and \(f_k(x)=-\frac{1}{k!}\nabla^kf(0)[x^{\otimes k}]\). Using this representation, we expand \(\log I(\lambda)\) in powers of \(t\). We show that the coefficients up to order \(2L\) in this expansion coincide with the coefficients in the formal expansion of \(\log \mathbb{E}[\exp(\sum_{k=1}^{2L-1}t^kf_{k+2}(Z))]\) in powers of \(t\). This latter expansion is formal since the expectation may not be finite. The reason the coefficients coincide is that the discrepancy between the coefficients of the two expansions is controlled by \(P(\|Z||\geq C\log(1/t))\sim\exp(-\log^2(1/t))=o(t^M)\) for any \(M\). Finally, we observe that the coefficients of the latter formal expansion can be expressed in terms of joint cumulants of \(f_k(Z)\), \(k=3,\dots,2L+1\).

4.2 Roadmap↩︎

In this section, we make the above proof outline precise, by breaking up the proof of Theorem 1 into three lemmas. These lemmas are stated here and proved in the subsequent sections, as indicated below the statement of each lemma.

We break up the proof of Theorem 1 into three key lemmas.

Lemma 1 is proved in Section 5. Let \(\alpha=(\alpha_1,\dots,\alpha_M)\) be a multiindex with \(\alpha_j\geq 0\) for all \(j\). Let \(|\alpha|=\alpha_1+\dots+\alpha_M\), and \(\partial_u^\alpha=\partial_{u_1}^{\alpha_1}\dots\partial_{u_M}^{\alpha_M}\).

Lemma 2 is proved in Section 6. Although we have used somewhat different notation, the above formula can be shown to coincide with that given in [15].

Lemma 3 is proved in Section 7. The three lemmas conclude the proof of Theorem 1. Indeed, they imply \(I(\lambda)>0\) when \(R^{2L}d^{L+1}/\lambda^L\) is sufficiently small, give the desired form for the coefficients, and imply \[\bigg|\log I(\lambda)-\sum_{k=1}^{L-1}b_k(f,g)\lambda^{-k}\bigg|\leq \bigg|\log I_{\mathrm{in}}(\lambda)-\sum_{k=1}^{L-1}b_k(f,g)\lambda^{-k}\bigg| + \left|\log\left(1+\frac{I_{\mathrm{out}}(\lambda)}{I_{\mathrm{in}}(\lambda)}\right)\right|.\]We now use 20 , 21 to conclude [main32res].

5 Proof of Lemma 1↩︎

Using Assumption 1, part [Ial], and using that \(R\sqrt\epsilon\leq r_0\), we have \[\max_{\|x\|\leq R\sqrt\epsilon}\lambda|\mathcal{R}_{2L+2}(x)| \lesssim\lambda(R^2\epsilon)^{L+1}=R^{2L+2}d\epsilon^L.\] Therefore, recalling [eq:Taylor-f], we have \[\label{I-in-1} I_{\mathrm{in}}(\lambda)=e^{\mathcal{O}(R^{2L+2}d\epsilon^L)}\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\|x\|\leq R\sqrt\epsilon}g(x)e^{-\lambda f_{2L+1}(x)}\mathrm{d}x.\tag{22}\] The integrand must be positive for 22 to hold. But indeed, \(g(0)=1\) and \(\log g\) has bounded first derivative, both by [Ial32g] of Assumption 1. Thus we can ensure \(g(x)>0\) for all \(\|x\|\leq R\sqrt\epsilon\) by choosing the upper bound \(c^+\) on \(R\sqrt\epsilon\) small enough in 19 . This proves our claim in Lemma 1 that \(I_{\mathrm{in}}(\lambda)>0\).

5.1 Iterative change of variables↩︎

Throughout this section, we assume all the conditions of Lemma 1 hold true without explicitly saying so. Before stating a key change of variables, we introduce the following two definitions.

Any time we write \(F_{k\to j}^{\epsilon}\) (or \(F_{k}^{\epsilon}\) when \(j=0\)) we mean a composite tensor according to this definition.

See Appendix 11 for the proof and construction of \(\varphi\). As an example, when \(L=1\), we take \(X(t)=t-\frac{1}{6}\nabla^3f(0)[t^{\otimes2}]\), where \(\nabla^3f(0)[t^{\otimes2}]\) is the vector such that \(\nabla^3f(0)[t^{\otimes2}]^\top u=\nabla^3f(0)[t,t,u]\) for all \(u\in\mathbb{R}^d\). It is then straightforward to check that _3(t-^3f(0)[t^])&=t̍-^3f(0)[t^]̍^2+^3f(0))^]
&=t̍̍^2/2+ (t̍̍^4). The basic observation that substituting \(x=t-F[t^{\otimes k-1}]\) into \(\|x\|^2/2+F[x^{\otimes k}]\) kills the order \(k\) polynomial is at the heart of all of our changes of variables.
We now show that \(X\) is bijective and characterize the set \(X^{-1}(\{\|x\|\leq R\sqrt\epsilon\})\). The following lemma states a slightly more general result which will be needed later on.

See Appendix 11 for the proof. Let \[\label{usus1} \mathcal{U}=\{\|x\|\leq R\sqrt\epsilon\},\quad \mathcal{U}_1=X^{-1}(\{\|x\|\leq R\sqrt\epsilon\}).\tag{23}\] Lemma 5 with \(C_1=C_2=1\), \(r=R\sqrt\epsilon\), and \(\mathcal{U}\) as above gives that if \(R\sqrt\epsilon\) is small enough then \(\{\|t\|\leq \frac{2}{3}R\sqrt\epsilon\}\subseteq\mathcal{U}_1\subseteq\{\|t\|\leq 2R\sqrt\epsilon\}\). Thus in particular, the conclusion of Lemma 4 holds for all \(t\in\mathcal{U}_1\). Combining 22 , Lemma 4, and the fact that \(X(t):X^{-1}(\mathcal{U})\to\mathcal{U}\) is a bijection by Lemma 5, we have \[\label{initchange} I_{\mathrm{in}}(\lambda)=e^{\mathcal{O}(R^{2L+2}d\epsilon^L)}\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathcal{U}_1}\exp\left(-\frac{\lambda}{2}\|t\|^2+\log g(X(t))+\log \mathrm{det}(X'(t))\right)\mathrm{d}t.\tag{24}\] Since \(\|\varphi'(t)\|_\mathrm{op}\ll 1\) for \(t\in\mathcal{U}_1\), we have \(\log\mathrm{det}X'(t)=\log \mathrm{det}( I_d+\varphi'(t))=\mathrm{tr}\log(I_d+\varphi'(t))\approx\mathrm{tr}\varphi'(t)\sim d\). This is made precise in the following lemma.

The proof follows from Lemma 16 with \(m=0\), and the fact that \(R\sqrt\epsilon\) can be made sufficiently small by choosing \(c^+\) appropriately in 19 . Next, we expand \(\log g(X(t))\).

We now use [logdet32v3] and [log-g] in 24 to get t̍̍^2&-g(X(t))-( X’(t))
&= (t̍̍^2+_k=1^2L-1F_k^[t^k]+_k=1^2L-1F_k^[t^k])+(t̍̍^2L),tU_1. We can combine \(\frac{1}{\lambda} F_{k}^{\epsilon}\) and \(\epsilon F_{k}^{\epsilon}\) into \(\epsilon F_{k}^{\epsilon}\). Also, since \(\|t\|\leq 2R\sqrt\epsilon\) for all \(t\in\mathcal{U}_1\) as discussed below 23 , we have \(\mathcal{O}\big(\left\|t\right\|^{2L}\big)=\mathcal{O}(R^{2L}\epsilon^L)\). We conclude that I_()&=e^(R^2L+2d^L)()^d/2_U_1(-E_1(t))t,
E_1(t)&:=t̍̍^2+_k=1^2L-1F_k^[t^k],
{t̍̍23R}&U_1{t̍̍2R}. Comparing \(E_1\) in [E1] with the original \(f\), we see that \(E_1\) is closer to being exactly quadratic. If \(L=1\), we stop here and estimate \(I_{\mathrm{in}}(\lambda)\) as in Section 5.2. If \(L\ge 2\), we show next that by iteratively changing variables, we can continue to increase the power of \(\epsilon\) in front of cubic and higher powers of \(t\).

See the end of Appendix 11 for the proof and the construction of \(T_m\). The value of \(C\) will be chosen below the proof of Corollary 1.

As we will see below, the big-\(\mathcal{O}\) term in 26 contributes the factor \(\exp\big(\mathcal{O}(\lambda(R\sqrt\epsilon)^{2L+2})\big)\) to \(I_{\mathrm{in}}\). Observe that this matches the desired order of magnitude in Lemma 1. Setting \(m=1\) in 25 gives the function \(E_1\) in [E1].

Comparing 25 to 26 , we see that the power of \(\epsilon\) and highest power \(k\) of \(s\) have changed. The effect of \(T_m(s)\) is only to increase the power of \(\epsilon\) from \(m\) to \(m+1\), not to decrease the highest power of \(s\) from \(2L-2m+1\) to \(2L-2m-1\). Terms with these higher powers of \(s\) do appear upon plugging in \(T_m(s)\) to \(E_m\). However, precisely because of the higher power of \(\epsilon\), any \(\epsilon^{m+1}F_{k}^{\epsilon}[s^{\otimes k}]\) with \(k>2L-2m-1\) can simply be thrown out, i.e. absorbed into the \(\mathcal{O}((R\sqrt\epsilon)^{2L+2})\) remainder.
In the following lemma, we study the change of variables \(T_m\).

The statements (1) and (2) follow from Lemma 5 with \(r=R\sqrt\epsilon\). Part (3) follows by setting \(N=2L-2m\) in Lemma 16 and using part (1) to characterize the diameter of \(\mathcal{U}_{m+1}\). We combine the above two lemmas in the following corollary, thereby completing one full iteration from \(E_m\) to \(E_{m+1}\).

Starting with [E1], we iteratively apply Lemma 9 and Corollary 1, stopping once we get to \(m+1=L\). We conclude that I_()&=e^(R^2L+2d^L)()^d/2_U_Le^-E_L(t)t,
E_L(t)&= F_1^[t]+F_2^[t^]+t̍̍^2, U_L=T_L-1^-1…T_1^-1(U_1). Furthermore, by the repeated application of Lemma 9, and using that \(c_1=2/3, C_1=2\) (recalling [E1]), we have \[\label{UL32size} \{\|t\|\leq (2/3)^LR\sqrt \epsilon\}\subseteq \mathcal{U}_{L}\subseteq \{\|t\|\leq 2^LR\sqrt\epsilon\}.\tag{28}\] Thus we see that in Lemma 8, we can take \(C=2^L\). Note that [IJ123] and 28 both also hold for \(L=1\), i.e. \(E_L\) from [IJ123] coincides in structure with \(E_1\) from [E1], upon setting \(F_{2}^{\epsilon}=0\) in [IJ123].

The following quantities will appear in Sections 5.2, 5.3, and 8.

Since \(F_{2}^{\epsilon}\) is by definition a symmetric bilinear form, the matrix \(B\) is symmetric.

5.2 Completing the square↩︎

Before presenting the main result of the section, we make the following observation: let \(A=F_{2}^{\epsilon}=F_{2\to 0}^{\epsilon}\) be a composite tensor given by a bilinear form, i.e. taking in two vectors and returning a scalar. Then \(A\) can also be viewed as \(A=F_{0\to 2}^{\epsilon}\), i.e. a constant mapping returning a matrix. The same is true for \(F_{1\to 0}^{\epsilon}\) also being \(F_{0\to 1}^{\epsilon}\). The operator norms are preserved under this identification.
The main result in this section is the following.

5.3 Structure of terms in the expansion↩︎

Recall from [aB] that \(a=\epsilon J_1\) and \(B=I_d+2\epsilon J_2\), and consider the two terms in 30 . We have \(\log\det(I_d+2\epsilon J_2)=\mathrm{tr}\log(I_d+2\epsilon J_2)\). Using that \(\|2\epsilon J_2\|_\mathrm{op}<1/2\) if \(\epsilon\) is sufficiently small (since \(J_2=F_{0\to 2}^{\epsilon}\)), we have (̍I_d+2J_2)^-1- _k=0^L-2 (- 2J_2)^k̍_&^L-1,
̍(I_d+2J_2) + _k=1^L-1 1k(-2J_2)^k̍_&^L. By Lemma 15 (see also Remark 9) the first sum is \(F_{0\to 2}^{\epsilon}\) and the second sum is \(\epsilon F_{0\to 2}^{\epsilon}\). Therefore, _1^(I_d+2J_2)^-1J_1=&J_1^F_0^J_1+(^L-1),
(I_d+2J_2)=&F_0^ +(d^L). From the Definition 3 and since \(J_1=F_{0\to 1}^{\epsilon}\), it is clear that \(J_1^\top F_{0\to 2}^{\epsilon}J_1=F_{0\to 0}^{\epsilon}=:\mu\), where \(\mu\) is just a polynomial in \(\epsilon\), i.e. \(\mu=\sum_{k\ge0}\mu_k\epsilon^k\). Its coefficients are bounded and depend only on \(d\), \(\nabla^kf(0)\), and \(\nabla^k(\log g)(0)\). We conclude that ^2J_1^&(I_d+2J_2)^-1J_1-(I_d+2J_2)
&=^2+F_0^ +(d^L). Finally, write \(F_{0\to 2}^{\epsilon}=\sum_{\ell\geq0}\epsilon^\ell B_\ell\), for bounded matrices \(B_\ell\) depending only on \(d\), \(\nabla^kf(0)\), and \(\nabla^k(\log g)(0)\). Let \(a_k=d^{-1}\mathrm{tr}B_{k-1}\), which are bounded. Then ^2+F_0^&=^2 +_^(B_)+(d^L)
&=_k=2^L _k^k +(^L+1)+_^L-2^(B_)+(d^L)
&=_k=2^L_k^k+_k=1^L-1da_k^k+(d^L)=_k=1^L-1(_k+1+a_k)d^k+(d^L). To get the second line, we moved the part of the polynomial \(\mu=\sum_{k\ge0}\mu_k\epsilon^k\) in which \(k\geq L-1\) into the remainder. Similarly, we moved the part of the second sum in which \(\ell\geq L-1\) into the remainder.

Let \(b_k(f,g,d)=(\mu_{k+1}+a_k)d^{k+1}\), so that \(|b_k(f,g,d)|\lesssim d^{k+1}\). Combining [lae] with [laef] and using the result in 30 , we finally have I_()=_k=1^L-1b_k(f,g,d)^-k+(R^2L+2d^L) This concludes the proof of Lemma 1.

6 Cumulants and proof of Lemma 2↩︎

In this section, we prove Lemma 2. First, we introduce the concept of cumulants.

6.1 Preliminaries on cumulants↩︎

When there is repetition among the \(Y_j\), there is an alternative expression for the cumulant. Specifically, suppose we have another set of random variables \(X_k\), \(k=1,2,3,\dots\). Let \(\alpha=(\alpha_1,\dots,\alpha_M)\) with \(\alpha_j\geq0\). Suppose that \(Y_1,\dots, Y_m\) consist of \(\alpha_1\) copies of \(X_1\), \(\alpha_2\) copies of \(X_2\), and so on, up to \(\alpha_M\) copies of \(X_M\). Let \(|\alpha|=\alpha_1+\dots+\alpha_M=m\). Then (Y_1,…,Y_m)&=(__1,__2,…,__M)
&=(-i)^||_u_1^_1…_u_M^_M )]_u_1=…=u_M=0. This follows from the fact that, for a function \(f(s_1,s_2)=g(s_1+s_2)\), it holds \(\partial_{s_1}\partial_{s_2}f(0,0) =g''(0)\). We will only explicitly use cumulants of order one and two, for which it holds (X)&=[X],
(X,Y)&=(X,Y)=[XY]-[X][Y]. See [21]. Let \[\label{pkdef} p_k(x)=\nabla^k(\log g)(0)[x^{\otimes k}] -\frac{1}{(k+1)(k+2)}\nabla^{k+2}f(0)[x^{\otimes k+2}],\quad k=1,2,3,\dots,\tag{33}\] be functions on \(\mathbb{R}^d\). For example, for \(p_1,p_2\), we have p_1(x) &= g(0)[x] -^3f(0)[x^],
p_2(x) &= (^2 g(0)-g(0)^)[x^] -^4f(0)[ x^].

Let \(X\in\mathbb{R}^d\) be a random variable. Recall that \(\alpha\) is a multiindex. We define (p_(X))&=(__1,__2,…,__M)
&=(-i)^||_u^)]_u_1=…=u_M=0. Here, the second line is by [cum32id]. Thus, for example, \[\mathrm{cum}(p_{3,1,0,1}(X))=\mathrm{cum}(p_1(X),p_1(X),p_1(X),p_2(X),p_4(X)).\] For convenience, we recall the formula in Lemma 2, which we will prove in Section 6.3. \[\label{eq:c} b_{\frac{M}{2}}(f,g,d)=b_{\frac{M}{2}}(f,g)=\sum_{\substack{\alpha_1,\dots,\alpha_M\geq0\\\sum_{i=1}^Mi\alpha_i=M}}\frac{\mathrm{cum}(p_\alpha(Z))}{\prod_{i=1}^M\alpha_i!(i!)^{\alpha_i}}.\tag{34}\] The fact that there is no explicit dependence on \(d\) is clear from the righthand formula. Indeed, we see from [cumpalph] that in \(\mathrm{cum}(p_\alpha(Z))\), the only appearance of \(d\) is in the functions \(p_1(Z),\dots,p_M(Z)\). But we see in 33 that the definition of \(p_1,\dots,p_M\) is agnostic to the number of arguments the functions \(f\) and \(g\) have.

6.2 More on the terms \(b_k\)↩︎

We first study \(b_{M/2}\) to determine the highest-order \(f\) derivative contributing to it. This is useful in Section 8 below, where we study the derivative order of \(b_{L-1}(f,g)-b_{L-1}(f,1)\). We then compute \(b_1(f,g)\) explicitly.

By 33 , the highest-order \(f\) derivatives in the formula 34 for \(b_{M/2}\) necessarily arise from \(\alpha\) such that \(\alpha_M>0\). But since \(\sum_{i=1}^Mi\alpha_i=M\), if \(\alpha_M>0\) then we must have \(\alpha_M=1\) and \(\alpha_i=0\) for all \(i\neq M\). Thus the highest-order \(f\) derivative contribution to \(b_{M/2}(f,g)\) is &==
&=-.

Next, we use 34 with \(M=2\) to compute \(b_1(f,g)\). Only \(\alpha=(2,0)\) and \(\alpha=(0,1)\) satisfy \(\alpha_1+2\alpha_2=2\). Thus \[\label{c1fg} b_1(f,g)=\frac{\mathrm{cum}(p_1(Z),p_1(Z))}{2!(1!)^{2}}+ \frac{\mathrm{cum}(p_2(Z))}{1!(2!)^{1}}= \frac{1}{2}\mathrm{Var}(p_1(Z)) + \frac{1}{2}\mathbb{E}[p_2(Z)],\tag{35}\] using [cum12]. Recall \(p_1,p_2\) from [p1p2]. Let \[T_1=\nabla g(0),\quad T_2 =\nabla^2g(0)-\nabla g(0)^{\otimes2},\quad T_3=-\nabla^3f(0)/6,\quad T_4=-\nabla^4f(0)/12.\]Then \(p_1(x)=T_1[x]+T_3[x^{\otimes3}]\) and \(p_2(x)=T_2[x^{\otimes2}]+T_4[x^{\otimes4}]\). The polynomial \(p_1\) is odd, therefore \(\mathbb{E}[p_1(Z)]=0\) and \(\mathrm{Var}(p_1(Z))=\mathbb{E}[p_1(Z)^2]\). We now rewrite \(p_1\) using Hermite polynomials. For \(x=(x_1,\dots,x_d)\) define the first- and third-order multivariate Hermite polynomials [53] \[H_i(x)=x_i,\qquad H_{ijk}(x) =x_i x_j x_k - x_i\mathbb{1}\{j=k\} - x_j\mathbb{1}\{i=k\} - x_k\mathbb{1}\{i=j\} .\] Thus \(H_{iii}(x)=x_i^3-3x_i\), \(H_{iij}(x)=(x_i^2-1)x_j\) for \(i\neq j\), and \(H_{ijk}(x)=x_i x_j x_k\) when \(i,j,k\) are distinct. Reordering indices in the subscript does not change the polynomial. It is straightforward to check that we may write \[p_1(x) = \sum_{i=1}^d\bigg(T_1^i+3\sum_{j=1}^dT_3^{ijj}\bigg)H_i(x) + \sum_{i,j,k=1}^dT_3^{ijk}H_{ijk}(x).\] The Hermite polynomials are orthogonal in \(L^2(\mathcal{N}(0,I_d))\) [53]. In other words, if \(Z\sim\mathcal{N}(0, I_d)\), then \(\mathbb{E}[H_i(Z)H_j(Z)]=0\) if \(i\neq j\), \(\mathbb{E}[H_i(Z)H_{jk\ell}(Z)]=0\) for any \(i,j,k,\ell\), and \(\mathbb{E}[H_{ijk}(Z)H_{\ell mn}(Z)]=0\) if \((i,j,k)\neq (\ell,m,n)\), viewed as unordered triplets. Furthermore, we have \(\mathbb{E}[H(Z_i)^2]=1\) and it is straightforward to show \(\mathbb{E}[H_{ijk}(Z)^2]=1=6/3!\) if \(i,j,k\) are distinct, \(\mathbb{E}[H_{iij}(Z)^2]=2=6/3\) if \(i,j\) are distinct, and \(\mathbb{E}[H_{iii}(Z)^2]=6=6/1\). Thus for general \(i,j,k\), the expectation \(\mathbb{E}[H_{ijk}(Z)^2]\) is given by \(6\) divided by the number of distinct ways to rearrange the indices \(i,j,k\).

Using these facts, and the symmetry of the tensor \(T_3\), we conclude that (p_1(Z))&=[p_1(Z)^2] = _i=1^d(T_1^i+3_j=1^dT_3^ijj)^2 +6T̍_3̍_F^2
=&̍g(0)-f(0)̍^2+̍^3f(0)̍_F^2. Next, a straightforward calculation gives [p_2(Z)] = &_i=1^dT_2^ii +3_i,j=1^dT_4^iijj =g(0)-̍g(0)̍^2 -^2f(0). We substitute [p1] and [p2] in 35 to get b_1(f,g) =&-f(0)^g(0)+ ̍f(0)̍^2+̍^3f(0)̍_F^2+g(0)-^2f(0).

6.3 Proof of Lemma 2↩︎

Let \(t=\lambda^{-1/2}\), and define \[\label{I032t} I_0(f, g, d, t)=\log I_{\mathrm{in}}(\lambda) = \log\left\{\Big(\frac{t^{-2}}{2\pi}\Big)^{d/2}\int_{\|x\|\leq t\log(1/t)\sqrt d} g(x)e^{-t^{-2} f(x)}\,\mathrm{d}x\right\}.\tag{36}\] Here, note that we have chosen \(R=R(t)=\log(1/t)\). Let \(c^\pm\) be as in Lemma 1. Then this lemma gives that for all \(f,g\) satisfying parts [Ial], [Ial32g] of Assumption 1, and if R=&(1/t)c^-,
Rt =&d(1/t)t (c^+,r_0), then \[\label{Gbt} \left|I_0(f,g,d,t)- \sum_{k=1}^{L-1}b_k(f,g,d)t^{2k}\right|\leq C(c_{f},c_{g}, d)t^{2L}\log^{2L+2}(1/t).\tag{37}\] Here, \(C(c_{f},c_{g},d)\) is a constant depending only on \(c_{f},c_{g},d\). By changing variables as \(y=x/t\) in 36 , it is easy to see that \(I_0(f, g, d, t)\) is a smooth function of \(t\) in a neighborhood of \(t=0\).

To compute \(b_k(f,g,d)\), we consider any arbitrary fixed \(f,g,d\) such that parts [Ial], [Ial32g] of Assumption 1 are satisfied, and take \(t\to0\). The condition [Rtd] is satisfied for all \(t\) small enough, so 37 implies that \(b_k(f,g,d)=\frac{1}{(2k)!}\partial_t^{2k}I_0(f,g,d,t)\vert_{t=0}\) for all \(k=1,\dots,L-1\). Here, we have used that \(t^{2L}\log^{2L+2}(1/t)=o(t^{2L-1})\), and \(\partial_t^{2k}\) means the partial derivative with respect to the fourth argument, keeping the first three frozen at fixed values. More precisely, if \(f\) and \(g\) depend on \(\lambda\), that \(\lambda\) is held fixed when computing the derivatives. Now that we have this expression for \(b_k\), we are free to substitute any \(f,g,d\) for which 37 is applicable, including \(\lambda\)-dependent \(f,g,d\).

Let \(I_0(t)\) be shorthand for \(I_0(f,g,d,t)\) for a fixed \(f,g,d\). We have \(b_{M/2}(f,g,d)=I_0^{(M)}(0)/(M)!\). To compute this derivative, we modify \(I_0\) to create new functions \(I_1,I_2\), each of which differs from the previous one by \(o(t^M)\). We will then show that \(I_2(t) = \sum_{k=1}^M \tilde{b}_kt^{k}+o(t^M)\) for explicit \(\tilde{b}_k\). This implies \(b_{M/2}=\tilde{b}_M\).
Define the set \[A=\{z\in\mathbb{R}^d\,:\,\|z\|\leq\log(1/t)\sqrt d \}.\]Using Assumption 1, we have that \(g(tz)>0\) for all \(z\in A\) provided \(t\) is small enough, since \(t\log(1/t)\sqrt d\to0\) as \(t\to0\). We can therefore define \[F_g(t,z)=\log g(tz) -\left(t^{-2}f(tz)-\|z\|^2/2\right),\quad z\in A.\]for \(t\) small enough. Thus \(I_0(t)=\log\mathbb{E}[ \exp( F_g(t,Z))\mathbb{1}_A(Z)]\) for a standard Gaussian \(Z\) in \(\mathbb{R}^d\). We now replace \(F_g\) by its Taylor expansion in \(t\) about \(t=0\).

Here and below in this section, the constant factors absorbed in small-\(o\), big-\(\mathcal{O}\), and \(\lesssim\) may depend on any parameter other than \(t\). Note that \(\tilde{F}_g(t,z)\) is precisely the \(M\)th Taylor polynomial of \(F_g(t,z)\) in \(t\) at \(t=0\), and this is how the \(p_k\) are constructed.

Combining Lemmas 11 and 12 gives \(I_0(t)=I_2(t)+o(t^M)\). This implies that the coefficients in front of \(t^k\), \(k=1,\dots,M\), in 37 coincide with the corresponding coefficients in Lemma 13. This concludes the proof of Proposition 2.

To prove Lemma 13, we need an auxiliary result.

See the end of the section for the proof of Lemma 14.

The third line of [G3deriv] gives the first \(M\) terms in the Taylor series expansion of \[\log\mathbb{E}[\exp(tp_1(X)+t^2p_2(X)/2!+\dots+t^Mp_M(X)/M!)],\] which is well-defined. Thus the fourth line of [G3deriv] gives the first \(M\) terms in the formal power series expansion of \(\log\mathbb{E}[\exp(tp_1(Z)+t^2p_2(Z)/2!+\dots+t^Mp_M(Z)/M!)]\), as claimed in Remark 2.

7 Proof of Lemma 3↩︎

We prove an upper bound on \(|I_{\mathrm{out}}(\lambda)|\) and a lower bound on \(I_{\mathrm{in}}(\lambda)\).

Lower bound on \(I_{\mathrm{in}}(\lambda)\).

Assume \(R\sqrt\epsilon\leq r_0\). From [eq:Taylor-f] and [Fk32norms] it follows that f(x)&(1+_u̍̍r_0̍^3f(u)̍_R)x̍̍^2/2
&(1+c_fR)x̍̍^2/2,x̍̍R. We let \(\mu=1+c_{f} R\sqrt\epsilon\). Also, since \(g(0)=1\) and \(\log g\) has bounded first derivative for \(\|x\|\leq R\sqrt\epsilon<r_0\) by [Ial32g] of Assumption 1, we can choose \(R\sqrt\epsilon\) small enough that \(g(x)\geq0.8\) for all \(\|x\|\leq R\sqrt\epsilon\). (Recall that the second inequality of 14 allows us to choose \(R\) small enough.) Thus, \[\label{inlb} I_{\mathrm{in}}(\lambda)\geq 0.8\Big(\frac{\lambda}{2\pi}\Big)^{d/2}\int_{\|x\|\leq R\sqrt\epsilon} e^{-\lambda\mu\left\|x\right\|^2/2}\,\mathrm{d}x=0.8\mu^{-d/2}\mathbb{P}\big(\|Z\|\leq R\sqrt{\mu}\sqrt d\,\big)\geq \frac{1}{2} \mu^{-d/2}.\tag{38}\] To get the last inequality we used that \(R\sqrt{\mu}>2\) and therefore \(\mathbb{P}(\|Z\|\leq R\sqrt{\mu}\sqrt d)\geq \mathbb{P}\big(\|Z\|\leq 2\sqrt d\,\big)\geq 1-\exp(-d/2)\geq 0.63\) when \(d\geq2L\geq2\).

Upper bound on \(|I_{\mathrm{out}}(\lambda)|\).

Let _1&=()^d/2_^d (1+x̍)̍^-d^-px,
T_2&=()^d/2_x̍̍R e^-x̍̍ x,
T_3&=()^d/2_x̍̍R e^-x̍̍^2/2 x. Using 11 , we know that \(|I_{\mathrm{out}}(\lambda)|\leq T_1+T_2+T_3\) for \(p=1/(2L)\). Thus it suffices to upper bound \(T_1\), \(T_2\), \(T_3\). We leave \(p\) unspecified for now to show where the choice \(p=1/(2L)\) comes from. Write \(a=\kappa \lambda d^{-p}\). Using spherical coordinates and assuming \(a>d\) gives \[T_1 \le \frac{\lambda^{d/2}}{2^{(d/2)-1}\Gamma(d/2)}\int_0^\infty (1+r)^{-a}r^{d-1}\,\mathrm{d}r = (\lambda/2)^{d/2}\,\frac{2\Gamma(d)}{\Gamma(d/2)}\cdot \frac{\Gamma(a-d)}{\Gamma(a)}.\] We have \[\frac{\Gamma(a-d)}{\Gamma(a)} = \frac{1}{(a-d)(a-d+1)\cdots(a-1)} \le \Big(\frac{2}{a}\Big)^d,\quad \frac{\Gamma(d)}{\Gamma(d/2)}\le (2d/e)^{d/2}.\] We assumed \(a\ge 2d\) in the first inequality and used the Stirling formula in the second one. Then \[T_1\leq 2\Big(\frac{\lambda}{2}\Big)^{d/2}\Big(\frac{2d}{e}\Big)^{d/2}\Big(\frac{2}{a}\Big)^d =2\Big(\frac{4}{e}\frac{d^{2p+1}}{\kappa^2 \lambda}\Big)^{d/2}.\] We now choose \(p\) so that \(d^{2p+1}/\lambda\) is small whenever the bound on \(I_{\mathrm{in}}\) from Lemma 1 is small. Thus we take \(2p+1=(L+1)/L\), i.e. \(p=1/(2L)\). This gives \[\label{T1} T_1 \leq 2\Big(\frac{4}{e}\frac{d^{(L+1)/L}}{\kappa^2 \lambda}\Big)^{d/2}.\tag{39}\] Returning to the condition \(a\geq 2d\), with \(a=\kappa\lambda d^{-1/(2L)}\), this is satisfied if \(d^{1+1/(2L)}/\lambda< \kappa/2\). But this is implied by the conditions \(R^2\kappa>R\kappa>2\) and \(R^{2(L+1)}d^{L+1}/\lambda^L<1/2\). These latter conditions are satisfied for appropriate choices of \(c^\pm\) in 14 .

Next, using spherical coordinates to compute \(T_2\), we have _2=&_Rd^1/2^e^-d^1/2rr^d-1r d^1/2e^d/2_R^r^d-1e^-drr. Let \(f(r)=(d-1)\log r -\kappa dr\), so that \(r^{d-1}e^{-\kappa dr}=e^{f(r)}\). Note that \(f''(r)<0\) for all \(r>0\), so \(f'(r)\) is decreasing. Therefore, \(f(r)\leq f(R)+f'(R)(r-R)\) for \(r\geq R\). We have \(e^{f(R)}=R^{d-1}e^{-\kappa dR}\) and \(f'(R)=(d-1)/R-\kappa d\). We thus obtain T_2 &d^1/2e^d/2e^f(R)_R^e^f’(R)(r-R)r = e^f(R)
& = e^-d[R-1/2-R]. To get the last inequality we used \(R\kappa d-(d-1)\geq d^{1/2}\), again by assuming \(R\kappa\geq2\).

Finally, using a Gaussian concentration inequality, we have T_3&=^-d_y̍̍Rde^-y̍̍^2/2y^-d(-(R-1)^2)
&(d)e^-(R)^2. To get the first inequality in the second line, we again used \(R\kappa \geq2\). To get the second inequality in the second line, we used \((R\kappa)^2\geq 16\log\frac{1}{\kappa}\), which is satisfied by choosing \(c^-\geq4\) in 14 .

Combining 39 , [T2], and [T3] gives |I_()| &(4e)^d/2 + e^-d[R-1/2-R]+e^-(R)^2
&(4e)^d/2+2e^-d[R-1/2-R]. To get the second line we used that \(\kappa R-1/2-\log R\leq \kappa R-1/2\leq (\kappa R)^2/16\), true for large enough \(\kappa R\). We now finish the proof of 21 using [outub] and 38 .

8 Sampling and expectations for Laplace-type densities↩︎

In statistical applications, the target of study is not a Laplace-type integral but a Laplace-type probability density, \[\pi(x)=\frac{e^{-\lambda f(x)}}{\int_{\mathbb{R}^d} e^{-\lambda f(x')}\mathrm{d}x'}.\] Specifically, one is interested in obtaining the following quantities:

1. i.i.d. samples \(X_1,\dots,X_N\sim\pi\),

2. Expectations \(\mathbb{E}_{X\sim\pi}[g(X)]\).

In Section 8.1 we present our general results for approximating these quantities to arbitrary order \(L\). In Section 8.2 we specialize to the case \(L=1\) and \(L=2\).

8.1 Results for general \(L\)↩︎

We start with the problem of computing expectations for smooth functions \(g\). We can do this using Theorem 1. Indeed, we have \(\mathbb{E}_{X\sim\pi}[g(X)]=\int ge^{-\lambda f}/\int e^{-\lambda f}\), and each integral can be approximated by [main32res].

We have formulated the result uniformly over a function class \(\mathcal{G}\) in order to compare it with our second result below about approximating \(\mathbb{E}_{X\sim\pi}[g(X)]\) for nonsmooth functions \(g\). Note that functions \(g\in\mathcal{G}\) have \(g(0)=1\). This is not restrictive; the theorem also applies for any \(cg\), \(g\in\mathcal{G}\), simply by multiplying [supgG] through by \(c\).

The proof of Theorem 4 is a straightforward application of Theorem 1.

Next, we construct an approximation \(\hat{\pi}_L\) of \(\pi\) which is easy to sample from, and which can be used to approximate \(\mathbb{E}_{X\sim\pi}[g(X)]\) for nonsmooth \(g\). The following theorem is our second main result.

The maps in the composition increase in complexity as one goes outward: \(Q\) is linear, and \(T_{L-m}\) is a polynomial of degree \(2m\), \(m=1,\dots, L\).

Theorem 6 indeed gives a tractable algorithm for approximately sampling from \(\pi\): simply draw \(Z_i\sim\mathcal{N}(0, \lambda^{-1}I_d)\) i.i.d. and return \(x_L(Z_i)\). To do this, we push \(Z_i\) through the sequence of \(L+1\) maps \(Q, T_{L-1},T_{L-2},\dots,T_0\). In the case \(L=1\), we explicitly construct \(x_1\) in Section 8.2.

See also further discussion of the uses for Theorems 4 and 6 in Section 9.4.

8.2 Special cases: \(L=1\) and \(L=2\)↩︎

We derive the approximation from Theorem 4 in the cases \(L=1\), \(L=2\), and the approximation from Theorem 6 in the case \(L=1\).

Theorem 4 with \(L=1\) gives that for all \(g\in\mathcal{G}_1(c_{g},\kappa_g)\), we have \[\mathbb{E}_{X\sim\pi}[g(X)]=1+\mathcal{O}\left(R^4\frac{d^2}{\lambda}\right),\]where \(\mathcal{O}\) suppresses dependence on \(c_{g}\) and \(c_{f}\). Here, recall that we assume \(g(0)=1\). To work out the case \(L=2\), we compute \(b_1(f,g)-b_1(f,1)\). Recall the formula for \(b_1(f,g)\) from [b1fg]. When we subtract \(b_1(f,1)\), all terms involving only \(f\) will cancel. The remaining expression is \(b_1(f,g)-b_1(f,1)=-\frac{1}{2}\nabla\Delta f(0)^\top\nabla g(0)+\frac{1}{2}\Delta g(0)\). We conclude that for all \(g\in\mathcal{G}_2(c_{g},\kappa_g)\), we have \[\mathbb{E}_{X\sim\pi}[g(X)]=\exp\left(-\frac{1}{2\lambda}\nabla\Delta f(0)^\top\nabla g(0)+\frac{1}{2\lambda}\Delta g(0)\right)+\mathcal{O}\left(R^6\frac{d^3}{\lambda^2}\right).\]

Next, we compute the map \(x_1=T_0\circ Q\) from Theorem 6. The map \(T_0\) is \(T_0=X\) from Lemma 4. Below this lemma, we showed that when \(L=1\) we have \[\label{T0-1}T_0(t)=t-\frac{1}{6}\nabla^3f(0)[t^{\otimes2}].\tag{48}\] To determine \(Q\), we need to derive \(E_1(t)\) and write it in the form 29 . The function \(E_1\) is given by adding the nonnegligible part of \(\frac{1}{\lambda}\log\det T_0'(t)\) to \(\|t\|^2/2\). We have \[\frac{1}{\lambda}\log\det T_0'(t)=\frac{1}{\lambda}\mathrm{tr}\log\left(I_d-\tfrac1{3}\nabla^3f(0)[t]\right)=-\tfrac1{3\lambda}\mathrm{tr}(\nabla^3f(0)[t])+\epsilon\mathcal{O}(\|t\|^2),\] as in Lemma 6 but explicitly computing that \(dF_{1}^{\epsilon}[t]=-\tfrac13\mathrm{tr}(\nabla^3f(0)[t])\). Here, \(\nabla^3f(0)[t]\) is the matrix with \((i,j)\)th entry given by \(\nabla^3f(0)[t, e_i, e_j]\). Throwing out \(\epsilon\mathcal{O}(\|t\|^2)\), we conclude \[\label{E1-1}E_1(t)=\frac{1}{2}\|t\|^2 + \frac{1}{3\lambda}\mathrm{tr}(\nabla^3f(0)[t])=\frac{1}{2}\|t\|^2+\frac{1}{3\lambda}\nabla\Delta f(0)^\top t.\tag{49}\] Comparing with  29 , we see that \(a=\frac{1}{3\lambda}\nabla\Delta f(0)\) and \(B=I_d\). Thus Theorem 6 gives \[\label{S-1} Q(s) = B^{-1/2}(s-B^{-1/2}a)=s-\frac{1}{3\lambda}\nabla\Delta f(0).\tag{50}\] Finally, \(x_1=T_0\circ Q\). Thus Theorem 6 with \(L=1\) gives that _X~[g(X)] &= )] + (R^4),
S&~N(-f(0), ^-1I_d) for all \(\|g\|_\infty\leq1\), where \(\mathcal{O}\) suppresses dependence on \(c_{f}\) only. Here, we have used that \(\hat{\pi}_1=(T_0\circ Q)_{\#}\mathcal{N}(0,\lambda^{-1}I_d)\), which is the pushforward under \(T_0\) of \(Q_{\#}\mathcal{N}(0,\lambda^{-1}I_d)=\mathcal{N}(-\frac{1}{3\lambda}\nabla\Delta f(0), \lambda^{-1}I_d)\).

The righthand expectation in [x1-g] typically cannot be evaluated in closed-form, but can easily be approximated by Monte Carlo using samples \(S_i\sim \mathcal{N}(-\frac{1}{3\lambda}\nabla\Delta f(0), \lambda^{-1}I_d)\).

9 Comparison to the literature and examples↩︎

In Section 9.1, we compare our integral expansion result to that of [4]. Then in Sections 9.2 and 9.3, we derive the integral expansion for a few examples. In Section 9.4, we discuss our results involving Laplace-type densities from Section 8, and compare them to related work in the literature.

Throughout the section, we assume \(d\geq\log\lambda\) to simplify formulas.

9.1 Comparison to integral expansions in the literature↩︎

In the below informal discussion of our results, we neglect any dependence on \(r_0,\kappa\). Since also \(d\geq\log\lambda\), a constant \(R\) can be used to satisfy 14 . Theorem 1 then shows that \[\label{curr} I(\lambda)=\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d}g(x)e^{-\lambda f(x)}\mathrm{d}x = \exp\left(\sum_{k=1}^{L-1}b_k\lambda^{-k} + \mathcal{O}(d^{L+1}/\lambda^L)\right),\quad |b_k|\lesssim d^{k+1},\tag{51}\] provided \(\max_{1\leq k\leq 2L}\|\nabla^k(\log g)(x)\|_\mathrm{op}< c_{g}\) and \(\max_{3\leq k\leq 2L+2}\|\nabla^kf(x)\|_\mathrm{op}< c_{f}\), uniformly over \(x\) in a small neighborhood of 0. As discussed in the introduction, in [4] it is shown that \[\label{prev} I(\lambda)=\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d}g(x)e^{-\lambda f(x)}\mathrm{d}x =\sum_{k=1}^{L-1}a_k\lambda^{-k} + \mathcal{O}(d^{2L}/\lambda^L),\quad |a_k|\lesssim d^{2k},\tag{52}\] provided \(\|\nabla^kg(x)\|_\mathrm{op}\lesssim d^{\lceil k/2\rceil}\), \(k=1,\dots,2L\) and \(\|\nabla^kf(x)\|_\mathrm{op}\lesssim d^{\lceil k/2\rceil-2}\), \(k=3,\dots,2L+2\), uniformly over \(x\) in a small neighborhood of 0. The suppressed constants in the big-\(\mathcal{O}\) and in the bound on \(|a_k|\) in 52 depend on the suppressed constants in the derivative bounds.

Both results also require a few other minor assumptions, but we have highlighted the most important ones for this discussion. We make a few comments on the difference between the two works.

1. For 51 , we have assumed the derivative operator norms are bounded independently of \(d\). This was not required in [4]. However, if the assumption does hold, then our result is strictly stronger than that of [4]. Specifically, suppose the operator norms are bounded and \(d^2\ll\lambda\). Then 52 can be derived from 51 .

2. As can be seen from [4], \(d^2\ll \lambda\) is necessary for the expansion of \(I(\lambda)\) even when the derivative operator norms are bounded. Indeed, consider the example in Section 9.2, also studied in Example 2.19 in [4]. We see that the derivative operator norms are indeed bounded, yet [4] proves that the expansion of \(I(\lambda)\) is valid only if \(d^2\ll \lambda\). This shows that our improved dimension dependence cannot simply be attributed to the stricter requirement we have imposed on the operator norms. Rather, it is due to the intrinsic difference between expanding \(I(\lambda)\) and expanding \(\log I(\lambda)\), as described in the introduction.

3. In the expansion of \(I(\lambda)\) in [4], the derivative operator norms beyond the fourth order were allowed to grow with \(d\) due to slack in the bound. A similar phenomenon may hold for the expansion of \(\log I(\lambda)\). Namely, there may be some slack which would permit derivative norm growth with \(d\). We leave this investigation to future work.

9.2 Example: quartic exponent↩︎

Let \(f(x)=\|x\|^2/2 + \|x\|^4/24\) and \(g(x)\equiv1\). Consider [eq:Taylor-f] and [Fk32norms]. For all \(L\geq2\) we have \(\mathcal{R}_{2L+2}(x)\equiv0\), and \(\nabla^3f(0)[x^{\otimes3}]=0\), \(\nabla^4f(0)[x^{\otimes4}]=\|x\|^4\), \(\nabla^kf(0)[x^{\otimes k}]=0\) for all \(k\geq5\). Thus [eq:Taylor-f] and [Fk32norms] are satisfied for any \(r_0\). Also, \(\log g\equiv 0\) trivially satisfies [eq:Taylor-g], [Gk32norms] for any \(r_0\). Furthermore, 11 holds with \(\kappa=1\), since \(|g(x)|e^{-\lambda f(x)}\leq e^{-\lambda\|x\|^2/2}\).

Theorem 1 therefore applies. Since we assume \(d\geq\log\lambda\) throughout the section, we can take \(R\) to be an absolute constant in 14 . We conclude that {()^d/2_^de^-(+)x} =_k=1^L-1 b_k^-k+(d^L+1/^L), for all \(d^{L+1}/\lambda^L\) small enough. Here, the constant suppressed by the big-\(\mathcal{O}\) depends on \(L\) only. Next, let us compute \(b_1\), given by [b1fg]. Plugging in \(f(x)=\|x\|^2/2 + \|x\|^4/24\) and \(g(x)\equiv1\) gives b_1 &=-^2f(0)=-_i,j=1^d_i^2_j^2(x̍̍^4)
&= -_ij_i^2_j^2(2x_i^2x_j^2)-_i=1^d_i^4(x_i^4)
&=-(d^2-d) -18d = -d^2 -d. Using this in [quartL] with \(L=2\), we conclude that \[\label{I1} I(\lambda)=\left(\frac{\lambda}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d}e^{-\lambda\left(\frac{\|x\|^2}{2}+\frac{\|x\|^4}{24}\right)}\mathrm{d}x = \exp\left(-\frac{d^2}{24\lambda} -\frac{d}{12\lambda} + \mathcal{O}(d^3/\lambda^2)\right),\tag{53}\] where the constant suppressed by the big-\(\mathcal{O}\) is absolute.

This example is also studied in [4], where it is shown that \(I(\lambda)= 1-(d^2/24+d/12)/\lambda+ \mathcal{O}(d^4/\lambda^2)\), but only if \(d^2/\lambda\ll1\). Our result improves on [4] in that it both tightens the remainder bound and broadens the range of applicability of the expansion into higher \(d\) regimes.

9.3 Example: logistic regression↩︎

Next, we consider an idealized statistical set-up, in which \(\int_{\mathbb{R}^d}e^{-\lambda f(x)}\mathrm{d}x\) is the normalizing constant of the posterior in a logistic regression model. See Section [ssec:statistics95motivation] for more details on how this quantity arises in statistics. We call the large parameter \(\lambda=n\) the sample size. Let \(X_1,\dots, X_n\in\mathbb{R}^d\) be the feature column vectors, and \(x^*\in\mathbb{R}^d\) be the ground truth we would like to estimate. We assume (1) \(X_1,\dots,X_n\) span \(\mathbb{R}^d\), (2) the lowest eigenvalue of \(\frac{1}{n}\sum_{i=1}^nX_iX_i^\top\) is bounded below away from zero by some \(c_X>0\), and (3) \(\max_{i=1,\dots,n}\|X_i\|\max(1, \|x^*\|)<C_X\). Next, let \(\psi:\mathbb{R}\to\mathbb{R}\) be a smooth strictly convex function, with \(\|\psi^{(k)}\|_\infty=\sup_{t\in\mathbb{R}}|\psi^{(k)}(t)|<\infty\) for all \(k=2,3,4,\dots\). Consider the function \[\ell(x) = \frac{1}{n}\sum_{i=1}^n\left[\psi(X_i^\top x)-\psi'(X_i^\top x^*)X_i^\top x\right]\] which has a unique global minimizer at \(x=x^*\). When \(\psi(t)=\log(1+e^{t})\), the function \(\ell\) is the negative, normalized population log likelihood for the logistic regression model. (The fact that \(\ell\) is the population rather than the sample log likelihood makes this an idealized set-up.)

Let \[H:=\nabla^2\ell(x^*)=\frac{1}{n}\sum_{i=1}^n\psi''(X_i^\top x^*)X_iX_i^\top.\] The matrix \(H\) is strictly positive definite. Indeed, let \(c_\psi=\inf_{|t|<C_X}\psi''(t)\), which is positive since \(\psi''(t)>0\) for all \(t\in\mathbb{R}^d\) and \(\psi\) is smooth. Then \(H\succeq c_\psi\,\frac{1}{n}\sum_{i=1}^nX_iX_i^\top\succeq c_\psi c_XI_d\). By a change of variables, we have \[\label{ell-f}\int_{\mathbb{R}^d}e^{-n\ell(x)}\mathrm{d}x =\frac{e^{-n\ell(x^*)}}{\sqrt{\det H}}\int_{\mathbb{R}^d}e^{-nf(x)}\mathrm{d}x,\qquad f(x) = \ell(H^{-1/2}x+x^*)-\ell(x^*).\tag{54}\] We verify the conditions of Theorem 1 for this \(f\). We have \(f(0)=0\), \(\nabla f(0)=0\) and \(\nabla^2f(0)=I_d\). Thus the Taylor expansion of \(f\) indeed takes the form in [eq:Taylor-f]. Furthermore, using the above lower bound on \(H\) and the definition of \(\ell\), we have for \(k\ge 3\) _x^d̍^kf(x)̍_&(c_c_X)^-k/2_x^d̍^k(x)̍_
&= (c_c_X)^-k/2_x^d̍1n_i=1^n^(k)(X_i^x)X_i^k̍_
&(c_c_X)^-k/2̍^(k)̍__i=1,…,nX̍_i̍^k
&(c_c_X)^-k/2̍^(k)̍_C_X^k. Thus [Fk32norms] is satisfied for any \(r_0\). The condition on \(g\equiv1\) in Assumption 1, part [Ial32g] is trivially satisfied. Finally, to verify 11 , we prove the sufficient condition 12 holds. Even more simply, it suffices to show \(f(x)\geq\kappa\min(\|x\|^2/2, \sqrt{\epsilon}\|x\|)\) for all \(x\in\mathbb{R}^d\). Let \(c=\max_{u\in\mathbb{R}^d}\|\nabla^3f(x)\|_\mathrm{op}\), which we know is bounded, by the above calculation. Then for all \(\|x\|\leq\sqrt\epsilon\), a Taylor expansion gives \(f(x)\geq (1-c\sqrt\epsilon/3)\|x\|^2/2 \geq\frac{1}{2}(\|x\|^2/2)\), if \(\epsilon\) is small enough. For all \(\|x\|\geq\sqrt\epsilon\), convexity of \(f\) implies that \(f(x)\geq \frac{\|x\|}{\|y\|}f(y)\), where \(\|y\|=\sqrt\epsilon\) and \(y\) lies on the line segment between \(0\) and \(x\). But now we use that \(f(y)\geq\|y\|^2/4\) since \(\|y\|=\sqrt\epsilon\), so \(f(x)\geq\|x\|\|y\|/4=\|x\|\sqrt\epsilon/4\). Therefore, \(f(x)\geq\kappa\min(\|x\|^2/2, \sqrt{\epsilon}\|x\|)\) is satisfied on \(\mathbb{R}^d\), with \(\kappa=1/4\).

Theorem 1 therefore applies to \(f\) from 54 . As in Section 9.2, we can satisfy the lower bound in 14 by taking \(R=c\) for some constant \(c\) depending on \(C_X\), \(c_x\), \(\psi\), and \(L\). Writing \[\int_{\mathbb{R}^d}e^{-n\ell(x)}\mathrm{d}x= \frac{e^{-n\ell(x^*)}}{\sqrt{\det H}}\left(\frac{2\pi}{n}\right)^{d/2}\times \left\{\left(\frac{n}{2\pi}\right)^{d/2}\int_{\mathbb{R}^d}e^{-nf(x)}\mathrm{d}x\right\},\] we conclude that _^de^-n(x)x = &-n(x^*)-H -d2
&++( ), whenever \(d^3/n^2\) is sufficiently small. Here, the constant suppressed by the big-\(\mathcal{O}\) depends on \(C_X,c_x,\psi\). The term \(b_1\) can be computed using the general formula [b1fg], and similar calculations as in [4], [5]. Note that the righthand side of the first line in [logreg-L] is the Bayesian information criterion (BIC) [54]. Thus \(b_1/n\) can be considered a higher-order correction to BIC, with the higher-order remainder explicitly controlled.

9.4 Comparison to Laplace-type density approximations in the literature↩︎

Here, we discuss the problem of approximating Laplace-type densities. We first compare the two approximation strategies from Theorems 4 and 6, showing that the former uses fewer \(f\) derivatives at the cost of requiring smoothness of \(g\), while the latter offers sampling flexibility. We then compare the combination of these two methods to other approaches in the literature. Our discussion in this section is focused on the statistical context.

Recall from Section 8 that Theorem 4 gives the approximation \[\label{Eg-expand} \mathbb{E}_{X\sim \pi}[g(X)]= \exp\left(\sum_{k=1}^{L-1}[b_k(f,g)-b_k(f,1)]\lambda^{-k}\right) + \mathcal{O}(d^{L+1}/\lambda^L),\tag{55}\] while Theorem 6 gives \[\label{Eg-expect} \mathbb{E}_{X\sim \pi}[g(X)]= \mathbb{E}_{X\sim \hat{\pi}_L}[g(X)]+ \mathcal{O}(d^{L+1}/\lambda^L),\qquad \hat{\pi}_L:=(x_L)_{\#}\mathcal{N}(0, \lambda^{-1}I_d).\tag{56}\] Here as above, we assume \(d\geq\log\lambda\) and neglect dependence on \(r_0,\kappa\). This allows us to take \(R\) to be constant, which is why \(R\) does not appear in the big-\(\mathcal{O}\)’s above.

When 55 is applicable, it is a more powerful method for computing expectations than 56 . The reason for this is two-fold. First, 56 cannot be implemented exactly, and requires a further sampling step, as seen in 8 . This incurs additional computational cost and loss of accuracy. On top of this, using \(\hat{\pi}_L\) requires computing higher-order \(f\) derivatives than are involved in \(\sum_{k=1}^{L-1}[b_k(f,g)-b_k(f,1)]/\lambda^k\). Indeed, Remarks 8 and 5 show that \(\hat{\pi}_L\) uses \(2L+1\) derivatives of \(f\) derivatives, while \(\sum_{k=1}^{L-1}[b_k(f,g)-b_k(f,1)]/\lambda^k\) uses only \(2L-1\). The reason for this gap is that the two approximations achieve the same accuracy, but \(\hat{\pi}_L\) does not exploit any smoothness of \(g\); it must compensate by extracting more information from \(f\). In practice, this difference in derivative count can be significant for computational cost. Note that the closed-form approximation 55 does require derivatives of \(g\), while the pushforward approximation 56 does not. However, in Bayesian statistics \(g\) is typically a simple function (e.g. linear or quadratic, corresponding to the mean or covariance of \(\pi\)), making its derivatives cheap to compute. In contrast, \(f\) encodes the data likelihood and involves a sum over \(\lambda=n\) terms, so each additional derivative of \(f\) carries substantial computational cost.

The number of derivatives of \(f\) and \(g\) used in the approximations 55 and 56 are summarized in Table 1 for the cases \(L=1\) and \(L=2\).

For convenience, we remind the reader of the explicit formulas for three of the approximations referenced in the table: \[\begin{align} {2} &\mathbb{E}_{X\sim\pi}[g(X)]\approx 1,\qquad &g\in \mathcal{G}_1(c_{g},\kappa_g),\tag{57}\\ &\mathbb{E}_{X\sim\pi}[g(X)] \approx \mathbb{E}\left[g\left(S-\tfrac16\nabla^3f(0)[S,S,\cdot]\right)\right] ,\qquad &\|g\|_\infty\leq1,\tag{58}\\ &S\sim\mathcal{N}\left(-\tfrac1{3\lambda}\nabla\Delta f(0), \lambda^{-1}I_d\right),&\notag\\ &\mathbb{E}_{X\sim\pi}[g(X)]\approx\exp\left(-\tfrac1{2\lambda}\nabla\Delta f(0)^\top\nabla g(0)+\tfrac1{2\lambda}\Delta g(0)\right),\qquad &g\in \mathcal{G}_2(c_{g},\kappa_g).\tag{59} \end{align}\] These were derived in Section 8.2. We have omitted the formula for the approximation via \(\mathbb{E}_{X\sim\pi_2}[g(X)]\) (i.e. the fourth column), which requires 5 derivatives of \(f\).

Although the density approximations \(\hat{\pi}_L\) to \(\pi\) are more derivative-intensive, they can be used to approximate expectations of nonsmooth functions \(g\), which the closed-form method 55 cannot do. Furthermore, approximately sampling from \(\pi\) is itself valuable, beyond just computing expectations. For example, in Bayesian statistics, samples can be used to construct approximate credible intervals to quantify uncertainty in the target parameter of inference. The algorithm to sample from \(\hat{\pi}_1\) is especially simple:

1. Draw \(S_i\sim\mathcal{N}(-\frac{1}{3\lambda}\nabla\Delta f(0),\lambda^{-1}I_d)\) i.i.d.

2. Return \(X_i=S_i-\frac{1}{6}\nabla^3f(0)[S_i,S_i,\cdot]\).

The combined ability to 1) easily generate approximate samples from \(\pi\) (via the above algorithm), 2) approximate expectations of nonsmooth \(g\) (via the Monte Carlo estimate), and 3) accurately and cheaply approximate expectations for smooth \(g\) (via 57 and 59 ) is extremely powerful compared to the state of the art.

Two noteworthy alternative methods in the literature for high-accuracy sampling and computing expectations are [37] and [40]. The two works construct approximations \(\hat{P}_{\mathrm{SKS}}\) and \(\hat{\gamma}_S\) to \(\pi\), respectively. They are similar to our \(\hat{\pi}_1\), involving only the second and third derivatives of \(f\). (The reason our \(\hat{\pi}_1\) only involves the third derivative is that we have assumed \(\nabla^2f(0)=I_d\). For a generic \(\nabla^2f(0)\), the second derivative will also arise.)

We argue that our combined approach takes the best elements of each of the approximation methods in the above works. For sampling accuracy, the relevant metric is TV distance. In  [37], the authors show only that \(\mathrm{TV}(\pi,\hat{P}_{\mathrm{SKS}})\lesssim d^3/\lambda\), whereas we show the tighter dimension dependence \(\mathrm{TV}(\pi,\hat{\pi}_1)\lesssim d^2/\lambda\). Furthermore, while \(\hat{P}_{\mathrm{SKS}}\) is easy to sample from, it cannot be integrated against in closed-form. As a result, Monte Carlo sampling is always needed for the purpose of computing expectations, even of smooth functions. As discussed above, this incurs extra computational cost as well as an additional source of error. But even if it were possible to compute \(\mathbb{E}_{X\sim \hat{P}_{\mathrm{SKS}}}[g(X)]\) exactly, the accuracy of the approximation remains \(d^3/\lambda\) at best. In contrast, by exploiting the smoothness of \(g\), our estimate 59 achieves the much higher accuracy \(d^3/\lambda^2\) while using the same number of derivatives of \(f\). This is a significant improvement by a factor of \(1/\lambda\).

The approximation \(\hat{\gamma}_S\) of [40] is a signed measure, not a true probability density. It has the advantage that expectations of polynomials against \(\hat{\gamma}_S\) can be computed in closed-form, unlike \(\hat{P}_{\mathrm{SKS}}\). However, our 56 gives a closed-form approximation for expectations of \(g\) in the even broader class of smooth functions, not just polynomials. Also, the fact that \(\hat{\gamma}_S\) is not a true probability density has disadvantages; for example, “sampling" from this signed measure is not well-defined. Our \(\hat{\pi}_1\) does not have this issue: it is a true probability density and can be easily sampled from.

Another significant advantage of our results, compared to those of [40] and [37], is that we give a method to approximate \(\pi\) (and expectations under \(\pi\)) to arbitrary order of accuracy. In contrast, the other two works focus only on a fixed order of approximation, and it is unclear whether their constructions (or proof techniques) can be extended to higher orders of accuracy.

Finally, it is natural to compare 55 to the analogous result from [4], which can be used to expand the Laplace integral in the numerator and denominator as follows: \[\mathbb{E}_{X\sim\pi}[g(X)]=\frac{\sum_{k=0}^{L-1}c_k(f,g)\lambda^{-k}+\mathcal{O}((d^2/\lambda)^L)}{\sum_{k=0}^{L-1}c_k(f,1)\lambda^{-k}+\mathcal{O}((d^2/\lambda)^L)}.\] But the drawback of this result, as already discussed in Section 9.1, is that it does not allow \(d\) to be larger than \(\sqrt\lambda\).

In summary, our approach combines the best features of these methods — closed-form expectations for smooth \(g\), easy sampling from a true density — and adds two more features: arbitrary-order accuracy, and validity up to the concentration threshold.

10 Structure-preserving tensor operations↩︎

Recall from Section 2 the notation of a tensor \(A_{k\to j}\): a multilinear mapping from \(k\) vectors in \(\mathbb{R}^d\) to either a scalar if \(j=0\), a vector in \(\mathbb{R}^d\) if \(j=1\), and a \(d\times d\) matrix if \(j=2\). Recall also the concept of a base tensor \(G_{k\to j}\) and a composite tensor \(F_{k\to j}^{\epsilon}\).

In the next lemma, we list some operations which preserve the structure of a composite tensor. Recall that all the composite tensors \(F_{k\to j}^{\epsilon}\) and all the base tensors \(G_{k\to j}\) they are composed of are symmetric.

11 Change of variable proofs↩︎

For the proof of Lemma 4, we state and prove a main auxiliary lemma. Recall that \(F_{k}^{\epsilon}\) is shorthand for \(F_{k\to 0}^{\epsilon}\).

For the proof of Lemma 8, we state and prove a main auxiliary lemma.

References↩︎