Fast computation and theoretical guarantees for the NPMLE in exponential family mixtures

Abstract: This work makes two advances in the study of the (approximate) nonparametric maximum likelihood estimator (NPMLE) for exponential family mixture models. First, we develop a data-compression strategy that reduces the cost of repeated likelihood evaluations in NPMLE computation to logarithmic order in the sample size. Second, we show that, for a broad class of approximate NPMLEs, the resulting marginal density estimator attains an almost parametric rate of convergence.

• The paper introduces a novel quadrature-based data compression method to reduce likelihood evaluation from linear to polylogarithmic complexity.
• It establishes near-parametric convergence rates for the NPMLE in exponential family mixtures under compact and unbounded support scenarios.
• Empirical studies demonstrate over two orders of magnitude speedup with maintained estimation accuracy compared to traditional interior-point and augmented Lagrangian methods.

Fast Computation and Theoretical Guarantees for the NPMLE in Exponential Family Mixtures

Introduction and Context

This paper addresses core methodological and theoretical challenges in the computation and analysis of the nonparametric maximum likelihood estimator (NPMLE) for mixture models where component densities belong to an exponential family. The context encompasses a broad range of statistical applications—such as empirical Bayes, principal component analysis, and large-scale multiple testing—where flexible yet principled modeling of heterogeneity via exponential family mixtures is essential. While the favorable adaptivity and asymptotic properties of the NPMLE are well-documented in specific settings (e.g., Gaussian and Poisson mixtures), the computational bottleneck posed by costly likelihood evaluations and the lack of a unified theory for near-parametric convergence rates in general exponential families have remained significant open problems.

Computational Framework: Data Compression for Likelihood Evaluation

A major contribution of this work is a generic data-compression strategy that dramatically reduces the computational complexity of repeated likelihood evaluations in the NPMLE optimization problem. The classical approach incurs at least linear cost in the sample size $n$ due to the necessity of evaluating the observed-data log-likelihood,

$\ell_n(f) = \sum_{i=1}^n \log f(X_i),$

repeatedly over candidate mixing distributions. The authors generalize the classical "grouping trick" for discrete models (where empirical histograms suffice) to arbitrary (possibly continuous) exponential families by proposing a moment-matching quadrature-based compression:

1. Gaussian Quadrature Compression: The empirical distribution $P_n$ of the data is replaced with a discrete measure $P_{J_n}$, supported on $J_n \ll n$ points, matching the first $2J_n-1$ moments of $P_n$ exactly.
2. Approximate NPMLE: The NPMLE is then computed with log-likelihood replaced by expectations under $P_{J_n}$, thereby reducing repeated likelihood evaluation costs to be polylogarithmic in $n$.

The procedure is encapsulated in Algorithm 1 (see main text for details).

Figure 1: Comparison of IPM, Ours, and ALM in terms of computation time, log-likelihood, and sum of squared errors in empirical Bayes estimation.

Figure 1 demonstrates substantial computational speedups—over two orders of magnitude—compared with both the canonical interior-point approach [koenker2014convex] and the more recent augmented Lagrangian method [zhang2024efficient]. Importantly, these speedups are achieved with statistically indistinguishable empirical Bayes estimation accuracy and log-likelihood values.

Theoretical Guarantees: Parametric-Like Marginal Density Estimation

The second main contribution is a unified, scenario-dependent theoretical analysis of the plug-in marginal density estimator $f_{\hat{g}}$ when either the exact or a quadrature-based approximate NPMLE is used:

• Scenarios: The theory distinguishes two settings:
  • (S1): Compact support for the mixing distribution; $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$0 may have exponential-like tails.
  • (S2): Unbounded support for $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$1, but with polynomial-type regularity and tail constraints on $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$2.
• Light Tail Condition: For near-parametric rates, the observation's marginal distribution must be sufficiently light-tailed, a condition made precise in assumption (LT) and automatically satisfied under (S1).

Likelihood Approximation and Choice of Quadrature Order

The core theoretical insight is a quantification (Theorem 1) of the number of quadrature points $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$3 needed so that the likelihood computed under $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$4 approximates the true empirical likelihood within a prescribed tolerance $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$5. For both scenarios, only polylogarithmic $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$6 in $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$7 suffices to ensure this, under moderate growth of the sample maxima.

Rate Results for Marginal Density Estimation

• Main Theorem (S1): If an approximate NPMLE achieves likelihood within $\ell_n(f) = \sum_{i=1}^n \log f(X_i),$8 of the optimum, then with high probability the squared Hellinger distance between the estimated and true mixture marginal densities satisfies

$\ell_n(f) = \sum_{i=1}^n \log f(X_i),$9

This rate is near-parametric, up to logarithmic factors, and is sharper than previous bounds in closely related settings.

• Main Theorem (S2): For more general (possibly unbounded) parameter spaces and suitably regular $P_n$0, the same structure of result holds,

$P_n$1

with a model-dependent logarithmic exponent $P_n$2, again attaining almost parametric rates provided logarithmic control of likelihood loss.

The analysis leverages complex-analytic smoothness arguments (effective polynomial approximation in parameter space), and bounding of covering numbers (entropy calculations) for the induced class of marginal densities or likelihood ratio functions. These arguments replace the more customary entropy-based approach with parametric bracketing numbers, which can be prohibitive in irregular settings.

Empirical Results

Simulation studies on large $P_n$3 Gaussian location mixture models empirically verify that the proposed method achieves substantial computation time reductions (more than two orders of magnitude), matching the log-likelihood and empirical Bayes risk of leading optimization-based solvers. The data compression step, which is the principal computational bottleneck, still markedly outperforms previous methods and remains feasible for very large datasets with sizes $P_n$4 (see Figure 1).

Extensions

The proposed methodology adapts straightforwardly to higher-dimensional and heteroscedastic settings, such as empirical Bayes analysis under the heteroscedastic Gaussian sequence model. For multivariate quadrature, Carathéodory–Tchakaloff subsampling is advocated as a practical mechanism to match all polynomial moments up to a prescribed degree with a compressed support, guaranteeing similar likelihood and risk bounds.

Implications and Prospects

This framework yields both practical and theoretical advances:

• Computational Implication: By reducing the computational complexity of likelihood evaluations from linear to polylogarithmic in $P_n$5, the method unlocks routine scalable NPMLE computation for previously prohibitive large-scale mixture modeling tasks.
• Statistical Implication: The theoretical results consolidate and generalize previously model-specific convergence rate analyses, demonstrating that (approximate) NPMLEs in exponential family mixtures adapt to the statistical regularity of the underlying component family and mixing distribution, inheriting near-parametric estimation rates even in high-dimensional, nonparametric settings.
• Future Directions: Open questions remain regarding extension to non-exponential families, optimal bounds for mixture classes with model misspecification, and the study of frequentist inference (e.g., likelihood-based testing) under compressed likelihoods. Further research may investigate refined adaptive bounds when the true mixing distribution is discrete or otherwise low-complexity, as well as computational gains in higher-dimensional inference problems.

Conclusion

This work rigorously addresses and resolves computational and theoretical difficulties in NPMLE for exponential family mixture models by (1) developing a provably fast, data-compressed likelihood evaluation strategy based on moment-matching quadrature and (2) establishing near-parametric convergence rates for approximate NPMLEs under weak regularity conditions. The results furnish both a substantial improvement in the scale at which NPMLE-based empirical Bayes and density estimation can be deployed, and a unified analytic framework for future theoretical investigations. The methodology is broadly applicable and is expected to facilitate further advances in adaptive inference under flexible models.