📊 statistics

A novel sequential method for building upper and lower bounds of moments of distributions

This paper introduces a novel sequential method based on the majorization-minimization framework and power diagrams to construct convergent upper and lower bounds for the moments of unnormalized distributions, addressing the critical need for guaranteed inequality preservation in numerical integration tasks across one and multi-dimensional settings.

Original authors: Solal Martin, Emilie Chouzenoux, Victor Elvira

Published 2026-05-06

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Solal Martin, Emilie Chouzenoux, Victor Elvira

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to measure the total amount of water in a very strange, irregularly shaped lake. You can't just pour it into a bucket to measure it (the math is too hard), and you can't see the bottom clearly. This is the problem statisticians face when they need to calculate "moments" (like averages or variances) of complex probability distributions. Usually, they use methods that give a single "best guess" number, but they don't know for sure if that guess is too high or too low.

This paper introduces a new method that acts like a smart, shrinking cage to catch the true answer. Instead of guessing a single number, the authors build a lower fence and an upper fence around the true value. They guarantee that the real answer is somewhere between these two fences, and they have a way to make the fences get closer and closer together until they are almost touching the truth.

Here is how the method works, using simple analogies:

1. The "Tangent Trampoline" (Building the Fences)

Imagine the shape of the lake (the distribution) is a bumpy hill. The authors want to cover this hill with a blanket (an upper bound) and lift a floor underneath it (a lower bound).

The Trick: They pick a specific point on the hill and place a Gaussian curve (a smooth, bell-shaped curve) there.
The Lower Fence: They place a bell curve under the hill that touches it at that one point but never goes above it. It's like a trampoline that supports the hill from below.
The Upper Fence: They place a bell curve over the hill that touches it at that same point but never goes below it. It's like a ceiling that hugs the hill from above.

2. The "Patchwork Quilt" (Refining the Fences)

One bell curve isn't enough to cover a bumpy hill perfectly; there will be big gaps between the curve and the hill.

The Solution: The authors pick many points on the hill. At each point, they place a new lower trampoline and a new upper ceiling.
The Envelope: They then stitch these together. The final "lower fence" is the highest point of all the trampolines combined (like a quilt made of the best parts of each trampoline). The final "upper fence" is the lowest point of all the ceilings combined.
The Result: This creates a "patchwork" shape that hugs the hill much more tightly than a single curve could.

3. The "Smart Gardener" (Iterative Improvement)

How do they know where to put the next point?

They look at the gaps between their current fences. Where is the gap the widest? That's where the approximation is worst.
They plant a new "tangency point" (a new trampoline/ceiling) right in the middle of that widest gap.
They repeat this process. Every time they add a new point, the fences get tighter, and the gap between the upper and lower bounds shrinks.

4. The "Guaranteed Catch" (Why It Matters)

Most computer methods give you a number and say, "I'm pretty sure this is right." This method says, "I guarantee the answer is at least X and at most Y."

Why is this useful? In the paper, they use this to check the reliability of a statistical tool called Importance Sampling.
The Analogy: Imagine you are trying to estimate how many fish are in a lake using a net. You want to know if your net is good. Usually, you just run the simulation and get a number. But sometimes, rare events (like a giant fish jumping) mess up your calculation.
The Paper's Contribution: By using their "shrinking cage" method, they can calculate the variance (the instability) of that net without running millions of expensive simulations. They can prove, mathematically, that the error is small. This helps them tune the net (the "proposal distribution") to catch fish more efficiently.

5. Going from 2D to 3D (The Multi-Dimensional Leap)

The paper also shows how to do this when the "lake" has more than one dimension (like a 3D volume instead of a 2D area).

Instead of simple bell curves, they use Power Diagrams. Imagine slicing a block of cheese with laser beams. The lasers cut the space into different polyhedral (multi-sided) cells.
Inside each cell, the math becomes simple enough to calculate. They sum up the results from all the cells to get the total bounds.

Summary

The paper presents a sequential, iterative method to build certified upper and lower bounds for difficult integrals.

It uses Gaussian curves to create initial fences.
It stitches them together into piecewise envelopes.
It iteratively adds points where the error is largest to tighten the bounds.
It proves mathematically that these bounds will eventually converge to the exact answer.
It demonstrates this by accurately estimating the variance of a statistical sampler, allowing researchers to tune their models for better performance without needing massive amounts of computing power.

The authors emphasize that this method is deterministic (it always gives the same result for the same inputs) and provides rigorous guarantees that the true value lies within the calculated range, which is a feature missing from many standard numerical methods.

Technical Summary: A Novel Sequential Method for Building Upper and Lower Bounds of Moments of Distributions

Problem Statement
The paper addresses the fundamental challenge of approximating intractable integrals of the form $I = \int_{-\infty}^{+\infty} f(x)\pi(x)dx$ , where $\pi(x) = \exp(-\phi(x))$ is an unnormalized distribution and $f(x)$ is a test function (specifically a monomial $x^k$ ). While standard numerical integration techniques (e.g., Gaussian quadrature, Monte Carlo) provide point estimates, they often fail to guarantee the preservation of majoration/minoration inequalities. This limitation is critical in applications such as Bayesian inference and importance sampling (IS), where the variance of an estimator is itself an intractable integral. In these contexts, point estimates may lack reliable uncertainty quantification, and the variance estimation is often unstable due to rare events. The authors propose a method to compute certified, deterministic upper and lower bounds for such integrals within a given computational budget, rather than merely seeking accurate point estimates.

Methodology
The proposed approach is a sequential, iterative framework based on the Majorization-Minimization (MM) principle and envelope construction.

Gaussian Tangent Envelopes:
The method relies on structural assumptions on the function $\phi(x)$ :
- Assumption 1 (Majorization): $\phi(x)$ is bounded above by a quadratic function at any point $t$ , allowing the construction of a tangent majorant for $\pi(x)$ in the form of an unnormalized Gaussian density.
- Assumption 2 (Minorization): $\phi(x)$ is bounded below by a quadratic function at any point $t$ , allowing the construction of a tangent minorant for $\pi(x)$ , also Gaussian.
  These assumptions hold for a wide class of functions, including log-concave distributions and those with Lipschitz-smooth gradients.
Piecewise Envelope Construction:
Instead of using a single tangent point, the algorithm maintains a set of tangency points $T = \{t_1, \dots, t_M\}$ .
- The lower envelope $C(x; T)$ is defined as the pointwise maximum of the tangent minorants.
- The upper envelope $\bar{C}(x; T)$ is defined as the pointwise minimum of the tangent majorants.
- These envelopes are piecewise Gaussian. The integration of these envelopes against the test function $f(x)$ (assumed to be a monomial) reduces to computing moments of (possibly truncated) Gaussian distributions, which can be calculated recursively with high precision.
Adaptive Refinement (Algorithm 2):
The algorithm iteratively refines the bounds by selecting new tangency points.
- At each iteration, the gap between the upper and lower integral bounds is computed over sub-intervals defined by the current tangency points.
- A new tangency point is added in the sub-interval where the approximation gap is maximal.
- To ensure theoretical convergence, the selection of new points is constrained to a dyadic candidate pool within a compact interval derived from the tails of the initial Gaussian majorant. This prevents the algorithm from neglecting regions of the domain.
Multidimensional Generalization:
For $d > 1$ , the method generalizes by assuming the curvature matrices of the quadratic bounds are diagonal (Assumption 4). This allows the Gaussian bounds to factorize coordinate-wise. The construction of the piecewise envelopes in higher dimensions is mapped to a Power Diagram (weighted Voronoi diagram) problem. The integration domain is partitioned into convex polyhedral cells, and the integral is computed as a sum of integrals over these cells using numerical quadrature.

Key Contributions

Sequential Bound Construction: Introduction of a novel algorithm that generates sequences of upper and lower bounds converging to the true integral value, rather than a single point estimate.
Theoretical Guarantees: Proof of convergence for the bound sequences under mild conditions (compact refinement of the tangency grid). The paper establishes that the gap between bounds vanishes asymptotically.
Controlled Accuracy: The method provides explicit error control, allowing users to stop when a desired precision $\epsilon$ is reached.
Variance Estimation for Importance Sampling: A specialized application demonstrating how to bound the variance of an IS estimator. This enables the tuning of proposal distributions (e.g., finding the optimal standard deviation $\theta$ ) without relying on unstable Monte Carlo variance estimates.
Efficient Multidimensional Implementation: A practical strategy for $d$ -dimensions using power diagrams to handle the geometric complexity of piecewise Gaussian envelopes.

Results
The authors validate the method through numerical experiments:

1D Case: In a Bayesian logistic regression setting, the method computes tight bounds for the normalization constant, first moment, and second moment. The results show rapid monotone convergence, with the bounds enveloping the ground truth (computed via MATLAB's integral function) with relative errors as low as $10^{-4}$ using a modest number of tangency points (e.g., ~100 points).
Variance Tuning: The method successfully identifies the optimal proposal variance $\theta$ that minimizes the IS estimator's variance, matching empirical results from $10^6$ Monte Carlo runs but without the computational cost or instability of MC variance estimation.
Comparison with State-of-the-Art: Compared to the envelope-based method by Evans et al. [35], the proposed method achieves comparable or superior accuracy with fewer evaluation points when using low-order approximations. While Evans' method can achieve higher accuracy with high-order derivatives ( $d=2$ ), it incurs significantly higher computational complexity due to root-finding for inflection points. The proposed method avoids high-order derivatives, offering a robust alternative.
2D Case: The method is extended to a 2D logistic regression problem. The use of power diagrams allows for efficient computation of bounds, demonstrating fast convergence of the integral gap and effective identification of optimal proposal parameters.

Significance
The paper claims significance primarily in providing deterministic guarantees for integral approximation, a feature often missing in standard numerical methods. By ensuring that the computed bounds strictly contain the true integral value, the method offers a reliable tool for:

Uncertainty Quantification: Providing explicit bounds on estimator quality (e.g., variance) in Bayesian inference, which is crucial when rare events dominate the error.
Hyper-parameter Tuning: Allowing for the systematic optimization of proposal distributions in Importance Sampling without expensive Monte Carlo simulations.
Robustness: Offering a convergence-guaranteed alternative to methods that may fail to preserve inequalities or require high-order derivatives that are difficult to compute.

The authors maintain a modest stance, noting that the method relies on specific structural assumptions (log-concavity or specific curvature bounds) and that future work could explore non-log-concave targets and other proposal families (e.g., Student's t).