Fast TV-Distance Estimation for Product Distributions

Fast TV-Distance Estimation for Product Distributions

Introduction

If you’ve ever compared two probability models and wondered “how different are they, really?”, you’ve run into total variation (TV) distance. It’s one of those core statistics that shows up everywhere—from measuring how hard it is to tell two worlds apart, to bounding errors in probabilistic reasoning. The catch is that computing (or even approximating) TV distance can get surprisingly nasty once you move beyond very small toy cases.

This blog post is based on new research from the original paper, where Konrad Anand, Alistair Benford, and Heng Guo tackle a key bottleneck: approximating TV distance between two product distributions in linear time (linear in the input size, for fixed accuracy parameters). The paper focuses on distributions of the form
- (P = \bigotimes{i=1}^n Pi) and (Q = \bigotimes{i=1}^n Qi) on ([q]^n),
meaning each coordinate is independent, but the marginals (Pi, Qi) can vary with (i).

Why is this nontrivial? Because TV distance does not tensorise nicely over products the way some other divergences do. In fact, the exact evaluation for product distributions is known to be #P-hard (cited in the paper as [BGM+25a]). So the field has moved toward approximation—especially relative-error approximation, since additive error can be too weak when the true distance is small.

The authors build on a line of work on randomized TV estimators (including the algorithm by Feng, Guo, Jerrum, and Wang [FGJW23], plus later improvements and derandomisations). But this note introduces a very different idea to get a faster runtime: an adaptation of filtered Monte Carlo (credited to [Gla93] in the paper), combined with a carefully designed unbiased estimator whose variance behaves nicely.

Why This Matters

This is significant right now because product distributions show up constantly in modern pipelines—even when people aren’t thinking in terms of “(q)-ary independent coordinates.” Think of any situation where you have (n) independent components (or you model them as such): per-token distributions in NLP approximations, independent feature models in recommender systems, compartmentalised uncertainty in sensor networks, or “factorised” generative models used for efficiency.

In those settings, TV distance is a natural knob: it’s directly tied to how distinguishable two probabilistic systems are. If TV distance between “model A” and “model B” is small, then downstream decisions that rely on sampling can’t reliably tell the systems apart. That’s exactly the kind of question you want to answer when you’re doing model auditing, robustness checks, or even privacy-style reasoning (“is an attacker’s view statistically close to the real one?”).

The practical kicker is runtime. The paper’s main result gives a randomized algorithm that outputs a nonnegative random variable (Z) such that the expectation matches TV distance, and the estimator concentrates well enough to get a multiplicative/relative-error approximation. For fixed (\varepsilon) and (\delta), the time is linear in the input size (qn). That matters when (n) is large and you’re given the distributions explicitly only through their marginals.

And compared to previous AI-adjacent research: AI communities often use relative-entropy or KL-style scores because of convenient decomposition properties. TV distance is harder precisely because it refuses to behave like KL. So what this work does is make TV distance computationally feasible again in a regime where independence holds—essentially giving TV the “fast approximation treatment” that KL had first (and that many ML workflows implicitly prefer).

Main Content: Linear-Time Relative-Error TV Estimation for Product Distributions

What TV Distance Means Here—and Why Independence Doesn’t Save You

Let’s ground the problem in intuition. For distributions (P) and (Q) over the same finite universe,
[
d{\mathrm{TV}}(P,Q)=\frac{1}{2}\sum{\omega} |P(\omega)-Q(\omega)|.
]
Equivalently, TV distance is the maximum advantage any (measurable) test can get in distinguishing samples from (P) versus (Q).

Now set up the model: you’re given product distributions on ([q]^n). Each coordinate (i) is sampled independently from (Pi) under (P), and independently from (Qi) under (Q). You might think this independence should make TV computation “manageable.” But the paper reminds us of the big structural obstacle: TV distance does not tensorise over product distributions. That means you can’t just combine coordinate-level differences in a straightforward dynamic-programming way.

Even worse, the paper points out that exact evaluation for product distributions is #P-hard, so unless you’re in small regimes (tiny (n) or tiny (q)), approximation is the right direction.

So the authors target relative error: they want an estimator (Z) such that (Z) approximates (d{\mathrm{TV}}(P,Q)) within a factor ((1\pm \varepsilon)) with high probability (failure probability (\delta)). Relative error is the right setting when (d{\mathrm{TV}}(P,Q)) might be small.

The Core Challenge: Unbiased Estimation with Bad Variance

A “first attempt” at unbiased estimation for TV often starts from mixture sampling. Consider the equal mixture:
[
\frac{P+Q}{2}.
]
If you sample (\omega) from ((P+Q)/2), you can construct an unbiased estimator using something proportional to (|P(\omega)-Q(\omega)|/(P(\omega)+Q(\omega))).

The paper explains this route and then immediately shows why it’s not enough: the relative variance of that estimator can be huge, especially when (d_{\mathrm{TV}}(P,Q)) is small. Intuitively, when the two distributions are very similar, “most samples” don’t carry informative difference signal; when the difference does show up, it’s rare—so your estimator becomes noisy.

This is a classic Monte Carlo pain point: unbiasedness alone doesn’t guarantee efficiency. You need an estimator whose randomness doesn’t explode right when the target quantity is small.

Filtered Monte Carlo: Turning One Noisy Estimator into Many Stable Increments

The main algorithmic idea is to replace the naïve estimator with something built from a filtered Monte Carlo framework (from [Gla93], adapted in this work).

Here’s the mental model the paper uses:

1. Introduce a hidden sign/hypothesis variable (H \in {+,-}) chosen fairly.
2. If (H=+), sample the product distribution (P); otherwise sample (Q).
   • Equivalently: the overall sampling is from ((P+Q)/2).
3. As you reveal the sample coordinate-by-coordinate (X=(X1,\dots,Xn)), you maintain “posterior” beliefs about whether the hidden world was (P) or (Q).
   • The paper uses quantities (\alphai) and (\betai) representing updated weights (after seeing coordinates up to (i)).
4. Define the “bias” (Ui = \alphai - \beta_i). This bias evolves as a martingale-like object (more precisely: the paper notes a bounded martingale / submartingale structure).
5. Instead of directly using (|Un|) (which would reintroduce high variance), you accumulate the predictable increases in (|Ui|) coordinate-by-coordinate:
   [
   A := \sum{i=1}^n Ci,
   ]
   where each (C_i) is computed from the information available after observing coordinates up to (i-1).

This choice—sum the “small informative increments” rather than measure the final magnitude directly—is the variance-control trick.

The paper states two crucial properties:
- Unbiasedness: (\mathbb{E}[A] = d{\mathrm{TV}}(P,Q)).
- A strong second-moment bound: (\mathbb{E}[A^2] \le 2\,d{\mathrm{TV}}(P,Q)^2).

That second bound is a big deal. It implies the relative variance is bounded (the paper derives (\mathrm{Var}(A)/(\mathbb{E}A)^2 \le 1)), which is exactly what you need for relative-error concentration.

How the estimator stays nonnegative and computable

The paper describes an implementation detail (its Algorithm 1) where the coordinate contribution (C_i) is computed using max(·,0) style logic. Practically: you can think of it as “how much the posterior bias is expected to move toward agreement/disagreement, but only counting the positive parts in a controlled way.”

That nonnegativity matters because the final estimator is meant to approximate a distance (which is nonnegative).

Putting the Concentration Layer on Top: Median of Means for Relative Error

Even with bounded relative variance, you still need a high-probability approximation. The paper uses the standard reliability tool:

• Run the estimator (s) times independently and average:
  [
  \overline{A} = \frac{1}{s}\sum_{k=1}^s A^{(k)}.
  ]
• Then repeat this group (r) times and take the median of the group means.

This “median of means” approach turns a constant-probability accuracy guarantee into a (1-\delta) guarantee.

The paper sets:
- (s = \left\lceil\frac{4}{\varepsilon^2}\right\rceil),
- and chooses (r) as the smallest odd integer with (r \ge 8\log(1/\delta)).

The outcome is: with probability at least (1-\delta), the reported estimate (Z) is within relative error (\varepsilon) of (d_{\mathrm{TV}}(P,Q)) (the paper frames it through the estimator properties plus Chebyshev/Chernoff style amplification).

The Main Runtime Result—and How It Beats Prior Work

Now for the part you probably care about most: how fast is this?

The theorem in the paper says: for product distributions (P=\bigotimes{i=1}^{n}Pi) and (Q=\bigotimes{i=1}^{n}Qi) on ([q]^n), with accuracy parameters (0<\varepsilon\le 1) and (0<\delta<1), there is a randomized algorithm producing a nonnegative estimator (Z) in time
[
O\left(qn\,\varepsilon^{-2}\log\frac{1}{\delta}\right).
]
Crucially, the input is given explicitly by the marginals, and for fixed (\varepsilon,\delta), this is linear in the input size (which scales like (qn)).

The paper also compares this to prior runtimes. Here’s the comparison, as presented in the introduction:

So the headline is: the authors’ filtered Monte Carlo estimator plus concentration gives a dramatic speedup over the earlier randomized estimator structure.

An alternate implementation that can be even faster in some regimes

The paper doesn’t stop at one implementation. It proposes an alternative way to compute the estimator increments using sorting and binary search per coordinate. That alternative runtime becomes:
[
O\left(nq\log q + \frac{n\log q}{\varepsilon^{2}}\log\frac{1}{\delta}\right),
]
and the paper notes this can be faster when (\varepsilon^{-2}\log q) and other terms behave in a favorable way relative to (q) (it states a condition comparing (\varepsilon^{-2}\log q) terms and ((\log q)/q) terms).

So, depending on your regime, you might pick:
- the “direct” linear-time estimator, or
- the “preprocess + search” version that trades time for faster per-run increments.

Limits: Why You Still Can’t Beat ( \Omega(qn))

A good theory result doesn’t just claim an algorithm is fast—it also answers: can anyone do asymptotically better?

The paper includes a lower bound (in a marginal-query model) showing that any randomized algorithm that achieves a multiplicative ((1\pm \varepsilon)) approximation with failure probability (\delta) must make at least (\Omega(qn)) queries in the worst case.

The proof idea is adversarial: they construct distributions where many coordinates look identical except one “needle” marginal differs. If your algorithm doesn’t query that marginal enough, it can’t reliably tell which case it’s in, forcing a lower bound on queries—and hence runtime.

The important takeaway for practitioners:
- linear time in input size is essentially optimal in the right model.

Key Takeaways

• TV distance is hard for products: even for independent (product) distributions, exact computation is #P-hard and TV doesn’t tensorise nicely like KL.
• Main contribution: a randomized linear-time relative-error approximation algorithm for TV distance between product distributions on ([q]^n), running in
  [
  O\left(qn\,\varepsilon^{-2}\log\frac{1}{\delta}\right).
  ]
• Estimator design matters: the algorithm uses a filtered Monte Carlo strategy that builds the estimator from stable, coordinate-by-coordinate increments rather than a single end value with huge variance.
• Concentration guarantee: by combining bounded second moment ((\mathbb{E}[A^2]\le 2\,d_{\mathrm{TV}}(P,Q)^2)) with a median-of-means wrapper, the estimator achieves relative error with probability at least (1-\delta).
• Practical implication: if your data or model is well-approximated by independent coordinates and you’re given marginals explicitly, you can now approximate TV distance efficiently even when (n) is large.
• Near-optimality: the paper proves an (\Omega(qn)) lower bound in a marginal-query model, meaning you shouldn’t expect asymptotically faster algorithms in that setting.

If you tell me your target regime (typical values of (q), (n), and whether (\varepsilon) is “loose” or “strict”), I can help you interpret which of the two implementations (direct vs sorted/binary-search) is likely the better fit in practice.