Sharp regret-Hellinger bounds for Gaussian empirical Bayes via polynomial approximation

Sharp regret--Hellinger bounds for Gaussian empirical Bayes via polynomial approximation

\ifpdf \DeclareGraphicsExtensions{.pgf} \graphicspath{{figures/}{plots/}} \fi

Introduction

Empirical Bayes (EB) is a powerful framework for large-scale statistical inference that applies Bayesian methods with a data-driven prior, thereby adapting to the latent structure in the data robbins1956empirical,zhang2003compound. For an account of the remarkable successes of EB in both theory and practice, we refer to the monograph and survey of Efron efron2012large,efron2024empirical. A key benchmark for an EB procedure is its regret: the excess risk relative to the Bayesian oracle that knows the true prior, or, in the compound setting, the empirical distribution of the parameters.

A standard setting for EB is the normal means model. Consider $\theta \sim G$ and $Y=\theta+Z$, where $Z\sim N(0,1)$ is independent of $\theta$. Tweedie's formula expresses the Bayes estimator in terms of the score of the marginal distribution of $Y$, namely, the Gaussian mixture $f_G := G * N(0,1)$:

This estimator minimizes the mean-squared error (MSE) among all estimators of $\theta$: $\min_{\widehat\theta(\cdot)} \mathbb{E}_G[(\theta-\widehat\theta(Y))^2] = \mathbb{E}_G[(\theta-m_G(Y))^2]$.

In the context of empirical Bayes, $m_G$ in \prettyref{eq:tweedie} is referred to as oracle Bayes, which requires the knowledge of the true prior $G$. When the Bayes rule $m_H$ with a misspecified prior $H$ is used, the regret, namely, the excess MSE over the oracle risk, is given by

the second equality applying the Tweedie's formula \prettyref{eq:tweedie}. This quantity is also known as relative Fisher information in information theory and optimal transport villani.topics,raginsky2013concentration and score matching loss in generative modeling hyvarinen2005estimation.

In practice---and especially in the $g$-modeling approach efron2014two---$H$ is a prior estimated from a training sample and $\mathsf{Regret}(H\|G)$ is the incurred regret when the learned Bayes rule $m_H$ is applied to a fresh test data point. Widely used approaches for learning the prior include maximum likelihood and minimum distance, both of which aim at approximating the true data-generating distribution $f_G$. As such, a central problem in EB theory is to control the regret $\mathsf{Regret}(H\|G)$ using statistical distance between $f_G$ and $f_H$, most commonly the squared Hellinger distance:

A foundational result in nonparametric EB of this type is due to Jiang and Zhang JiangZhang2009. A substantial body of subsequent work in statistics and econometrics builds on this result; see, for instance, jiang2020general,saha2020nonparametric,soloff2025multivariate,SW22,adusumilli2025empirical,ghosh2025stein,chen2026empirical,ignatiadis2026compound,kim2026empirical. In addition, this result also finds applications beyond the classical realm of EB; for example, it plays a crucial role in constructing optimal kernel-based estimators of score functions for diffusion models WWY24.

A key ingredient of Jiang--Zhang's result JiangZhang2009 controls the following regularized version of the regret \prettyref{eq:regret}: \[ \mathsf{Regret}_{\rho}(H \Vert G) := \int_\R \pr{ \frac{f'_G}{f_G \vee \rho} - \frac{f'_H}{ f_H \vee \rho} }^2 f_G dy. \addtocounter{equation}{1}\tag{\theequation}\label{eq:regularized_regret} \] This quantity is related to the regret of the regularized Bayes estimator

where the regularization parameter $\rho$ ensures the denominator is not too small. Using a delicate recursive argument, Jiang and Zhang showed that \prettyref{eq:regularized_regret} can be bounded by the Hellinger distance $\epsilon := d_{\rm H}(f_G,f_H)$ as follows:

This quantity can be further related to the regret of the regularized estimator \prettyref{eq:bayes-rho} by quantifying the regularization error on the true Bayes rule, under suitable tail assumptions on the prior $G$. In practice, $\rho$ is typically chosen as a small polynomial of $\epsilon$, leading to a regret bound of the order $O(\epsilon^2 (\log \frac{1}{\epsilon})^3)$.

The key difficulty of proving a regret--Hellinger inequality like \prettyref{eq:JZ09} is that regret compares the scores (derivative of log density) while Hellinger only compares the densities themselves. To bridge this gap, Jiang and Zhang introduced a sequence $\Delta_k$ defined in terms of the $k$th derivative of $f_G-f_H$, where $\Delta_0$ and $\Delta_1$ are approximately the Hellinger and the regularized regret, respectively. They then proved a remarkable recursive inequality $\Delta_k^2 \leq \Delta_{k-1}(\Delta_{k+1} + C \sqrt{\log(1/\rho)} \Delta_k)$ and bounded $\Delta_k$ separately using Fourier analysis, concluding the desired \prettyref{eq:JZ09} by backward induction.

Since its introduction, the Jiang--Zhang method has remained essentially the only route for bounding the regret of $g$-modeling estimators,\footnote{One exception is JPW21 for the Poisson EB, where the regret can be bounded optimally by the Hellinger thanks to the discreteness of the Poisson model and the light tail of the prior. Nevertheless, this method fails for the moment class and an argument \`a la Jiang--Zhang is still needed SW22.} with extensions to multiple dimensions saha2020nonparametric, higher-order estimands ghosh2025stein (see also the earlier version of chen2026empirical), heteroskedastic variances jiang2020general,soloff2025multivariate,chen2026empirical, and the Poisson model SW22. At the same time, it is long suspected that this argument may not be tight,\footnote{JiangZhang2009 remarked that “Still, the cubic power of the logarithmic factors [...] is crude.” and that “we do not believe that mathematical induction is sharp in the argument with higher and higher order of differentiation in the proof of Lemma 1.” } partially contributing to a $\operatorname{polylog}(n)$ gap between the best known regret upper bound and information-theoretic limits polyanskiy2021sharp,kang2026function. Moreover, it is also not obvious whether the regularization in ((ref)) is merely a proof technique or plays a more fundamental role for empirical Bayes.

\paragraph{Main contributions.} We present a new technique for bounding regret in terms of Hellinger distance without regularization. Compared to the recursive argument of Jiang and Zhang, the proposed method is conceptually simpler and yields sharper, sometimes optimal, regret bounds. The major contributions of this paper are as follows.

\paragraph{Proof ideas.} As discussed earlier, controlling regret in terms of Hellinger distance entails bounding a derivative in terms of the original function. Such a bound is clearly impossible in general, but it becomes possible for certain classes of regular functions, most notably polynomials. This crucial insight is gained through interactions with GPT-5.4 using OpenAI Codex. It motivates the following program:

Putting these ingredients together and optimizing the degree yield a regret--Hellinger inequality.

There are multiple ways to execute this program. To fix ideas, let us focus on the class of compactly supported priors. One may choose $g = \frac{f_G-f_H}{f_G+f_H}$ and $w=(f_G+f_H)/2=f_{(G+H)/2} * \varphi$, in order to capitalize on the classical Bernstein inequality for the standard Gaussian density $\varphi$ (an immediate consequence of properties of Hermite polynomials; see \prettyref{rmk:bernstein-perturb}):

and hence (by linearity) for any Gaussian mixture. In addition, the best degree-$k$ polynomial approximation error of this choice of $g$ is at most $\mathop{\rm exp}(-\Omega(\sqrt{k}))$. Choosing $k =\Theta((\log \frac{1}{\epsilon})^2)$ leads to $\mathsf{Regret}(H\|G) = O(\epsilon^2 (\log \frac{1}{\epsilon})^2)$, which already improves over the state of the art of JiangZhang2009 and does so without regularization.

To obtain the sharp inequality, we instead choose \[ g = \frac{f_G-f_H}{\varphi}, \quad w= \frac{\varphi^2}{(f_G+f_H)/2}. \] This $g$ admits a substantially better polynomial approximation error of $\mathop{\rm exp}(-\Omega(k \log k))$. In addition, we prove a Bernstein-type inequality \prettyref{eq:bernstein-general} for the non-standard weight $w$ with $a_k = O(\sqrt{k})$. Setting $k =\Theta(\log \frac{1}{\epsilon}/ \log\log \frac{1}{\epsilon})$ yields the desired \prettyref{eq:regret-compact-intro}. We note that there is a sizable literature on Bernstein-type inequalities for weighted norms (see, e.g., LevinLubinsky2001,lubinsky2007survey); however these results require additional assumptions that appear difficult to verify in our setting. Instead, we present self-contained proofs using a linear algebraic argument, as the inequality \prettyref{eq:bernstein-general} boils down to bounding the operator norm of the matrix representing the differentiation operator on the space of polynomials.

Main results

New Regret--Hellinger inequality

The following result shows that the regret upper bound in \prettyref{thm:ub-compact-intro} is sharp up to constants:

We now extend \prettyref{thm:ub-compact-intro} to priors with exponential tails. For $\alpha, \sigma>0$, define

In the special case of $\alpha=2$, this is the class of subgaussian distributions with proxy variance proportional to $\sigma^2$; for $\alpha = 1$, it is the subexponential class.

For heavy-tailed priors---for instance, those satisfying only a bounded $p$th-moment condition---we show in \prettyref{app:moment} that regret can be polynomially larger than the squared Hellinger distance. Hence the regularization of JiangZhang2009 is not merely a proof device but is in fact necessary in order to obtain a regret nearly linear in $d_{\rm H}^2$.

Statistical applications

Next we discuss implications of these regret--Hellinger inequalities for empirical Bayes. To simplify presentation and compare with existing results, we consider priors that are either compactly supported or subgaussian.

One of the most popular methods for learning a prior is the nonparametric maximum likelihood estimator (NPMLE) KW56:

where $y_1,\ldots,y_n$ are the given data points and the maximization is over all probability measures on $\mathbb{R}$. It is known that \prettyref{eq:NPMLE} admits a unique and discrete maximizer, which has at most $n$ atoms in the worst case and $O(\log n)$ atoms in the typical case when $y_i$'s are drawn independently from a subgaussian distribution simar1976maximum,lindsay1983geometry1,PW20-npmle.

Suppose the NPMLE is computed based on a training set $y_1,\ldots,y_n$ and then applied to a fresh test data point. The resulting regret (excess over the Bayes squared error) is bounded by the following result, which improves over the state of the art $O((\log n)^5/n)$ in JiangZhang2009 for regularized NPMLE.

The proof of \prettyref{cor:NPMLE} is a direct consequence of \prettyref{thm:ub-exp}, once $\widehat G_{\rm NPMLE}$ is shown to be $O(1)$-subgaussian with high probability. This is done in \prettyref{app:NPMLE}.

Several further remarks are in order.

Compactly supported priors

We prove a more precise version of \prettyref{thm:ub-compact-intro} that is applicable to growing $M$.

A first reduction

We first bound $\mathsf{Regret}$ by quantities $\delta, \Delta$, defined below.

Note that $\delta$ is on par with the squared Hellinger distance $\epsilon^2$:

Thus, \prettyref{lmm:reduction} reduces the problem to proving \[ \Delta \leq C \delta \frac{\log \frac{1}{\delta}}{\log\log \frac{1}{\delta}}. \]

Polynomial approximation and a Bernstein-style inequality

Let $f \equiv (f_G+f_H)/2= f_{(G+H)/2}$. Define the density ratio

where $\varphi$ is the standard normal density. Define the weight function

Crucially, we have

To estimate the differential operator on $g$, we first approximate $g$ by a polynomial, bound the operator norm of the differential operator on polynomials by a linear algebraic argument, then quantify the approximation error bound. These are carried out in the following two lemmas.

Complete the proof

With Lemmas (ref)--(ref), we now prove \prettyref{thm:ub-compact}.

Proof of lemmas

\paragraph{Proof of \prettyref{lmm:bernstein}}

Define the inner product weighted by $w$ as: \[ \langle f, g \rangle \equiv \int_\mathbb{R} f g w, \quad f,g\in L_2(w). \] Note that $w$ has at least a Gaussian tail so $L_2(w)$ includes all polynomials. Denote by $\{q_j: j\geq 0\}$ the system of orthonormal polynomials for $w$ such that \[ \langle q_i, q_j \rangle = {1\{{i=j}\}}. \] In particular, $\langle q_i, p \rangle=0$ for all polynomial of degree at most $i-1$. To be definitive, assume that the leading coefficient of each $q_i$ is (strictly) positive, denoted by $\kappa_j>0$.

Expand a degree-$k$ polynomial $p$ as $p=\sum_{j=0}^k a_j q_j$ and $p'= \sum_{j=0}^k b_j q_j$, where the coefficients are related by vector-matrix multiplication: \[ b = L a \quad \text{ and }\quad \|p'\|_{L_2(w)} \leq \|L\|_{\mathrm{op}} \|p\|_{L_2(w)}. \] where $ \|L\|_{\mathrm{op}}$ is the largest singular value of the $(k+1)\times (k+1)$ matrix $L$.

We thus analyze $L$ and bound ${\left\Vert L \right\Vert}_{\mathrm{op}}$. Note that $p'=\sum_{j=0}^k a_j q_j'$, and so $L= (L_ {ij})$ is given by expanding each $q_j'$, namely \[ q_j' = \sum_{i=0}^k L_{ij} q_i. \] By orthogonality, $L_{ij} = \langle q_i, q_j' \rangle$. Clearly, $L_{ij} = 0$ for all $i\geq j$; in other words, $L$ is upper triangular. Note that, by integration by parts, we have

for all $f,g$ with at most exponential growth. Therefore, \[ L_{ij} =

, \] since $\ip{q_i', q_j} = 0$ when $ i < j$.

To bound the operator norm of $L$, we decompose it into a dense part and a sparse part. Since $w = \varphi^2/f$, $f= f_\nu$ with $\nu=\frac{G+H}{2}$ supported on $[-M,M]$, we have

Then $L=A+B$, where \[ A_{ij} =

\quad B_{ij} =

\] where we used the fact that, by orthogonality, $\langle q_i, q_j y \rangle= \langle q_i y, q_j \rangle =0$ if $i<j-1$. It remains to bound the operator norms of $A$ and $B$ separately.

For $A$, we compute its Frobenius norm: \[ \|A\|_F^2= \sum_{j=1}^k \sum_{i=0}^{j-1} \langle q_i, q_j m_\nu \rangle^2 \leq \sum_{j=1}^k \|q_j m_\nu\|_{L_2(w)}^2 \leq M^2 \sum_{j=1}^k \| {q_j}\|_{L_2(w)}^2 = M^2 k. \] Thus $\|A\|_{\mathrm{op}} \leq M \sqrt{k}$.

Note that $B$ contains a single superdiagonal. Then \[ \|B\|_{\mathrm{op}} = \max_{j=1,\ldots,k} |a_j| \] where \[ a_j \equiv B_{j-1,j} = \langle q_{j-1}, q_j y \rangle = \ip{q_{j-1}y, q_j}. \] Equivalently, $a_j$ appears in the three-term recursion:

for some $b_j$. By assumption $q_j(y)=\kappa_j y^j +\ldots$ and $q_{j-1}(y)=\kappa_{j-1} y^j +\ldots$ with $\kappa_j>0$. Thus $q_j = \frac{\kappa_j}{\kappa_{j-1}}y q_{j-1} + p_{j-1}$ for some polynomial of degree at most $j-1$. Comparing leading coefficients, we thus get $a_j =\langle q_j , q_{j-1} y \rangle = \frac{\kappa_{j-1}}{\kappa_j}>0$. By a similar argument, $\langle q_j', q_{j-1} \rangle = \frac{j\kappa_j}{\kappa_{j-1}}$. Applying \prettyref{eq:adjoint} and \prettyref{eq:Vprime}, \[ \underbrace{\langle q_j', q_{j-1} \rangle}_{j/a_j} = -\underbrace{\langle q_j, q_{j-1}' \rangle}_{=0} + \langle q_j , q_{j-1} V' \rangle = \underbrace{\langle q_j , q_{j-1} m_\nu \rangle}_{\equiv \beta_j} + \underbrace{\langle q_j , q_{j-1}y \rangle}_{a_j} \] Thus we get $a_j(a_j+\beta_j) = j$, that is, $a_j = (-\beta_j + \sqrt{\beta_j^2 + 4 j})/2$. Applying Cauchy-Schwarz, \[ |\beta_j| = |\langle q_j , q_{j-1} m_\nu \rangle| \leq \|q_j\|_{L_2(w)} \|q_{j-1}\|_{L_2(w)} \|m_\nu\|_\infty \leq M. \] and hence $a_j \leq \sqrt{j} + M$. This shows \[ \|B\|_{\mathrm{op}} \leq \sqrt{k}+M \] and hence $\|L\|_{\mathrm{op}} \leq \sqrt{k}+M + M\sqrt{k}$, as desired.

\paragraph{Proof of \prettyref{lmm:poly-f}} Since $\nu=(G+H)/2$ is centered, we have by Jensen's inequality \[ f(y) = \mathbb{E}_{U\sim \nu}[\varphi(y-U)] \geq \mathop{\rm exp}(-\sigma^2/2) \varphi(y) \] where $\mathbb{E}_{U\sim \nu}[U^2] \leq \sigma^2$ by the assumption that the second moments of $G$ and $H$ are both at most $\sigma^2$. Thus

So it suffices to bound the approximation error of $g$ and $g'$ in $L_2(\varphi)$, under standard Gaussian.

It is well-known that the likelihood ratio between a Gaussian mixture and a standard Gaussian admits the following Hermite expansion (see, e.g., WY18): \[ \frac{f_G(y)}{\varphi(y)} = \sum_{j\geq 0} \frac{m_j(G)}{j!} H_j(y) \] where $H_j$ is the degree-$j$ Hermite polynomial satisfying $\mathbb{E}[H_i(Z)H_j(Z)]=i! {1\{{i=j}\}}$, and $m_j(G)$ is the $j$th moment of $G$. Hence \[ g(y) = \frac{f_G(y)-f_H(y)}{\varphi(y)} = \sum_{j\geq 0} \frac{m_j(G)-m_j(H)}{j!} H_j(y). \] Truncate this series at $k$ and denote it by $p(y)$. Then we also have $g'(y)-p'(y) = \sum_{j>k} \frac{m_j(G)-m_j(H)}{j!} H_j'(y) = \sum_{j>k} \frac{m_j(G)-m_j(H)}{(j-1)!} H_{j-1}(y)$. Thus \[ \|g-p\|_{L^2(\varphi)}^2 = \sum_{j>k} \frac{|m_j(G)-m_j(H)|^2}{j!}, \quad \|g'-p'\|_{L^2(\varphi)}^2 = \sum_{j>k} \frac{|m_j(G)-m_j(H)|^2}{(j-1)!} \] and it suffices to bound the latter. Note that $|m_j(G)|,|m_j(H)| \leq M^j$ and $j! \geq (j/e)^j$, we get

completing the proof by summing the geometric series.

Priors with exponential tails

In this section we prove \prettyref{thm:ub-exp}, extending the regret bound in \prettyref{thm:ub-compact} for compactly supported priors to those in the class ${\mathcal{G}}_\alpha(\sigma)$ defined in \prettyref{eq:Galpha}, having exponential tails. Throughout this section, all constants (including those implicit in $\lesssim$) depend on $\alpha,\sigma$.

Note that one may apply \prettyref{thm:ub-compact} by approximating the true prior with a compact-supported one. This yields a worse logarithmic factor. For example, for the subgaussian class, applying \prettyref{thm:ub-compact} with $M \asymp \sqrt{\log \frac{1}{\epsilon}}$ yields an $O\pr{\epsilon (\log \frac{1}{\epsilon})^2}$ upper bound, which is inferior to \prettyref{thm:ub-exp}. Instead, we take a more direct approach.

Similar to the proof of \prettyref{thm:ub-compact}, the first step is again to reduce $\mathsf{Regret}(H\|G)$ to a weighted derivative quantity $\Delta$. Recall the definitions of $g$, $w$, $\delta$, and $\Delta$ in ((ref))-- ((ref)) from \prettyref{sec:compact}.

The next lemma is an analogue of the polynomial approximation in \prettyref{lmm:poly-f}.

Finally, the following result is an analogue of \prettyref{lmm:bernstein}. The proof is deferred to the end of \prettyref{app:bernstein}.

Putting these together, we prove \prettyref{thm:ub-exp}.