arXiv CSopen access

On Computing Total Variation Distance Between Mixtures of Product Distributions

On Computing Total Variation Distance Between Mixtures of Product Distributions

arXiv:2605.03839v1 [cs.DS] 5 May 2026

Weiming Feng ∗

Yucheng Fu ∗

Minji Yang ∗

Anqi Zhang †

Abstract We study the problem of approximating the total variation distance between two mixtures of product distributions over an n-dimensional discrete domain. Given two mixtures P and Q with k1 and k2 product distributions over [q]n , respectively, we give a randomized algorithm that approximates dTV (P, Q) within a multiplicative error of (1 ± ε) in time poly((nq)k1 +k2 , 1/ε). We also study the special case of mixtures of Boolean subcubes over {0, 1}n . For this class, we give a deterministic algorithm that exactly computes the total variation distance in time poly(n, 2O(k1 +k2 ) ), and show that exact computation is #P-hard when k1 + k2 = Θ(n).

1

Introduction

Let P and Q be two discrete distributions over a sample space Ω. The total variation distance (TV-distance) between P and Q is defined as dTV (P, Q) =

X 1 X |P(ω) − Q(ω)| = max{0, P(ω) − Q(ω)}. 2 ω∈Ω ω∈Ω

The TV-distance is one of the most fundamental metrics for measuring the difference between distributions, as it characterizes the optimal distinguishability between P and Q. Computing the TV-distance is therefore a basic problem arising in learning and testing. A line of research has investigated the problem of approximating the TV-distance between two high-dimensional distributions, which is particularly interesting when the two distributions admit succinct representations. Early works [SV03, CR14, CK14, Kie18, BGMV20] studied the algorithms and complexity of approximating the TV-distance with additive error. More recently, [BGM+ 23] established that even for two product distributions, the exact computation of the TVdistance is #P-hard. In the same paper, the authors initiated the study of approximating the TV-distance with relative error, a more challenging task than the classical additive-error approximation. Subsequently, polynomial-time approximation algorithms have been developed for product distributions [BGM+ 23, FGJW23, FLL24], as well as for many high-dimensional distributions defined by local interactions, including Markov chains [FLL24], Bayesian networks [BGM+ 24], and undirected graphical models [FLY25]. ∗

School of Computing and Data Science, The University of Hong Kong. Emails: [email protected], [email protected], and [email protected]. † The Institute for Interdisciplinary Information Sciences, Tsinghua University. Email: [email protected].

1

Beyond local interactions, it is natural to consider more general high-dimensional distributions with inherently non-local structure. A fundamental and extensively studied example is given by mixtures of product distributions. Let [q] = {1, 2, . . . , q} and let Ω = [q]n denote the n-dimensional N (s) sample space. Fix two integers k1 , k2 ≥ 1. For each s ∈ [k1 ], let P(s) = ni=1 Pi be a product (s) distribution over Ω, where Pi denotes the marginal on the i-th coordinate, so that the coordinates are P 1 mutually independent under P(s) . Let α(1) , α(2) , . . . , α(k1 ) be mixing weights satisfying ks=1 α(s) = 1. The resulting mixture distribution is written as P=

k1 X

α(s)

s=1

n O (s)

Pi .

i=1

More explicitly, P is the distribution over Ω = [q]n given by ∀ω ∈ Ω,

P(ω) =

k1 X

α(s) P(s) (ω) =

s=1

k1 X

α(s)

s=1

n Y (s)

Pi (ωi ),

(1)

i=1

Mixtures of product distributions form a fundamental class of latent-variable models and have been extensively studied in learning theory [FOS08, JO14, CM19, GMRS21, GJM+ 24]. To draw a random sample X ∼ P, one first samples a hidden component index S ∈ [k1 ] according to N (S) Pr [S = s] = α(s) , and then samples X ∼ P(S) = ni=1 Pi from the product distribution P(S) . Although each component P(s) is itself a product distribution, the hidden index S induces a highly non-local dependency that simultaneously couples all coordinates of the sample X = (X1 , X2 , . . . , Xn ), which cannot be captured by any sparse local structure. In this paper, we study the problem of approximating the TV-distance between two mixtures of product distributions to within an arbitrarily small relative error, formally stated as follows. Problem 1.1. Approximate the TV-distance between two mixtures of product distributions. • Input: An error bound ε > 0 and two mixtures of product distributions over Ω = [q]n : P=

k1 X

α(s)

s=1

n O (s)

Pi

and Q =

k2 X t=1

i=1

β (t)

n O

(t)

Qi .

i=1

The distribution P is specified by the k1 mixing weights α(1) , α(2) , . . . , α(k1 ) together with the k1 product distributions P(1) , P(2) , . . . , P(k1 ) , where each P(s) is described by its n marginals (s) (s) (s) P1 , P2 , . . . , Pn over [q]. The distribution Q is given analogously. The total input size of the two mixtures is O((k1 + k2 )nq). • Output: A number dˆ such that (1 − ε)dTV (P, Q) ≤ dˆ ≤ (1 + ε)dTV (P, Q) . Due to the non-local nature of these dependencies among all coordinates, existing techniques for TV-distance approximation cannot be directly applied to mixtures of product distributions. The only prior progress is a polynomial-time algorithm for deciding whether dTV (P, Q) = 0 or dTV (P, Q) > 0 [BGM+ 25], leaving the approximation problem open. We partially resolve this open question by giving an FPRAS for the TV-distance between two mixtures of product distributions, provided that the numbers of components k1 and k2 are constants. 2

Theorem 1.2. Let k1 , k2 ≥ 1 be constants. There exists an FPRAS for Problem 1.1: given any P 1 N (s) ε > 0 and any pair of mixtures of product distributions P, Q over [q]n , where P = ks=1 α(s) ni=1 Pi Pk2 (t) Nn (t) has k1 components and Q = t=1 β i=1 Qi has k2 components, it solves the problem in time 2

2(k1 +k2 −1)

O( q (nq) ε2

) with success probability at least 99%.

Our algorithm is a Monte Carlo algorithm based on a carefully designed estimator, and we provide an overview of it in Section 2. Since the number of components k1 + k2 appears as an exponent in the running time, we require k1 + k2 to be a constant. A similar dependence on the number of components also arises in the running time of certain learning tasks for mixtures of product distributions [FOS08, CM19, GMRS21]. We leave it as an open question to understand the computational complexity when the number of components grows with n. Our second result concerns mixtures of Boolean subcubes, a special case of mixtures of product P N (s) distributions. A distribution P = ks=1 α(s) ni=1 Pi over {0, 1}n is called a mixture of Boolean (s) (s) subcubes if Pi (1) ∈ {0, 12 , 1} for every s ∈ [k] and i ∈ [n], where Pi denotes the marginal distribution of the i-th coordinate in the s-th component. Equivalently, for each component s, the (s) product distribution P(s) fixes every coordinate i with Pi (1) = 1 to 1, fixes every coordinate i with (s) Pi (1) = 0 to 0, and is uniform on the remaining coordinates. Mixtures of Boolean subcubes have been extensively studied in learning theory [CM19, BLMT23, BLST25] and capture many natural distributions; for instance, the uniform distribution over the satisfying assignments of a decision tree is a mixture of Boolean subcubes [FOS08]. In contrast to the general case, where we only obtain an FPRAS when k1 + k2 = O(1), we show that the TV-distance between two mixtures of Boolean subcubes P and Q can be computed exactly in time polynomial in the dimension n, provided that k1 + k2 = O(log n). Problem 1.3. Compute exactly the total variation distance between two mixtures of Boolean P 1 N P 2 (t) Nn (s) (t) subcubes P = ks=1 α(s) ni=1 Pi and Q = kt=1 β i=1 Qi . • Input: two mixtures of Boolean subcubes, given in the same format as in Problem 1.1. • Output: the exact value of dTV (P, Q) Theorem 1.4. There exists a deterministic algorithm for Problem 1.3: given any pair of mixtures P 1 N (s) of Boolean subcubes P, Q over {0, 1}n , where P = ks=1 α(s) ni=1 Pi has k1 components and P 2 (t) Nn (t) k1 +k2 · n) time. Q = kt=1 β i=1 Qi has k2 components, it solves the problem in O((k1 + k2 )3 The running time in Theorem 1.4 is polynomial in n whenever k1 + k2 = O(log n), yielding a deterministic polynomial-time algorithm in this regime. This does not contradict the #P-hardness of exact TV-distance computation established in [BGM+ 23]: their hardness already holds for two product distributions with arbitrary marginals, whereas our algorithm crucially exploits the restriction that each marginal of a Boolean subcube takes value in {0, 12 , 1}. The idea of the algorithm is as follows. For each component P(s) (and similarly for each Q(t) ), (s) the value of P(s) (ω) is either 0 or 2−r , where r is the number of coordinates i ∈ [n] with Pi (1) = 12 . Indeed, if some coordinate is fixed to be 1 (or 0) but ωi = 0 (or 1), then P(s) (ω) = 0; otherwise, each free coordinate contributes a factor of 1/2, so P(s) (ω) = 2−r . Thus, although there are 2n configurations ω ∈ {0, 1}n , the possible values of |P(ω) − Q(ω)| are at most 2k1 +k2 . Our algorithm enumerates these possible values of |P(ω)−Q(ω)| and counts the number of configurations ω ∈ {0, 1}n that achieve each one. The details are given in Section 5. 3

Finally, we complement this positive result by showing that the dependence on the number of components cannot be removed entirely: when the number of components grows linearly with n, exact computation becomes #P-hard. Theorem 1.5. If the number of mixture components is k1 + k2 = Θ(n), then the exact computation of the total variation distance in Problem 1.3 is #P-hard. The theorem is proved by a reduction showing that if we can exactly compute the TV-distance between two mixtures of Boolean subcubes with k1 + k2 = Θ(n) in polynomial time, then we can also exactly compute #3SAT in polynomial time. The detailed proof is in Section 6. Open problems As discussed after Theorem 1.2, the dependence on the number of components in the running time is not yet well understood. [BGM+ 25] gave an algorithm that decides whether dTV (P, Q) = 0 in time poly(n, k1 + k2 ), whereas the running time of our approximation algorithm depends exponentially on k1 + k2 . The analogous regime k1 + k2 = ω(1) is also not well understood for learning tasks on mixtures of product distributions, where obtaining algorithms with polynomial dependence on the number of components remains a central challenge [FOS08, CM19, GMRS21]. A natural question is whether one can approximate dTV (P, Q) in time polynomial in k1 + k2 , or prove a hardness result ruling this out when the number of components grows with n. As a first step toward the latter, Theorem 1.5 already shows that exact computation is #P-hard when k1 + k2 = Θ(n), and one direction is to extend its proof to rule out approximation as well. One can ask analogous TV-distance approximation questions for other structured mixture models. Natural candidates include mixtures of Gaussians [ABH+ 18, LS17], mixtures of Markov chains and MDPs [KTT23]. These models also give succinct descriptions of high-dimensional distributions and are widely studied in machine learning, but their dependence structures are more complex than that of mixtures of product distributions. More generally, it would be interesting to understand the complexity of computing or approximating the TV-distance between trajectory distributions induced by reward-mixing Markov decision processes (RMMDPs). RMMDPs use latent mixture structure to model the reward context, and TV-distance bounds between trajectory distributions have been used to learn near-optimal policies [KECM23]. The problem of computing the TV-distance between RMMDPs appears to require new ideas. At last, our algorithm for mixtures of product distributions is randomized. It remains open to design a deterministic algorithm for approximating the TV-distance between mixtures of product distributions. Currently, deterministic algorithms are only known for approximating the TV-distance between product distributions [FLL24, BGM+ 23].

2

Technical Overview

We present the main ideas behind the algorithm in Theorem 1.2. An obstacle for approximating the TV-distance is that the TV-distance itself can be very small. To see the reason, for two distributions P and Q over the sample space Ω, where Q is absolutely continuous w.r.t. P, one can write 

Q(ω) dTV (P, Q) = P(ω) max 0, 1 − P(ω) ω∈Ω X

n

o



= Eω∼P [g(ω)] ,

where g(ω) := max 0, 1 − Q(ω) P(ω) . We can efficiently sample ω and compute g(ω), so that a direct Monte Carlo estimate of dTV (P, Q) = Eω∼P [g(ω)] can only achieve additive error in polynomial 4

time. However, since the expectation Eω∼P [g(ω)] itself can be exponentially small, this simple approach is not sufficient to achieve a relative error approximation required by Problem 1.1.

2.1

Approximating the TV-distance via a coarse coupling

A framework for approximating the TV-distance was proposed in [FGJW23] for two product distributions (i.e., k1 = k2 = 1 in our setting). It can be generalized to any two distributions P and Q. Instead of approximating dTV (P, Q) directly, they estimated the normalized quantity ˜ where d˜ is a polynomial-time computable coarse estimator. Moreover, if the ratio dTV (P, Q) /d, dTV (P, Q) /d˜ is lower bounded by 1/poly(n), one can try to achieve a relative-error approximation to the ratio in polynomial time. Hence, we can recover dTV (P, Q) by multiplying d˜ by the estimate of the ratio. The coarse estimator d˜ can be constructed via a coupling. For two distributions P and Q over Ω, a coupling is a pair of joint random variables (X, Y ) such that X ∼ P and Y ∼ Q. It is well known that dTV (P, Q) = PrO [X ′ ̸= Y ′ ], where O = (X ′ , Y ′ ) is the optimal coupling between P and Q. Suppose we can construct a coupling C = (X, Y ) that is not far from the optimal coupling O. Formally, PrO [X ′ ̸= Y ′ ] = γ · PrC [X ̸= Y ] for some γ = 1/poly(n). Then d˜ = PrC [X ̸= Y ] is a ˜ define a distribution π: coarse estimator of dTV (P, Q). To estimate the ratio dTV (P, Q) /d, ω ∈ Ω,

π(ω) = PrC [X = ω | X ̸= Y ]

(2)

and define function f such that ω ∈ Ω,

f (ω) =

max{0, P(ω) − Q(ω)} . PrC [X = ω ∧ X ̸= Y ]

dTV (P,Q) 1 We can verify that Eπ [f ] = Pr = dTV (P,Q) = γ ≥ poly(n) . Using the approach described d˜ C [X̸=Y ] above, if we can perform the following operations in polynomial time

(a) sample ω ∼ π; (b) compute f (ω) for any ω; (c) compute PrC [X ̸= Y ]; then we can get a relative-error approximation of dTV (P, Q). Specifically, we draw T = poly( γ1 ) = P poly(n) independent samples ω1 , . . . , ωT ∼ π and compute f¯ = T1 Ti=1 f (ωi ). By Chebyshev’s inequality, with high probability, f¯ · PrC [X ̸= Y ] is an estimate of dTV (P, Q) within relative error.

2.2

Recursive coupling for mixtures of product distributions

It remains to construct a coupling that supports the three operations above. For product distributions, [FGJW23] used a coordinate-wise greedy coupling C = (X, Y ), which couples each coordinate optimally and independently, and showed that 1/n ≤ Eπ [f ] ≤ 1. This approach does not directly extend to mixtures of product distributions, where all the coordinates are correlated, so one cannot optimally couple all coordinates simultaneously. Our main technical contribution in this part is a new explicit recursive coupling for mixtures of product distributions, together with an algorithmic implementation of the three operations 5

above. At a high level, the recursive structure draws inspiration from [KECM23], who introduced a recursive argument to bound the TV-distance between two mixtures of product distributions. Their analysis, however, is non-constructive: it yields a fixed-error bound on the TV-distance, but does not by itself provide the algorithmic interface required here. To obtain a relative-error approximation with arbitrary ε, we need a coupling whose failure probability can be computed, and whose failure-conditioned distribution can be sampled and evaluated. We design such a coupling tailored to mixtures of product distributions, adapt the recursive analysis of [KECM23] to this explicit construction, and show that the three operations in (a), (b), and (c) can all be implemented in polynomial time. P 1 N (s) α(s) ni=1 Pi We construct a coupling between two mixtures of product distributions, P = ks=1 Pk2 (t) Nn (t) and Q = t=1 β i=1 Qi . Although coupling the first coordinate would bias the marginal (s) (t) distribution of the remaining coordinates, we can still extract the common parts of P1 and Q1 , and reduce the problem to coupling another mixture of product distributions on the remaining coordinates. For each c ∈ [q], we first compute a lower bound for each active mixture component: (

ℓ (c) ≜ min (s)

)

min

s∈[k1 ]:α(s) >0

(s) P1 (c),

min

t∈[k2 ]:β (t) >0

(t) Q1 (c)

(3)

,

(t)

where P1 (resp. Q1 ) is the marginal distribution of the first coordinate of the s-th component of P(s) (resp. Q(t) ). For ω ∈ [q]n with ω1 = c, we can write the probability of ω as P(ω) = ℓ (c)  α(s) ≜

α(s)

c





(s) P1 (c) − ℓ (c)

(P1 (c) − ℓ(c))



= ℓ (c)

k1 X

n Y (s)

k1 X

s=1

i=2

s=1

k1 X

n Y

α

(s)

α(s)

s=1

(s)



{z

AP

n Y

(s)

α(s) P1 (c) − ℓ (c)

(s)

Pi (ωi )

i=2

(s)

Pi (ωi ) + (P1 (c) − ℓ(c))

k1 X

αc(s)

s=1

i=2

|

where αc

Pi (ωi ) +

}

n Y

(s)

Pi (ωi ),

i=2

|

{z

BP,c

}

is a valid distribution over s ∈ [k1 ]. Note that P1 in the denominator  is the marginal  (t)

(t)

distribution of the mixture P over the first coordinate. Similarly, let βc ≜ Q(ω) = ℓ (c)

k2 X

β (t)

n Y

t=1

i=2

|

{z

AQ

(t)

Qi (ωi ) + (Q1 (c) − ℓ(c))

k2 X t=1

|

}

βc(t)

n Y

β (t) Q1 (c)−ℓ(c) (Q1 (c)−ℓ(c))

and

(t)

Qi (ωi ) .

i=2

{z

BQ,c

}

Each of AP , AQ , BP,c , and BQ,c is a mixture of product distributions over the remaining coordinates. The component distributions in AP and BP,c are obtained from the components of P by removing the first coordinate; the same relation holds between AQ , BQ,c and Q. However, the mixing weights may be different from those in P and Q. We use the above decomposition to construct a coupling (X, Y ) between P and Q. We first couple the first coordinate of X and Y . For each c ∈ [q], with probability ℓ(c), we set X1 = Y1 = c and then couple AP and AQ recursively for the remaining coordinates. For each c ∈ [q], with probability min{P1 (c) − ℓ(c), Q1 (c) − ℓ(c)}, we set X1 = Y1 = c and couple BP,c and BQ,c recursively 6

for the remaining coordinates. For the remaining probability mass1 of the first coordinate, there is no common part with X1 = Y1 = c for any c ∈ [q]. We arbitrarily couple this remaining mass; it must hold that X1 ̸= Y1 , and hence X = ̸ Y . We call this case the coupling failure case. Suppose that after coupling the first coordinate, we have X1 = c = ̸ c′ = Y1 for some c, c′ ∈ [q]. We then independently sample the remaining coordinates according to BP,c and BQ,c′ , respectively. The above recursive process terminates either when coupling fails at some coordinate or when it reaches the last coordinate. Denote the resulting recursive coupling by CRC . Adapting the recursive TV-distance analysis of [KECM23] to this explicit coupling, we show that CRC is not far from the optimal coupling: PrCRC [X ̸= Y ] ≤ (4nq)k1 +k2 −1 · dTV (P, Q) . Therefore, we set γ = (4nq)−k1 −k2 +1 , which is inverse-polynomial in n and q when k1 and k2 are constants.

2.3

Implementation of three operations

We now give polynomial-time algorithms for the three operations in (a), (b), and (c). We create a state space S to keep track of the recursive coupling process. A state S ∈ S is a tuple S = (i, ᾱ, β̄), where i ∈ [n] means we need to couple the remaining coordinates i, i + 1, . . . , n, and ᾱ = (ᾱ(s) )s∈[k1 ] and β̄ = (β̄ (t) )t∈[k2 ] are the current mixing weights for each mixture component. Section 2.2 describes the recursive coupling starting from the root state (1, α, β), and hence defines the quantities AP , AQ , BP,c , BQ,c and ℓ(c) in that setting. For the implementation, we use the same decomposition at an arbitrary state S = (i, ᾱ, β̄), replacing the first coordinate by the current coordinate i and the original mixing weights by the current weights ᾱ, β̄. Thus, (

ℓ(c) = min

)

min

s∈[k1 ]:ᾱ(s) >0

(s) Pi (c),

min

t∈[k2 ]:β̄ (t) >0

(t) Qi (c)

.

(s)

Similarly, the updated weights ᾱc = (ᾱc )s∈[k1 ] in the BP,c branch are defined by (s)

∀s ∈ [k1 ], ᾱc(s) = ᾱ(s) ·

X Pi (c) − ℓ(c) (s) ᾱ(s) Pi (c). where P̄i (c) = P̄i (c) − ℓ(c) s∈[k ] 1

The β̄c for the BQ,c is analogous. There is also a failure state ⊥. If we couple Xi = ̸ Yi (coupling failure case), the state S transitions to the failure state ⊥. If we couple Xi = Yi and then recursively couple AP and AQ for the remaining coordinates (we call it the good case A), the state S transitions to S ′ = (i + 1, ᾱ, β̄). If we couple Xi = Yi and then recursively couple BP,c and BQ,c for the remaining coordinates (we call it the good case B), the state S transitions to Sc′ = (i + 1, ᾱc , β̄c ). The whole recursive coupling process is a random walk on the state space S starting from the root state Sroot = (1, α, β). Every non-failure state S = (i, ᾱ, β̄) ̸=⊥ can transition to a state in {⊥} ∪ {S ′ } ∪ {Sc′ }c∈[q] , and it transitions to some Sc′ if and only if the good case B occurs. The key observation is that case B can occur at most k1 + k2 − 2 times in the whole recursive coupling (s) (t) (s) process. By the definition of ℓ(c), ℓ(c) is equal to some Pi (c) or Qi (c). Suppose ℓ(c) = Pi (c). 1

Specifically, the remaining probability for X1 = c is max{0, P1 (c) − Q1 (c)} and the remaining probability for Y1 = c is max{0, Q1 (c) − P1 (c)}.

7

(s)

(s)

(c)−ℓ(c)

Then, the new mixing weights in Sc′ satisfy ᾱc = ᾱ(s) · Pi i (c)−ℓ(c) = 0. Thus, whenever case B occurs, the total number of positive mixing weights decreases by at least one. Since the root state has k1 + k2 positive mixing weights, case B can occur at most k1 + k2 − 2 times. Using this property, we can show that the total number of states is |S| = (nq)O(k1 +k2 ) , which is polynomial in n and q. Given a state S, it is easy to compute the transition probabilities to the next states. The whole transition graph forms a DAG with polynomial size. To sample ω ∼ π in (2), it suffices to sample a trajectory starting from the root state, conditioned on the random walk reaching the failure state ⊥. We first run a dynamic programming algorithm to compute the probability of reaching ⊥ from each state in S. Then, we use these probabilities to sample the desired conditional trajectory. Similarly, using dynamic programming, we can compute f (ω) and PrCRC [X ̸= Y ] in polynomial time. P

Remark 2.1 (Comparison with the coupling in [BGM+ 24]). Finally, it is worth comparing our coupling with the coupling used in [BGM+ 24]. They proposed a simpler coupling and applied it to approximating the TV-distance between two Bayesian networks. Given two n-dimensional distributions P and Q over [q]n , they couple the coordinates one by one using the optimal coupling of the conditional marginals. Suppose the first i − 1 coordinates of X and Y have been generated and coupled consistently such that X1:i−1 = Y1:i−1 = σ ∈ [q]i−1 . For the i-th coordinates Xi and Yi , they use the optimal coupling of the conditional marginals Pi and Qi given σ. The coupling terminates if Xi ̸= Yi at some coordinate or when it reaches the last coordinate. Their coupling is different from ours. Unlike Bayesian networks, for general mixtures of product distributions, the conditional marginals at the i-th coordinate require the full information of the prefix condition σ ∈ [q]i−1 . Since there are O(q n ) possible prefix conditions σ, this coupling needs to consider exponentially many possible distributions for our problem. In our coupling, even when the good event Xi = Yi occurs, we further split it into two sub-cases, A and B. This splitting allows us to exploit the structure of mixtures of product distributions, so the total number of possible distributions appearing in our coupling is bounded by a polynomial.

3

A General Approach for Approximating the TV-Distance

When k1 = k2 = 1, the problem reduces to approximating the TV-distance between two product distributions P and Q over Ω = [q]n . This special case was studied in [FGJW23]. Although the algorithm there is tailored to product distributions, its underlying idea extends to a more general framework for approximating TV-distances for general discrete distributions. Let P and Q be arbitrary discrete distributions over Ω. A coupling of P and Q is a pair of jointly distributed random variables (X, Y ) such that X ∼ P and Y ∼ Q. For every coupling C of P and Q, the coupling inequality states that dTV (P, Q) ≤ PrC [X ̸= Y ] . A coupling C is called optimal if equality holds, namely if PrC [X ̸= Y ] = dTV (P, Q). An optimal coupling exists for every pair of distributions P and Q. The following fact holds for every coupling C of P and Q. Fact 3.1. For any coupling C of P and Q, it holds that ∀ω ∈ Ω,

PrC [X = ω ∧ X ̸= Y ] ≥ max{0, P(ω) − Q(ω)}.

If C is an optimal coupling, then for all ω ∈ Ω, PrC [X = ω ∧ X ̸= Y ] = max{0, P(ω) − Q(ω)}. 8

Proof. By definition, PrC [X = ω ∧ X ̸= Y ] = P(ω) − PrC [X = Y = ω] and PrC [X = Y = ω] ≤ min{P(ω), Q(ω)}. Therefore, PrC [X = ω ∧ X ̸= Y ] ≥ max{0, P(ω) − Q(ω)}. When C is optimal, we have dTV (P, Q) = PrC [X ̸= Y ] =

X

(∗) X

PrC [X = ω ∧ X ̸= Y ] ≥

ω∈Ω

max{0, P(ω) − Q(ω)} = dTV (P, Q) ,

ω∈Ω

which implies that the equality in (∗) must be achieved by the optimal coupling, i.e., for any ω ∈ Ω, PrC [X = ω ∧ X ̸= Y ] = max{0, P(ω) − Q(ω)}. The inequality (∗) itself is due to the fact that for any valid coupling, Pr [X = ω] = P(ω) and Pr [Y = ω] = Q(ω), and thus PrC [X = ω ∧ X ̸= Y ] = Pr [X = ω] − Pr [X = ω ∧ Y = ω] ≥ P(ω) − min{P(ω), Q(ω)} = max{0, P(ω) − Q(ω)}. We now state an abstract version of the algorithm from [FGJW23], formulated for two arbitrary discrete distributions P and Q. Let γ, T0 , T1 , T2 , T3 > 0 be parameters. Assumption 3.2. There exists a probability mass oracle that, for any ω ∈ Ω, returns P(ω) and Q(ω) in time T0 . dTV (P,Q) Assumption 3.3. There exists a coupling C between P and Q such that Pr ≥ γ, together C [X̸=Y ] with an oracle that can be preprocessed in time T1 and supports the following query:

• Discrepancy query: return the value of PrC [X ̸= Y ] in time O(1). Furthermore, if PrC [X ̸= Y ] > 0, then the oracle also supports the following two queries: • Sampling query: draw an independent sample of X from the conditional distribution of X given X = ̸ Y under C in time T2 . Equivalently, it returns ω ∈ Ω with probability PrC [X = ω | X ̸= Y ]. • Evaluation query: given any ω ∈ Ω, return the value of PrC [X = ω ∧ X ̸= Y ] in time T3 . Under Assumption 3.2 and Assumption 3.3, the following theorem gives an algorithm for approximating the TV-distance between two arbitrary discrete distributions P and Q. Theorem 3.4 ([FGJW23]). Suppose Assumption 3.2 and Assumption 3.3 hold for distributions P and Q with parameters γ, T0 , T1 , T2 , T3 . Then there exists a randomized algorithm such that, given any ε > 0, the parameters γ, and access number dˆ in h to the two oracles above, it outputs a random i T0 +T2 +T3 time O(T1 + ) satisfying Pr (1 − ε)dTV (P, Q) ≤ dˆ ≤ (1 + ε)dTV (P, Q) ≥ 99%. γε2 Proof. If PrC [X ̸= Y ] = 0, then 0 ≤ dTV (P, Q) ≤ PrC [X ̸= Y ] implies dTV (P, Q) = 0, so we may directly return dˆ = 0. We therefore assume that PrC [X ̸= Y ] > 0. By Assumption 3.3, after 9

preprocessing time T1 we can obtain the value of PrC [X ̸= Y ] in constant time. Let π denote the conditional distribution of X given the event X ̸= Y under the coupling C, namely, ∀ω ∈ Ω,

π(ω) = PrC [X = ω | X ̸= Y ] =

PrC [X = ω ∧ X ̸= Y ] . PrC [X ̸= Y ]

For all ω ∈ Ω with π(Ω) ̸= 0, it holds that PrC [X = ω ∧ X ̸= Y ] > 0, and we define f (ω) =

max{0, P(ω) − Q(ω)} . PrC [X = ω ∧ X ̸= Y ]

By Fact 3.1, 0 ≤ f (ω) ≤ 1 for all ω ∈ Ω with π(ω) > 0, hence Eπ f 2 ≤ Eπ [f ], and 

h

i

h

i



Varπ [f ] = Eπ f 2 − (Eπ [f ])2 ≤ Eπ f 2 ≤ Eπ [f ] . Now let Ω+ = {ω ∈ Ω | π(ω) > 0}. Then Eπ [f ] =

X

π(ω)

ω∈Ω+

X max{0, P(ω) − Q(ω)} max{0, P(ω) − Q(ω)} = . PrC [X = ω ∧ X ̸= Y ] ω∈Ω PrC [X ̸= Y ] +

For ω ∈ Ω with π(ω) = 0, then either PrC [X = ω] = 0 (hence P(ω) = 0), or PrC [X ̸= Y | X = ω] = 0 (which implies that whenever X = ω, we also have Y = ω, and therefore Q(ω) ≥ P(ω)). Hence, π(ω) = 0

=⇒

max{0, P(ω) − Q(ω)} = 0.

(4)

By (4), it holds that max{0, P(ω) − Q(ω)} = ̸ 0 implies π(ω) > 0. Hence, Ω+ contains all ω ∈ Ω such that max{0, P(ω) − Q(ω)} ̸= 0. Therefore, dTV (P, Q) =

X

max{0, P(ω) − Q(ω)} =

ω∈Ω

X

max{0, P(ω) − Q(ω)},

ω∈Ω+

we obtain Eπ [f ] =

dTV (P, Q) ∈ [γ, 1]. PrC [X ̸= Y ]

By the sampling query in Assumption 3.3, we can draw a sample from π in time T2 . Moreover, using Assumption 3.2 and Assumption 3.3, for any ω ∈ Ω we can compute P(ω),

Q(ω),

PrC [X = ω ∧ X ̸= Y ] ,

f (ω)

100 (1) , ω (2) , . . . ω (m) ∼ π and compute in time T0 + T3 . Let m = γε 2 . We draw m independent samples ω P f¯ = 1 m f (ω (i) ). By Chebyshev’s inequality, m

i=1

h i Varπ [f ] 1 1 Pr f¯ − Eπ [f ] ≥ εEπ [f ] ≤ 2 ≤ mε2 E [f ] ≤ 100 2 = 0.01. 2 ε γ mε Eπ [f ] π γε2

Thus, with probability at least 99%, f¯ ∈ [(1 − ε)Eπ [f ] , (1 + ε)Eπ [f ]]. The proof is completed by outputting dˆ = f¯ · PrC [X ̸= Y ], where PrC [X ̸= Y ] is the value returned by the discrepancy query in Assumption 3.3 in constant time. 10

4

Algorithm for Mixtures of Product Distributions

To apply Theorem 3.4 to the problem of approximating the TV-distance between two mixtures of product distributions, we need to verify Assumption 3.2 and Assumption 3.3. By the definition of mixture of product distributions in (1), the first assumption holds with T0 = O(n(k1 + k2 )). Our main task is to construct a coupling to verify the second assumption. Lemma 4.1. Let k1 , k2 ≥ 1 be two constants. For any two mixtures of product distributions P 1 N P 2 (t) Nn (s) (t) α(s) ni=1 Pi and Q = kt=1 β P = ks=1 i=1 Qi , there exists a coupling C between P and Q such that Assumption 3.3 holds with parameters γ = (4nq)−k1 −k2 +1 , T1 , T2 , T3 = O(q 2 (nq)k1 +k2 −1 ). Theorem 1.2 is a simple corollary of Theorem 3.4 and Lemma 4.1. The rest of this section is devoted to the proof of the lemma.

4.1

The recursive coupling for mixtures of product distributions

In this subsection, we construct a recursive coupling between two mixtures of product distributions P and Q in two steps. In Section 4.1.1, we first introduce a recursive procedure for drawing a sample from a single mixture P. Building on this procedure, in Section 4.1.2 we then construct a recursive coupling between P and Q. Our construction is inspired by the technique of [KECM23], who develop a recursive analysis that bounds the total variation distance between two mixtures of product distributions up to a fixed error. Our main contribution is a new recursive coupling process, which is crucial for approximating the TV-distance with arbitrary relative error. In Section 4.2, we use this process to implement an oracle that answers the queries in Assumption 3.3 within polynomial time. 4.1.1

The recursive sampling process

We now describe a recursive sampling procedure that generates a sample X from the mixture of product distributions P=

k1 X

α(s)

s=1

n O (s)

Pi .

i=1

This sampling procedure will later serve as the basis for our coupling construction. Throughout the recursion, the component product distributions P(s) remain fixed, while the mixing weights are 1 updated from one recursive step to the next. We use (ᾱ(s) )ks=1 to denote the mixing weights at the 1 current recursive step, to distinguish them from the original mixing weights (α(s) )ks=1 . At recursion depth j, let P̄ =

k1 X

ᾱ(s)

s=1

n O (s)

Pi

i=j

be the current mixture over [q]n−j+1 . We consider the recursive coordinate-wise procedure RecursiveSample (j, ᾱ) , where j ∈ [n] is the current coordinate. At this point, the coordinates X1 , . . . , Xj−1 have already been fixed, and the task is to generate Xj from the current mixture P̄ and then update the mixing weights for the next recursive step. To describe the procedure, we define a lower-bound selector. 11

Definition 4.2 (lower-bound selector). A lower-bound selector ℓP is a function ℓP (j, ᾱ, c) defined 1 for any coordinate j ∈ [n], any mixing-weight vector ᾱ = (ᾱ(s) )ks=1 , and any c ∈ [q], satisfying: (s)

min

s∈[k1 ]:ᾱ(s) >0

Pj (c) ≥ ℓP (j, ᾱ, c) ≥ 0.

(5)

Here the minimum ranges over all components s ∈ [k1 ] with ᾱ(s) > 0, namely over the components that are still active at the current recursive step. Thus ℓP (j, ᾱ, c) is required to be a nonnegative lower bound on the marginal probability of seeing value c in coordinate j across all active components. The above definition does not specify a unique lower-bound selector ℓP . Any function satisfying (5) is valid. Let us fix an arbitrary lower-bound selector ℓP satisfying (5). Our recursive sampling process in Section 4.1.1 works for any lower-bound selector ℓP . In Section 4.1.2, when constructing the coupling, we will specify a particular choice of ℓP . 1 For the current recursive state (j, (ᾱ(s) )ks=1 ), we first compute the lower bounds ℓP (j, ᾱ, c) for all c ∈ [q]. Fix a value c ∈ [q] and an assignment ω ∈ [q]n−j+1 satisfying ωj = c. The probability of ω under the current mixture P̄ can be written as P̄(ω) =

k1 X

ᾱ

n Y (s)

k1 X

i=j

s=1

Pi (ωi ) =

(s)

s=1

(s)

ᾱ(s) Pj (c)

n Y

(s)

Pi (ωi ),

i=j+1 (s)

where the second equality uses the assumption ωj = c. By (5), it holds that Pj (c) ≥ ℓP (j, ᾱ, c) (s)

for all s ∈ [k1 ] with ᾱ(s) > 0. Hence, for each s with ᾱ(s) > 0, we can decompose Pj (c) as two (s)

non-negative terms: ℓP (j, ᾱ, c) and Pj (c) − ℓP (j, ᾱ, c). We have the following decomposition: P̄(ω) = ℓP (j, ᾱ, c)

k1 X

ᾱ(s)

n Y

s=1

i=j+1

|

{z

(s)

Pi (ωi ) +

k1 X



n  Y

(s)

ᾱ(s) Pj (c) − ℓP (j, ᾱ, c)

s=1

(s)

Pi (ωi ).

(6)

i=j+1

}

A

The above summation enumerates all s ∈ [k1 ], which may contain inactive components s with ᾱ(s) = 0. However, those inactive components only contribute 0 to the whole sum. In (6), the first term is the product of ℓP (j, ᾱ, c) and A, where A is exactly the probability of the suffix (ωi )ni=j+1 (s) (s) Nn i=j+1 Pi . To write the second term in a similar form, we (s) introduce a modified mixing weight ᾱj,c for all s ∈ [k1 ] as follows:

under the current mixture

Pk1

 ᾱ(s) ,

s=1 ᾱ

(s) (s) P ᾱj,c ≜ ᾱ(s) Pj (c)−ℓ (j,ᾱ,c) , P P̄j (c)−ℓ (j,ᾱ,c)

P̄j (c) − ℓP (j, ᾱ, c) = 0; otherwise,

where P̄j (c) =

k1 X

(s)

ᾱ(s) Pj (c).

s=1 (s)

We next verify that ᾱj,c defines a valid mixing-weight vector. By Definition 4.2, Pj (c)−ℓP (j, ᾱ, c) ≥ (s)

(s)

1 0, and P̄j (c) = ks=1 ᾱ(s) Pj (c) is an average of the active values Pj (c), therefore P̄j (c) ≥ ℓP (j, ᾱ, c). If P̄j (c) − ℓP (j, ᾱ, c) = 0, ᾱj,c = ᾱ is a valid mixing-weight vector as long as ᾱ is valid. Otherwise, when P̄j (c) > ℓP (j, ᾱ, c),

P

k1 X (s)

k1 X

s=1

s=1

ᾱj,c =

(s)

ᾱ

(s)

Pj (c) − ℓP (j, ᾱ, c) P̄j (c) − ℓP (j, ᾱ, c)

(s) P(s) (c) ℓP (j, ᾱ, c) j = − = 1. P̄j (c) − ℓP (j, ᾱ, c) P̄j (c) − ℓP (j, ᾱ, c)

Pk1

s=1 ᾱ

12

Hence, for ω ∈ [q]n−j+1 with ωj = c, we can rewrite P̄(ω) as P̄(ω) = ℓ (j, ᾱ, c) P

k1 X s=1

ᾱ

(s)

n Y

(s) Pi (ωi ) +



k1 X

P̄j (c) − ℓ (j, ᾱ, c) P

(s) ᾱj,c

s=1

i=j+1

n Y

(s)

Pi (ωi ).

(7)

i=j+1

1 Equation (7) naturally suggests a recursive sampling rule. Given the state (j, (ᾱ(s) )ks=1 ), we first sample the current coordinate Xj , and then reduce the task of sampling the suffix Xj+1:n to sampling from one of two smaller mixtures. Concretely, once the current coordinate is set to P 1 N P 1 (s) (s) N (s) c, the suffix distribution is either ks=1 ᾱ(s) ni=j+1 Pi or ks=1 ᾱj,c ni=j+1 Pi . The component distributions remain the same (except that their dimension is reduced by one), while the mixing weights may change from ᾱ to ᾱj,c . Applying this rule recursively yields the sampling procedure shown in Algorithm 1.

Algorithm 1: Recursive sampling process for mixtures of product distributions 1 RecursiveSample (j, ᾱ);

Parameter : an index j ∈ [n + 1] and mixing weights ᾱ(1) , ᾱ(2) , . . . , ᾱ(k1 ) Static Input : product distributions P(1) , P(2) , . . . , P(k1 ) over [q]n ; an arbitrary lower-bound selector function ℓP satisfying (5) P 1 N (s) Output : a random sample X = (Xj , Xj+1 , . . . , Xn ) ∼ P̄ = ks=1 ᾱ(s) ni=j Pi 2 if j > n then return ∅; 3 Compute the probability lower bound ℓP (j, ᾱ, c) for all c ∈ [q]; 4 Sample a random pair (c, e) ∈ [q] × {0, 1} such that for each c ∈ [q], the probability of (c, 0) is ℓP (j, ᾱ, c) and the probability of (c, 1) is P̄j (c) − ℓP (j, ᾱ, c), where P 1 (s) P̄j (c) = ks=1 ᾱ(s) Pj (c) is the marginal probability of coordinate j under P̄; 5 if e = 0 then 6 Sample (Xj+1 , Xj+2 , . . . , Xn ) ← RecursiveSample(j + 1, ᾱ); 7 else 8

(s)

(s)

P

(c)−ℓP (j,ᾱ,c)

Compute mixing weights ᾱj,c = ᾱ(s) P̄j (c)−ℓP (j,ᾱ,c) for all s ∈ [k1 ]; j

9

Sample (Xj+1 , Xj+2 , . . . , Xn ) ← RecursiveSample(j + 1, ᾱj,c );

10 return (Xj = c, Xj+1 , Xj+2 , . . . , Xn );

In Algorithm 1, the lower-bound selector function ℓP is treated as part of the input. The following lemma shows that the algorithm outputs a correct sample from P as long as ℓP satisfies (5). (s)

1 Lemma 4.3. Given a mixture of product distributions P = ks=1 α(s) ni=1 Pi and an arbitrary lower-bound selector function ℓP satisfying (5), the procedure RecursiveSample(1, α) in Algorithm 1 outputs a sample X ∼ P.

P

N

Lemma 4.3 follows by induction on the recursion depth using (7). The proof is in Appendix A. The same recursive procedure applies to any mixture of product distributions. In particular, we P 2 (t) Nn (t) (s) k1 may use it to sample from Q̄ = kt=1 β̄ i=j Qi by replacing k1 with k2 , replacing (ᾱ )s=1 with k 2 (β̄ (t) )t=1 , replacing P(s) with Q(t) , and using an arbitrary lower-bound selector ℓQ (j, β̄, c) satisfying min

t∈[k2 ]:β̄ (t) >0

(t)

Qj (c) ≥ ℓQ (j, β̄, c) ≥ 0. 13

(8)

Lemma 4.3 also shows that by substituting P with Q and α with β, the procedure RecursiveSample(1, β) outputs a sample Y ∼ Q. 4.1.2

The recursive coupling process

We now define a coupling between two mixtures of product distributions P=

k1 X

α

s=1

(s)

n O (s)

and

Pi

Q=

k2 X

β

(t)

t=1

i=1

n O

(t)

Qi .

i=1

At a high level, the coupling runs Algorithm 1 on both P and Q to generate samples X ∼ P and Y ∼ Q, while sharing randomness between the two executions so that the resulting pair (X, Y ) is a coupling of P and Q. We use RecursiveSample(1, α) to generate a sample X ∼ P and RecursiveSample(1, β) to generate a sample Y ∼ Q, and couple the two recursive procedures at the same time. Suppose that at some point the execution for P is at recursion state RecursiveSample(j, ᾱ) and the execution 1 2 for Q is at recursion state RecursiveSample(j, β̄), where ᾱ = (ᾱ(s) )ks=1 and β̄ = (β̄ (t) )kt=1 denote the current mixing weights and may differ from the original vectors α and β. We define a shared lower-bound selector ℓ(j, ᾱ, β̄, c) as follows. 1 Definition 4.4. For any j ∈ [n], c ∈ [q], and valid mixing-weight vectors ᾱ = (ᾱ(s) )ks=1 and k (t) 2 β̄ = (β̄ )t=1 , define





(

ℓ j, ᾱ, β̄, c ≜ min

)

min

s∈[k1 ]:ᾱ(s) >0

(s) Pj (c),

min

t∈[k2 ]:β̄ (t) >0

(t) Qj (c)

.

The one-sided lower-bound selectors used by the two recursive sampling procedures are then defined by ℓP (j, ᾱ, c) = ℓQ (j, β̄, c) = ℓ(j, ᾱ, β̄, c). By construction, ℓP (j, ᾱ, c) satisfies (5) for P, and ℓQ (j, β̄, c) satisfies (8) for Q. In Line 4 of Algorithm 1, each process draws a random pair (c, e) ∈ [q] × {0, 1}. Let (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ) denote the pairs drawn by RecursiveSample(j, ᾱ) and RecursiveSample(j, β̄), respectively. Since both procedures use the same lower bounds ℓ(j, ᾱ, β̄, c), we define a coupling C¯ between these two pairs as follows: h

i

h

i

• For each c ∈ [q], PrC¯ cP̄ = cQ̄ = c ∧ eP̄ = eQ̄ = 0 = ℓ(j, ᾱ, β̄, c). • For each c ∈ [q], PrC¯ cP̄ = cQ̄ = c ∧ eP̄ = eQ̄ = 1 = min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c), where (s)

1 ᾱ(s) ni=j Pi and Q̄ = P̄ = ks=1 at the current recursive step.

P

N

Pk2

t=1 β̄

(t) Nn Q(t) are the mixtures of product distributions i=j i

• Given the above two points, for each c ∈ A = {x ∈ [q] | P̄j (x) > Q̄j (x)}, the random variable (cP̄ , eP̄ ) has a remaining probability P̄j (c) − Q̄j (c) of taking value (c, 1); and for each c ∈ B = {x ∈ [q] | P̄j (x) < Q̄j (x)}, the random variable (cQ̄ , eQ̄ ) has remaining probability Q̄j (c) − P̄j (c) of taking value (c, 1). The total remaining probabilities on the two sides are P P equal: c∈A (P̄j (c) − Q̄j (c)) = c∈B (Q̄j (c) − P̄j (c)). We therefore couple the remaining mass arbitrarily. Since A and B are disjoint, in this case we necessarily have cP̄ ̸= cQ̄

and 14

eP̄ = eQ̄ = 1.

Using the coupling C¯ between (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ), we now couple two sampling processes RecursiveSample(j, ᾱ) and RecursiveSample(j, β̄) as follows: ¯ • Sample (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ) according to C. • If cP̄ = cQ̄ = c for some c ∈ [q] and eP̄ = eQ̄ = 0, then we set Xj = Yj = c and recursively couple the suffixes via RecursiveCouple(j + 1, ᾱ, β̄). • If cP̄ = cQ̄ = c for some c ∈ [q] and eP̄ = eQ̄ = 1, then we set Xj = Yj = c and recursively couple the suffixes via RecursiveCouple(j + 1, ᾱj,c , β̄j,c ). • In the remaining case, we have Xj = cP̄ ̸= cQ̄ = Yj and eP̄ = eQ̄ = 1. Hence the coupling already fails at coordinate j, and we simply sample the suffixes independently using RecursiveSample(j + 1, ᾱj,cP̄ ) and RecursiveSample(j + 1, β̄j,cQ̄ ). The pseudocode of the coupling process is given in Algorithm 2. Algorithm 2: Recursive coupling process for mixtures of product distributions 1 RecursiveCouple (j, ᾱ, β̄); 1 2 Parameter : an index j ∈ [n + 1] and mixing weights ᾱ = (ᾱ(s) )ks=1 , β̄ = (β̄ (t) )kt=1 (1) (2) (k ) (1) (2) (k ) Static Input : product distributions P , P , . . . , P 1 and Q , Q , . . . , Q 2 over [q]n ; a fixed shared probability lower-bound selector ℓ as in Definition 4.4 Output : joint random variables (X, Y ) such that P 1 N (s) X = (Xj , Xj+1 , . . . , Xn ) ∼ P̄ = ks=1 ᾱ(s) ni=j Pi and

(t)

2 Y = (Yj , Yj+1 , . . . , Yn ) ∼ Q̄ = kt=1 β̄ (t) ni=j Qi 2 if j > n then return (∅, ∅); 3 Compute ℓ(j, ᾱ, β̄, c) for all c ∈ [q], and then use the coupling C¯ to jointly sample (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ); 4 if eP̄ = eQ̄ = 0 then 5 Sample (Xj+1:n , Yj+1:n ) ← RecursiveCouple(j + 1, ᾱ, β̄);

P

N

6 else 7

(s)

P̄

8 9 10 11 12

P

(s)

(c )−ℓ(j,ᾱ,β̄,c )

(t)

(t)

Q

(c )−ℓ(j,ᾱ,β̄,c )

Let ᾱj,c = ᾱ(s) P̄j (c P̄)−ℓ(j,ᾱ,β̄,c P̄) for s ∈ [k1 ] and β̄j,c = β̄ (t) Q̄j (c Q̄)−ℓ(j,ᾱ,β̄,c Q̄) for t ∈ [k2 ]; j

P̄

Q̄

P̄

j

Q̄

Q̄

if cP̄ = cQ̄ then Sample (Xj+1:n , Yj+1:n ) ← RecursiveCouple(j + 1, ᾱj,c , β̄j,c ), where c = cP̄ = cQ̄ ; else Sample Xj+1:n ← RecursiveSample(j + 1, ᾱj,cP̄ ) with ℓP = 0; Sample Yj+1:n ← RecursiveSample(j + 1, β̄j,cQ̄ ) with ℓQ = 0;

13 Let X = (Xj = cP̄ , Xj+1:n ) and Y = (Yj = cQ̄ , Yj+1:n ) and return (X, Y );

Lemma 4.5 (Validity of recursive coupling). Given two mixtures of product distributions P = Pk1 (s) Nn P(s) and Q = Pk2 β (t) Nn Q(t) , fix a shared lower-bound selector ℓ as in Definii=1 i i=1 i s=1 α t=1 tion 4.4. Then RecursiveCouple(1, α, β) in Algorithm 2 outputs joint random variables (X, Y ) such that X ∼ P and Y ∼ Q. 15

Lemma 4.5 follows by a straightforward induction. We defer the proof to Appendix A. 4.1.3

Discrepancy of recursive coupling

We now analyze the discrepancy of the recursive coupling process. Let P = P 2 (t) Nn (t) Q = kt=1 β i=1 Qi be two mixtures of product distributions. Define

Pk1

s=1 α

(s) Nn P(s) and i=1 i

(X, Y ) ← RecursiveCouple(1, α, β). Let CRC denote the joint distribution of (X, Y ). By the coupling inequality, the discrepancy PrCRC [X ̸= Y ] is always at least dTV (P, Q). The following lemma provides an upper bound on the discrepancy of our recursive coupling. Lemma 4.6 (Discrepancy of recursive coupling). Let [q] be a finite domain. For two mixtures of P 1 N P 2 (t) Nn (s) (t) n product distributions P = ks=1 α(s) ni=1 Pi and Q = kt=1 β i=1 Qi over [q] , the recursive coupling satisfies PrCRC [X ̸= Y ] ≤ (4nq)k1 +k2 −1 · dTV (P, Q) . The proof of Lemma 4.6 follows the known techniques of [KECM23]. We include the proof in Appendix B for completeness.

4.2

Implementation of oracles

In this section, we implement the oracle queries in Assumption 3.3 based on Algorithm 2. Consider P 1 N (s) the execution of Algorithm 2 on two mixtures of product distributions P = ks=1 α(s) ni=1 Pi and P 2 (t) Nn (t) Q = kt=1 β i=1 Qi . Define the following set of states and transitions between them. For each recursion instance RecursiveCouple(j, ᾱ, β̄), we create a state S with label (j, ᾱ, β̄), where 1 ≤ j ≤ n + 1. A state S is said to be in layer j if j is the first coordinate of the label. We add a special failure state S =⊥ to denote that the coupling fails. Let S denote the set of all possible states that can be reached by the coupling process with a positive probability, including the failure state ⊥. For each state S ∈ S with label (j, ᾱ, β̄) for j ∈ [n], we define the following transitions according to the coupling process. The coupling process samples (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ) for P̄ and Q̄, respectively, and then makes the transition. • Type-I: For any c ∈ [q] with ℓ(j, ᾱ, β̄, c) > 0, with probability ℓ(j, ᾱ, β̄, c), S moves to a new state S ′ with label (j +1, ᾱ, β̄). Define a transition t from S to S ′ with weight w(t) = ℓ(j, ᾱ, β̄, c) and label L(t) = (c, c). This means that S moves to S ′ through t with probability ℓ(j, ᾱ, β̄, c) and sets Xj = Yj = c. • Type-II: For any c ∈ [q] with pc ≜ min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) > 0, with probability pc , S moves to a new state S ′ with label (j + 1, ᾱj,c , β̄j,c ). Define a transition t from S to S ′ with weight w(t) = pc and label L(t) = (c, c). This means that S moves to S ′ through t with probability pc and sets Xj = Yj = c. • Type-III: With the remaining probability,i the coupling fails (Xj ̸= Yj ). For any c, c′ ∈ [q] h with c ̸= c′ and pc,c′ ≜ Pr cP̄ = c ∧ cQ̄ = c′ > 0, define a transition t from S to the failure 16

state ⊥ with weight w(t) = pc,c′ and label L(t) = (Xj:n , Yj:n ), where Xj = c, Yj = c′ , and Xj+1:n , Yj+1:n are sampled independently in Lines 11 and 12 of Algorithm 2. Hence, the label contains a pair of random configurations, both of length n − j + 1. The label L(t) is random because we further sample the pair (Xj+1:n , Yj+1:n ). We remark that there are multiple transitions from S to the failure state ⊥ with different weights and labels. Lemma 4.7. The total number of states is bounded by |S| ≤ (nq + 1)k1 +k2 −1 + 1. Proof. For any non-failed state S (i.e., S ̸=⊥), recall that it is labeled by (j, ᾱ, β̄). Define the active component count m(S) ≜ {s ∈ [k1 ] : ᾱ(s) > 0} + {t ∈ [k2 ] : β̄ (t) > 0} . At the root Sroot = (1, α, β), it holds that m(Sroot ) ≤ k1 + k2 . Consider any transition from a non-failed state S with label (j, ᾱ, β̄) in layer j. Consider a type-I transition t from S to some state SI , where the label on t is L(t) = (c, c). Then SI keeps the same mixing weights, hence m(SI ) = m(S). Consider a type-II transition from S to some state SII with label L(t) = (c, c). Note that in this case, we must have min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) > 0. The mixing weights in SII are updated to (s) Pj (c) − ℓ(j, ᾱ, β̄, c) (s) (s) ᾱj,c = ᾱ · P̄j (c) − ℓ(j, ᾱ, β̄, c) (t) Qj (c) − ℓ(j, ᾱ, β̄, c) (t) (t) β̄j,c = β̄ · Q̄j (c) − ℓ(j, ᾱ, β̄, c)

for all s ∈ [k1 ], for all t ∈ [k2 ].

By the definition of ℓ(j, ᾱ, β̄, c) in Definition 4.4, there exists at least one active index s ∈ [k1 ] (s) (t) or t ∈ [k2 ] such that (Pj (c) − ℓ(j, ᾱ, β̄, c) = 0 and ᾱ(s) > 0) or (Qj (c) − ℓ(j, ᾱ, β̄, c) = 0 and β̄ (t) > 0). Note that the two denominators P̄j (c) − ℓ(j, ᾱ, β̄, c) and Q̄j (c) − ℓ(j, ᾱ, β̄, c) are both positive. Therefore, whenever the process moves to SII , at least one previously-active mixture component becomes inactive, and m(SII ) ≤ m(S) − 1.

(9)

Finally, for a type-III transition, the next state is SIII =⊥, and the whole process terminates. All transitions except type-III transitions move to a new state in the next layer j + 1. We bound the number of states by an encoding argument. Now fix a layer j ∈ {1, 2, . . . , n + 1}. Any non-failed state at layer j is encoded by the choices made along a path from the root to that state. Along this path, at each previous layer we either take a type-I transition (which we encode by the symbol 0) or take a type-II transition with label (c, c) for some value c ∈ [q]. Hence, given the root state Sroot = (1, α, β), we can use a string d of length j − 1 with characters in [q] ∪ {0} to encode the label of the current state. To decode such a string, start from Sroot at layer 1. For each h = 1, 2, . . . , j − 1, if dh = 0, update the current state (h, ᾱ, β̄) to (h + 1, ᾱ, β̄); if dh = c for some c ∈ [q], update it to (h + 1, ᾱh,c , β̄h,c ). Let r be the number of nonzero symbols in the code d. Using (9), m(Sroot ) ≤ k1 + k2 , and the fact that the two mixing-weight vectors remain valid probability distributions throughout the recursive coupling process (so each vector has at least one active component), we have r ≤ m(Sroot ) − 2 = k1 + k2 − 2. 17

We count the number of possible strings d of length j − 1 that encode a non-failed state at layer j−1
1 ≤ j ≤ n + 1. For a fixed j and r, there are at most r ways to choose the r positions at which a non-zero symbol occurs, and then q r ways to choose the corresponding nonzero symbols. Hence the number of non-failed states over all layers is at most n+1 X k1 +k X2 −2 j=1

r=0

k +k −2

1X 2 j−1 r q ≤ (n + 1) (nq)r ≤ (nq + 1)k1 +k2 −1 . r r=0

!

Finally, there is one failure state. Hence the total number of states is at most (nq + 1)k1 +k2 −1 + 1. The oracle queries in Assumption 3.3 can be implemented by the following algorithm. One execution of the coupling process can be represented as a random walk on the state space S. The random walk starts from the root state S = Sroot = (1, α, β). At each step, the random walk samples a transition t = (S → S ′ ) from S with probability w(t) and moves to the next state S ′ . The random walk ends either at a state in layer n + 1 (coupling succeeds) or at the failure state ⊥ (coupling fails). The coupling process outputs a pair (X1:n , Y1:n ), which is determined by the labels of all transitions along the random walk. Preprocessing step We build the DAG with all states in S and all transitions between them. Using Lemma 4.7, the total number of states and transitions in the DAG is at most O(q 2 (nq + 1)k1 +k2 −1 ). For each state S ∈ S, compute pfail (S), the probability that a random walk starting from S reaches the failure state ⊥. By definition, pfail (⊥) = 1 and pfail (S) = 0 for all S in layer n + 1. For each non-failed state S ∈ S in layer j ∈ [n], we have the following recursive relation: pfail (S) =

X

w(t)pfail (S ′ ).

t=(S→S ′ )∈transitions from S

Using dynamic programming, we can compute pfail (S) for all S ∈ S in time O(q 2 (nq + 1)k1 +k2 −1 ). By the definition of pfail (S), the failure probability of the recursive coupling is PrCRC [X ̸= Y ] = pfail (Sroot ). Sampling query Next, we give an algorithm that generates a random sample of X conditional on X = ̸ Y in the recursive coupling. Equivalently, we can sample a random walk on the DAG starting from the root state Sroot , conditional on the random walk reaching the failure state ⊥. If pfail (Sroot ) = 0, then dTV (P, Q) = 0 by the coupling inequality, and this sampling query is not needed. Otherwise, the query can be implemented by the following algorithm. • Initialize S = Sroot . ′

fail (S ) • Sample a random transition t = (S → S ′ ) from S with probability w(t)p and then update pfail (S) ′ the current state S to S .

• The process terminates when S reaches the failure state ⊥. The above process gives a random sequence of transitions on the DAG. By concatenating the labels in the sequence, we obtain a random pair (X, Y ) ∈ [q]n × [q]n , and we output the sample X ∈ [q]n . The correctness of the conditional sample X follows directly from the definition of pfail (·). The running time to generate a random sample X is O(q 2 (nq + 1)k1 +k2 −1 ). 18

Evaluation query

Given a fixed configuration σ ∈ [q]n , we want to evaluate the probability PrCRC [X = σ ∧ X ̸= Y ] .

For each non-failed state S ∈ S in layer j ∈ [n], let pσfail (S) be the probability that a random walk starting from S reaches the failure state ⊥ and that the label Xj:n of the random walk is exactly σj:n . By definition, since Sroot is the layer-1 state, it holds that pσfail (Sroot ) = PrCRC [X = σ ∧ X ̸= Y ] . For each S ∈ S in layer j = n + 1, by definition, we have pσfail (S) = 0. Fix any non-failed state S ∈ S in layer j ∈ [n] with S = (j, ᾱ, β̄). Define T = {t | t is a transition from S to S ′ ̸=⊥ and L(t) = (σj , σj )}. The above definition considers all transitions from S to non-failed states S ′ . We still need to consider the type-III transitions from S to the failure state ⊥. By definition, there are at most q(q − 1) such transitions. Let T⊥ (S) denote the set of all type-III transitions from S to the failure state ⊥. Every such transition can be specified by a pair (c, c′ ) with c = ̸ hc′ . We use tc,c′ to i denote the corresponding ′ ′ type-III transition. Note that the weight w(tc,c ) is Pr cP̄ = c ∧ cQ̄ = c in the coupling at Line 3 of Algorithm 2. The label of tc,c′ is a random pair (Xj:n , Yj:n ), where Xj = c, Yj = c′ , and Xj+1:n and Yj+1:n are sampled independently in Lines 11 and 12 of Algorithm 2. We define τ (tc,c′ , σ) ≜ Pr [Xj+1:n = σj+1:n ] =

k1 n X (s) Y

ᾱj,c

s=1

(s)

Pi (σi ).

i=j+1

Now, we have the following recursive relation: pσfail (S) =

w(t)pσfail (S ′ ) +

X t=(S→S ′ )∈T

X

w(tc,c′ )τ (tc,c′ , σ).

tc,c′ ∈T⊥ (S):c=σj

Again, we can use dynamic programming on the DAG from layer n + 1 to layer 1 to compute all values pσfail (S) for all S ∈ S in time O(q 2 (nq + 1)k1 +k2 −1 ). The final result is pσfail (Sroot ).

5

Algorithm for Mixtures of Boolean Subcubes

5.1

Reduction to a counting problem (s)

(t)

1 2 Consider two mixtures of Boolean subcubes P = ks=1 α(s) ni=1 Pi and Q = kt=1 β (t) ni=1 Qi . (s) A component P(s) over {0, 1}n is a Boolean subcube if, for any i ∈ [n], Pi (1) ∈ {0, 12 , 1}. Thus, for

P

N

P

N

(s)

any ω ∈ {0, 1}n , P(s) (ω) = ni=1 Pi (ωi ) is a product of values in {0, 12 , 1}. We define the following partition Λone , Λzero , and Λhalf of [n]: P(s) P(s) P(s) Q

(s) (s) zero half Λone P(s) = {i ∈ [n] | Pi (1) = 1}, ΛP(s) = {i ∈ [n] | Pi (1) = 0}, ΛP(s) =



1 (s) i ∈ [n] | Pi (1) = 2



.

In words, for the product distribution P(s) , the set Λone contains all dimensions whose value is fixed P(s) zero to 1, the set ΛP(s) contains all dimensions whose value is fixed to 0, and the set Λhalf contains all P(s) 19

dimensions whose value is sampled uniformly at random. Similarly, for a component Q(t) , we define the sets Λone , Λzero , and Λhalf in the same way. Q(t) Q(t) Q(t) The probability of each ω ∈ {0, 1}n can be written as follows: P(s) (ω) = Q(t) (ω) =

  Λhalf

1 2

zero · 1 ∀i ∈ Λone P(s) , ωi = 1 · 1 ∀i ∈ ΛP(s) , ωi = 0

P(s)

h

  Λhalf

1 2

i

h

i

zero · 1 ∀i ∈ Λone Q(t) , ωi = 1 · 1 ∀i ∈ ΛQ(t) , ωi = 0 ,

Q(t)

h

i

h

i

where the indicator function 1[∗] returns 1 if statement ∗ holds and 0 otherwise. For the distribution P(s) , we say ω ∈ {0, 1}n is feasible if for all i ∈ Λone , ωi = 1 and for all i ∈ Λzero , ωi = 0. Each P(s) P(s) |Λhalf |

feasible ω ∈ {0, 1}n has probability (1/2) P(s) under the distribution P(s) , while each infeasible ω ∈ {0, 1}n has probability 0. A similar structural property holds for the distribution Q(t) . The feasibility of a configuration ω ∈ {0, 1}n under the distribution P(s) can be expressed by a Boolean formula FP(s) : {0, 1}n → {0, 1} as follows: FP(s) (ω) =

^

^

ωi ∧

i∈Λone P(s)

(10)

¬ωi .

i∈Λzero P(s)

Similarly, for the distribution Q(t) , define the Boolean formula FQ(t) : {0, 1}n → {0, 1} as follows: FQ(t) (ω) = |Λhalf | (s)

By definition, P(s) (ω) = (1/2)

P

^

^

ωi ∧

i∈Λone Q(t)

(11)

¬ωi .

i∈Λzero Q(t) |Λhalf | (t)

· FP(s) (ω) and Q(t) (ω) = (1/2)

· FQ(t) (ω). Nn (s) (s) i=1 Pi s=1 α Q

Pk1

variation distance between mixtures of Boolean subcubes P = Pk2 (t) Nn (t) i=1 Qi can be expressed as follows: t=1 β k

dTV (P, Q) =

The total and Q =

k

n n 2 1 X Y Y 1 X 1 X X (t) (s) Qi (ωi ) Pi (ωi ) − β (t) |P(ω) − Q(ω)| = α(s) 2 ω∈{0,1}n 2 ω∈{0,1}n s=1 t=1 i=1 i=1 k

1 1 X X = α(s) · 2 ω∈{0,1}n s=1

  Λhalf (s) 1 P

2

· FP(s) (ω) −

k2 X

β

(t)

·

  Λhalf (t) 1 Q

2

t=1

· FQ(t) (ω) .

Define the characteristic function F : {0, 1}n → {0, 1}k1 +k2 as follows: 



F (ω) = FP(1) (ω), FP(2) (ω), . . . , FP(k1 ) (ω), FQ(1) (ω), FQ(2) (ω), . . . , FQ(k2 ) (ω) . By the above calculation, the value of |P(ω)−Q(ω)| is fully determined by F (ω). The key observation is that the image space of F has constant size 2k1 +k2 , so we can reorganize all ω ∈ {0, 1}n according to the value of the characteristic function F (ω). Formally, 1 dTV (P, Q) = 2

X χ∈{0,1}k1 +k2

Nχ ·

k1 X s=1

α

(s)

·

  Λhalf (s) 1 P

2 20

· χs −

k2 X t=1

β

(t)

·

  Λhalf (t) 1 Q

2

· χt+k1 ,

(12)

where Nχ counts the number of ω ∈ {0, 1}n such that F (ω) = χ, i.e., Nχ = |{ω ∈ {0, 1}n | F (ω) = χ}| . By (12), computing the total variation distance between two mixtures of Boolean subcubes P and Q reduces to counting Nχ for all χ ∈ {0, 1}k1 +k2 .

5.2

The counting algorithm

Lemma 5.1. Given any χ ∈ {0, 1}k1 +k2 , the value of Nχ can be computed in time O(2|S0 | ·n(k1 +k2 )), where S0 = { j ∈ [k1 + k2 ] | χj = 0 } is the set of indices j with χj = 0. Fix χ ∈ {0, 1}k1 +k2 , and consider the set of ω ∈ {0, 1}n such that F (ω) = χ. We need to count the number of such ω ∈ {0, 1}n . To simplify the notation, define a sequence of Boolean functions F1 , F2 , . . . , Fk1 +k2 : {0, 1}n → {0, 1} as follows: FP(j) (ω) if j ∈ [k1 ], FQ(j−k1 ) (ω) if j ∈ [k1 + 1, k1 + k2 ],

(

∀j ∈ [k1 + k2 ],

Fj (ω) =

where FP(s) and FQ(t) are defined in (10) and (11), respectively. For a configuration ω ∈ {0, 1}n with F (ω) = χ, we have Fj (ω) = χj for all j ∈ [k1 + k2 ]. For a subset S ⊆ [k1 + k2 ] of indices, define the following set of configurations: Φ(S) = {ω ∈ {0, 1}n | Fj (ω) = 1 for all j ∈ S}. Let S1 = {j ∈ [k1 + k2 ] | χj = 1} and S0 = {j ∈ [k1 + k2 ] | χj = 0}. By the definition of Nχ and the inclusion-exclusion principle, we have Nχ = |Φ(S1 )| − |{ω ∈ Φ(S1 ) | ∃ j ∈ S0 , Fj (ω) = 1}| =

X

(−1)|S| · |Φ(S1 ∪ S)|.

(13)

S⊆S0

Observation 5.2. Given any S ⊆ [k1 + k2 ], the size |Φ(S)| can be computed in time O(n · |S|). Proof. By the definition of FP(s) and FQ(t) , we have FP(s) (ω) = 1 iff ωi = 1 for all i ∈ Λone and P(s) one and ω = 0 for all i ∈ Λzero . ωi = 0 for all i ∈ Λzero . Similarly, F (ω) = 1 iff ω = 1 for all i ∈ Λ (t) i i Q P(s) Q(t) Q(t) We can go through all j ∈ S and fix the values of the corresponding dimensions i ∈ [n]. If there is a contradiction, then |Φ(S)| = 0. Otherwise, |Φ(S)| = 2r , where r is the number of dimensions whose values are not fixed. The running time is O(n · |S|). We now prove Lemma 5.1. Combining (13) and Observation 5.2, we can exactly compute the value of Nχ in time !

X S⊆S0

O(n|S1 ∪ S|) ≤

X 0≤k≤|S0 |

|S0 | · O(n(k1 + k2 )) = O(2|S0 | · n(k1 + k2 )). k

21

5.3

The total variation distance algorithm

The algorithm for computing the total variation distance between two mixtures of Boolean subcubes P and Q is obtained by combining (12) and Lemma 5.1. Let K = k1 + k2 . The overall running time is !

X 0≤m≤K

K · O(2m · nK) = O(3K K · n) = O(3k1 +k2 (k1 + k2 ) · n). m

Here m enumerates the number of indices j with χj = 0. This proves the main result in Theorem 1.4.

6

Proof of the hardness result

Proof of Theorem 1.5. We show that if there is a polynomial-time algorithm for computing the TV-distance between two mixtures of Boolean subcubes with k1 + k2 = Θ(n), then there is a polynomial-time algorithm for #3SAT, restricted to 3-CNF formulas. Let φ(x1 , x2 , . . . , xr ) = C1 ∧ C2 ∧ . . . ∧ Cm be a 3-CNF formula over the used variables x1 , x2 , . . . , xr , and m ≥ 1. Set N = max{r, m}. If N > r, we add N − r dummy variables xr+1 , . . . , xN . Let S denote the number of satisfying assignments of φ over x1 , . . . , xN . As the dummy variables do not appear in any clause, the number of satisfying assignments for φ is given by #SAT(φ) =

S 2N −r

.

We establish the proof by showing the relation between S and the TV-distance between certain mixtures of Boolean subcubes. Let n = N + 1. We construct two mixtures of subcubes P and Q over {0, 1}n . A point in this space can be represented as (b, ω1 , . . . , ωN ), where b ∈ {0, 1} is an additional bit, and (ω1 , . . . , ωN ) is an assignment to (x1 , . . . , xN ). For each j ∈ [m], let xj,a , xj,b , xj,c be the three variables appearing in the clause Cj . Define (ωj,a , ωj,b , ωj,c ) ∈ {0, 1}3 to be the unique assignment to xj,a , xj,b , xj,c that falsifies Cj . Let Uj be a subcube such that PrUj [(b, xj,a , xj,b , xj,c ) = (0, ωj,a , ωj,b , ωj,c )] = 1, and the other N − 3 coordinates are uniform over {0, 1}. Define P=

m 1 X Uj , m j=1

is a mixture of k1 = m subcubes. Next we define Q. Let V0 and V1 be two subcubes such that PrV0 [b = 0] = PrV1 [b = 1] = 1, 1 and that are uniform over (x1 , . . . , xN ). Let λ = 2m , and define

Q = λV0 + (1 − λ)V1 , 22

which is a mixture of k2 = 2 subcubes. We now compute dTV (P, Q). For any assignment ω = (ω1 , . . . , ωN ), let F (ω) = |{j ∈ [m] | ω falsifies Cj }| denote the number of clauses falsified by ω. For every ω ∈ {0, 1}N , we have P(0, ω) =

m m 1 X 1 X F (ω) −(N −3) 8F (ω) −N Uj (0, ω) = 1[ω falsifies Cj ] · 2−(N −3) = ·2 = ·2 , m j=1 m j=1 m m

and P(1, ω) = 0 by definition of Uj . Remark that Q(0, ω) = λ · 2−N =

1 8F (ω) −N · 2−N < ·2 = P(0, ω) 2m m

iff F (ω) ≥ 1;

Q(1, ω) = (1 − λ) · 2−N > 0 = P(1, ω). Therefore dTV (P, Q) =

max {0, P(b, ω) − Q(b, ω)} =

X (b,ω)∈{0,1}N +1

=

X

P(0, ω) − Q(0, ω)

ω∈{0,1}N :F (ω)≥1

X ω∈{0,1}N :F (ω)≥1

8F (ω) −N ·2 − m

X ω∈{0,1}N :F (ω)≥1

1 · 2−N . 2m

(14)

For the first term, we consider counting all pairs (ω, Cj ) such that ω falsifies Cj . On the one hand, for each ω, it falsifies F (ω) clauses; on the other hand, for each Cj , it is falsified by 2N −3 assignments. Thus X X F (ω) = F (ω) = m · 2N −3 . ω∈{0,1}N :F (ω)≥1

ω∈{0,1}N

For the second term, there are S assignments ω satisfying all clauses, that is, F (ω) = 0. So there are 2N − S assignments satisfying F (ω) ≥ 1. Substituting back to (14), we have dTV (P, Q) =

X ω∈{0,1}N :F (ω)≥1

=

2−N +3 m

8F (ω) −N ·2 − m

· (m · 2N −3 ) −

2−N 2m

X ω∈{0,1}N :F (ω)≥1

· (2N − S) = 1 −

1 · 2−N 2m

1 2−N + S. 2m 2m

Thus the number of satisfying assignments is given by #SAT(φ) =

S 2N −r

= 2r (2m · dTV (P, Q) − 2m + 1).

The constructed distributions live on the n = N + 1 dimensional subcube, with k1 + k2 = m + 2 total components. Since N = max{r, m} and r ≤ 3m, we have m ≤ N ≤ 3m, hence k1 + k2 = m + 2 = Θ(N + 1) = Θ(n). Therefore, an algorithm for Problem 1.3 with k1 +k2 = Θ(n) yields an exact algorithm for #3SAT. 23

Acknowledgements Weiming Feng acknowledges the support of ECS grant 27202725 from Hong Kong RGC. We thank Arnab Bhattacharyya and Guy Van den Broeck for bringing the problem to our attention.

References [ABH+ 18] Hassan Ashtiani, Shai Ben-David, Nicholas J. A. Harvey, Christopher Liaw, Abbas Mehrabian, and Yaniv Plan. Nearly tight sample complexity bounds for learning mixtures of Gaussians via sample compression schemes. In NeurIPS, pages 3416–3425, 2018. [BGM+ 23] Arnab Bhattacharyya, Sutanu Gayen, Kuldeep S. Meel, Dimitrios Myrisiotis, A. Pavan, and N. V. Vinodchandran. On approximating total variation distance. In IJCAI, pages 3479–3487. ijcai.org, 2023. [BGM+ 24] Arnab Bhattacharyya, Sutanu Gayen, Kuldeep S. Meel, Dimitrios Myrisiotis, A. Pavan, and N. V. Vinodchandran. Total variation distance meets probabilistic inference. In ICML. OpenReview.net, 2024. [BGM+ 25] Arnab Bhattacharyya, Sutanu Gayen, Kuldeep S. Meel, Dimitrios Myrisiotis, Aduri Pavan, and N. V. Vinodchandran. Computational explorations of total variation distance. In ICLR. OpenReview.net, 2025. [BGMV20] Arnab Bhattacharyya, Sutanu Gayen, Kuldeep S. Meel, and N. V. Vinodchandran. Efficient distance approximation for structured high-dimensional distributions via learning. In NeurIPS, 2020. [BLMT23] Guy Blanc, Jane Lange, Ali Malik, and Li-Yang Tan. Lifting uniform learners via distributional decomposition. In STOC, pages 1755–1767. ACM, 2023. [BLST25] Guy Blanc, Jane Lange, Carmen Strassle, and Li-Yang Tan. A distributional-lifting theorem for PAC learning. In Nika Haghtalab and Ankur Moitra, editors, COLT, Proceedings of Machine Learning Research, pages 375–379. PMLR, 2025. [CK14] Taolue Chen and Stefan Kiefer. On the total variation distance of labelled markov chains. In LICS, pages 33:1–33:10. ACM, 2014. [CM19] Sitan Chen and Ankur Moitra. Beyond the low-degree algorithm: mixtures of subcubes and their applications. In STOC, pages 869–880. ACM, 2019. [CR14] Clément L. Canonne and Ronitt Rubinfeld. Testing probability distributions underlying aggregated data. In ICALP, volume 8572 of Lecture Notes in Computer Science, pages 283–295. Springer, 2014. [FGJW23] Weiming Feng, Heng Guo, Mark Jerrum, and Jiaheng Wang. A simple polynomialtime approximation algorithm for the total variation distance between two product distributions. TheoretiCS, 2, 2023. [FLL24] Weiming Feng, Liqiang Liu, and Tianren Liu. On deterministically approximating total variation distance. In SODA, pages 1766–1791. SIAM, 2024. 24

[FLY25] Weiming Feng, Hongyang Liu, and Minji Yang. Approximating the total variation distance between spin systems. In COLT, volume 291 of Proceedings of Machine Learning Research, pages 1974–2025. PMLR, 2025. [FOS08] Jon Feldman, Ryan O’Donnell, and Rocco A. Servedio. Learning mixtures of product distributions over discrete domains. SIAM, 37(5):1536–1564, 2008. [GJM+ 24] Spencer L. Gordon, Erik Jahn, Bijan Mazaheri, Yuval Rabani, and Leonard J. Schulman. Identification of mixtures of discrete product distributions in near-optimal sample and time complexity. In COLT, Proceedings of Machine Learning Research, pages 2071–2091. PMLR, 2024. [GMRS21] Spencer Gordon, Bijan H. Mazaheri, Yuval Rabani, and Leonard J. Schulman. Source identification for mixtures of product distributions. In COLT, Proceedings of Machine Learning Research, pages 2193–2216. PMLR, 2021. [JO14] Prateek Jain and Sewoong Oh. Learning mixtures of discrete product distributions using spectral decompositions. In COLT, JMLR Workshop and Conference Proceedings, pages 824–856. JMLR.org, 2014. [KECM23] Jeongyeol Kwon, Yonathan Efroni, Constantine Caramanis, and Shie Mannor. RewardMixing MDPs with few latent contexts are learnable. In ICML, volume 202 of Proceedings of Machine Learning Research, pages 18057–18082. PMLR, 2023. [Kie18] Stefan Kiefer. On computing the total variation distance of hidden Markov models. In ICALP, volume 107 of LIPIcs, pages 130:1–130:13. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2018. [KTT23] Chinmaya Kausik, Kevin Tan, and Ambuj Tewari. Learning mixtures of Markov Chains and MDPs. In ICML, Proceedings of Machine Learning Research, pages 15970–16017. PMLR, 2023. [LS17] Jerry Li and Ludwig Schmidt. Robust and proper learning for mixtures of Gaussians via systems of polynomial inequalities. In COLT, Proceedings of Machine Learning Research, pages 1302–1382. PMLR, 2017. [SV03] Amit Sahai and Salil Vadhan. A complete problem for statistical zero knowledge. J. ACM, 50(2):196–249, 2003.

A

Validity of Sampling and Coupling Processes

A.1

Validity of the sampling process

Proof of Lemma 4.3. We prove by induction on j that for any j ∈ {1, 2, . . . , n + 1} and any mixing 1 weights ᾱ = (ᾱ(s) )ks=1 , RecursiveSample(j, ᾱ) outputs a random sample Xj:n ∼ P̄ ≜

k1 X s=1

25

ᾱ(s)

n O (s)

Pi .

i=j

The base case j = n + 1 is immediate since the algorithm returns ∅. For the induction step, fix any ωj:n ∈ [q]n+1−j . Let c = ωj . In Line 4, the algorithm chooses (c, 0) with probability ℓP (j, ᾱ, c) and chooses (c, 1) with probability P̄j (c) − ℓP (j, ᾱ, c). Conditioned on e = 0, by the induction P 1 Q (s) hypothesis the recursive call returns ωj+1:n with probability ks=1 ᾱ(s) ni=j+1 Pi (ωi ). Conditioned (s)

1 on e = 1, the algorithm updates the mixing weights to ᾱj,c = (ᾱj,c )ks=1 and again by the induction

hypothesis the recursive call returns ωj+1:n with probability Pr [Xj:n = ωj:n ] = ℓP (j, ᾱ, c)

k1 X

ᾱ(s)

s=1

n Y

(s) Qn (s) i=j+1 Pi (ωi ). Therefore, s=1 ᾱj,c

Pk1

k1 n
X (s) Y

(s)

Pi (ωi ) + P̄j (c) − ℓP (j, ᾱ, c)

ᾱj,c

s=1

i=j+1

(s)

Pi (ωi )

i=j+1

= P̄(ωj:n ), where the second equality is exactly the decomposition in (7). Note that (5) guarantees P̄j (c) − ℓP (j, ᾱ, c) =





(s)

ᾱ(s) Pj (c) − ℓP (j, ᾱ, c) ≥ 0,

X s∈[k1 ]:ᾱ(s) >0

which completes the induction.

A.2

Validity of the coupling

Proof of Lemma 4.5. We prove a slightly more general statement: for any j ∈ {1, 2, . . . , n + 1} and 1 2 any mixing-weight vectors ᾱ = (ᾱ(s) )ks=1 and β̄ = (β̄ (t) )kt=1 , RecursiveCouple(j, ᾱ, β̄) outputs joint random variables (Xj:n , Yj:n ) such that Xj:n ∼ P̄ ≜

k1 X s=1

ᾱ(s)

n O (s)

Yj:n ∼ Q̄ ≜

Pi ,

k2 X t=1

i=j

β̄ (t)

n O

(t)

Qi .

i=j

1 2 Taking j = 1 with the initial weights α = (α(s) )ks=1 and β = (β (t) )kt=1 gives the lemma. We now prove Xj:n ∼ P̄ by induction on j. The proof for Yj:n ∼ Q̄ is symmetric. We proceed by induction on j. The base case is j = n + 1. The algorithm returns (∅, ∅), which matches the (trivial) distributions on the empty product space. Fix any j ≤ n. In Line 3 of Algorithm 2, the algorithm samples (cP̄ , eP̄ ) and (cQ̄ , eQ̄ ) using the ¯ By construction of C, ¯ the marginal distribution of (c , e ) is exactly as follows: for coupling C. P̄ P̄ each c ∈ [q], Pr [(cP̄ , eP̄ ) = (c, 0)] = ℓ(j, ᾱ, β̄, c) and Pr [(cP̄ , eP̄ ) = (c, 1)] = P̄j (c) − ℓ(j, ᾱ, β̄, c). All the above probabilities are non-negative due to Definition 4.4. Hence, ignoring the other side, the pair-sampling rule on the P̄ side matches exactly the rule in RecursiveSample (Line 4 of Algorithm 1). Now fix any ωj:n ∈ [q]n+1−j and write c = ωj . We compute Pr [Xj:n = ωj:n ] by conditioning on eP̄ ∈ {0, 1}:

Pr [Xj:n = ωj:n ] = Pr [(cP̄ , eP̄ ) = (c, 0)] · Pr [Xj+1:n = ωj+1:n | (cP̄ , eP̄ ) = (c, 0)] + Pr [(cP̄ , eP̄ ) = (c, 1)] · Pr [Xj+1:n = ωj+1:n | (cP̄ , eP̄ ) = (c, 1)] . If eP̄ = 0, the algorithm recurses with RecursiveCouple(j + 1, ᾱ, β̄); by the induction hypothesis, P 1 Q (s) the marginal distribution of Xj+1:n is given by ks=1 ᾱ(s) ni=j+1 Pi (ωi ). If eP̄ = 1, the algorithm (s)

1 updates the weights to ᾱj,c = (ᾱj,c )ks=1 . If cP̄ = cQ̄ , then the algorithm recurses via RecursiveCouple;

26

(s) Q

(s)

1 by the induction hypothesis, the marginal law of Xj+1:n = ωj+1:n is ks=1 ᾱj,c ni=j+1 Pi (ωi ). Otherwise, the algorithm samples Xj+1:n via RecursiveSample; by Lemma 4.3, its law is again Pk1 (s) Qn (s) i=j+1 Pi (ωi ). Plugging these together with the above marginal distribution of (cP̄ , eP̄ ) s=1 ᾱj,c into the decomposition, we obtain Pr [Xj:n = ωj:n ] = P̄(ωj:n ) by the same argument as in the proof of Lemma 4.3. Therefore, Xj:n ∼ P̄.

P

B

Discrepancy of the Recursive Coupling

Proof of Lemma 4.6. Consider a step (j, ᾱ, β̄) in the RecursiveCouple(j, ᾱ, β̄) process with P̄ = Pk1 (s) Nn P(s) and Q̄ = Pk2 β̄ (t) Nn Q(t) . Let k(ᾱ) and k(β̄) denote the number of nonzero i=j i i=j i s=1 ᾱ t=1 1 2 weights in (ᾱ(s) )ks=1 and (β̄ (t) )kt=1 , respectively. Formally, k(ᾱ) = |{ᾱ(s) > 0 | 1 ≤ s ≤ k1 }| and k(β̄) is defined similarly. For j, ᾱ, β̄, we define a parameter δ(j, ᾱ, β̄) as follows: 



δ j, ᾱ, β̄ ≜

max

Λ∈S,ω∈[q]Λ

P̄Λ (ω) − Q̄Λ (ω) =

k1 X

max

Λ∈S,ω∈[q]Λ s=1

ᾱ(s)

Y (s)

k2 X

i∈Λ

t=1

Pi (ωi ) −

β̄ (t)

Y

(t)

Qi (ωi )

i∈Λ

where the maximum is taken over all subsets Λ ∈ S such that n

o

S ≜ Λ ⊆ {j, j + 1, . . . , n} | |Λ| ≤ k(ᾱ) + k(β̄) − 1 .

(15) 

Note that both P̄ and Q̄ are defined over the indices {j, j + 1, . . . , n}. The definition of δ j, ᾱ, β̄



enumerates all possible subsets Λ of size at most k(ᾱ) + k(β̄) − 1 and finds the configuration ω ∈ [q]Λ that the  maximizes   absolute  difference between P̄Λ (ω) and Q̄Λ (ω). Hence, for every such Λ, we have δ j, ᾱ, β̄ ≤ dTV P̄Λ , Q̄Λ , and by the data processing inequality for total variation distance, we have 







(16)

δ j, ᾱ, β̄ ≤ dTV P̄, Q̄ .

The algorithm RecursiveCouple defines a coupling C¯RE between P̄ and Q̄. Let Xj:n ∼ P̄ and Yj:n ∼ Q̄ denote the output of the coupling algorithm. Let 



F j, ᾱ, β̄ ≜ PrC¯RE [Xj:n ̸= Yj:n ] ,

where Xj:n ∼ P̄ and Yj:n ∼ Q̄.

1 2 We next prove that for the recursive coupling step RecursiveCouple(j, (ᾱ(s) )ks=1 , (β̄ (t) )kt=1 ),









F j, ᾱ, β̄ ≤ (4(n − j + 1)q)k(ᾱ)+k(β̄)−1 · δ j, ᾱ, β̄ .

(17)

Note that (17) implies the lemma by plugging in j = 1, α = ᾱ, and β = β̄, and using the fact that k(α) + k(β) ≤ k1 + k2 and δ(1, α, β) ≤ dTV (P, Q) by (16). We apply mathematical induction on j from n + j = n + 1, RecursiveCouple  1 to 1. When  always returns (Xj:n , Yj:n ) = (∅, ∅), which implies F n + 1, ᾱ, β̄ = 0, and (17) holds trivially. For the inductive step, fix 1 ≤ j ≤ n and we assume that (17) holds for all j ′ with j + 1 ≤ j ′ ≤ n + 1. We need to show that (17) holds for j. We can decompose the probability as follows: 



F j, ᾱ, β̄ = PrC¯RE [Xj:n ̸= Yj:n ] =

X

Pr [(cP̄ , eP̄ ) = (c, 0)] PrC¯RE [Xj:n ̸= Yj:n | (cP̄ , eP̄ ) = (c, 0)]

c∈[q]

+

X

Pr [(cP̄ , eP̄ ) = (c, 1)] PrC¯RE [Xj:n ̸= Yj:n | (cP̄ , eP̄ ) = (c, 1)] .

c∈[q]

27

In the above equation, we partition the probability by enumerating all possible values of cP̄ ∈ [q] and eP̄ ∈ {0, 1}. The equation follows from the law of total probability. The event (cP̄ , eP̄ ) = (c, 0) happens with probability ℓ(j, ᾱ, β̄, c). By the definition of the coupling, cQ̄ = c and eQ̄ = 0 in this case, and the algorithm proceeds recursively to Line 5. For the case (cP̄ , eP̄ ) = (c, 1), we can further divide it into two subcases: (1) cP̄ = cQ̄ = c, which happens with probability min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c), and the algorithm proceeds recursively to Line 9; (2) cP̄ ̸= cQ̄ , which happens with probability max{P̄j (c) − Q̄j (c), 0}, and the coupling fails (Xj ̸= Yj implies X ̸= Y ). Overall, we have the following recursion: 



F j, ᾱ, β̄ ≤

X



ℓ(j, ᾱ, β̄, c) · F j + 1, ᾱ, β̄



c∈[q]

X

+





min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) · F j + 1, ᾱj,c , β̄j,c



c∈[q]

+

max{P̄j (c) − Q̄j (c), 0}.

X c∈[q]

We analyze the three terms on the right-hand side one byone. By the induction hypothesis, the first term is at most (4(n − j)q)k(ᾱ)+k(β̄)−1 · δ j + 1, ᾱ, β̄ . For the second term, the induction 







hypothesis gives F j + 1, ᾱj,c , β̄j,c ≤ (4(n − j)q)k(ᾱj,c )+k(β̄j,c )−1 · δ j + 1, ᾱj,c , β̄j,c . Moreover, to upper bound F (j, ᾱ, β̄), we only consider c ∈ [q] with min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) ̸= 0. The  third term is the total variation distance between P̄j and Q̄j . By the definition of δ j, ᾱ, β̄ , 



taking Λ = {j} and ωj = c implies |P̄j (ω) − Q̄j (ω)| ≤ δ j, ᾱ, β̄ . Hence, the third term is at most 



q · δ j, ᾱ, β̄ . Overall, we have the following recursion: 



F j, ᾱ, β̄ ≤

X



ℓ(j, ᾱ, β̄, c) · (4(n − j)q)k(ᾱ)+k(β̄)−1 · δ j + 1, ᾱ, β̄



c∈[q]

X





min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) · (4(n − j)q)k(ᾱj,c )+k(β̄j,c )−1 δ j + 1, ᾱj,c , β̄j,c

+



c∈[q]





+ q · δ j, ᾱ, β̄ . Recall that the updated parameters when min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) > 0: (s)

(s)

ᾱj,c ≜ ᾱ(s)







(t)  and β̄j,c ≜ β̄ (t)

Pj (c) − ℓ j, ᾱ, β̄, c P̄j (c) − ℓ j, ᾱ, β̄, c

(t)







 .

Qj (c) − ℓ j, ᾱ, β̄, c Q̄j (c) − ℓ j, ᾱ, β̄, c (s)

(t)

Remark that if ᾱ(s) = 0 or β̄ (t) = 0, then the updated mixing weight ᾱj,c = 0 or β̄j,c = 0, that is k(ᾱj,c ) ≤ k(ᾱ) and k(β̄j,c ) ≤ k(β̄). By Definition 4.4, there must exist some s ∈ [k1 ] with ᾱ(s) > 0 or (s) (t) t ∈ [k2 ] with β̄ (t) > 0, such that Pj (c) = ℓ(j, ᾱ, β̄, c) or Qj (c) = ℓ(j, ᾱ, β̄, c), and thus there must (s)

(t)

exist ᾱ(s) > ᾱj,c = 0 or β̄ (t) > β̄j,c = 0. Hence, if the recursion goes to Line 9 with cP̄ = cQ̄ = c, then it must hold that min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) > 0 and we have k(ᾱj,c ) + k(β̄j,c ) ≤ k(ᾱ) + k(β̄) − 1. Next, we claim the following inequality: 

 







min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) δ j + 1, ᾱj,c , β̄j,c ≤ 3δ j, ᾱ, β̄ . 28

(18)

Let us first assume that (18) holds. Then we have 













F j, ᾱ, β̄ ≤ (4(n − j)q)k(ᾱ)+k(β̄)−1 · δ j + 1, ᾱ, β̄ + 3q (4(n − j)q)k(ᾱ)+k(β̄)−2 δ j, ᾱ, β̄ + qδ j, ᾱ, β̄ 



≤ (4(n − j + 1)q)k(ᾱ)+k(β̄)−1 δ j, ᾱ, β̄ , where the last inequality holds because δ(j + 1, ᾱ, β̄) ≤ δ(j, ᾱ, β̄). This finishes the induction step. Finally, we prove (18). Note that (18) holds trivially when min{P̄j (c), Q̄j (c)} − ℓ(j, ᾱ, β̄, c) = 0. We assume the case ℓ(j, ᾱ, β̄, c) < P̄j (c) ≤ Q̄j (c), and the other case, when ℓ(j, ᾱ, β̄, c) < Q̄j (c) ≤ P̄j (c), can be proved similarly. By the definition of δ(j + 1, ᾱj,c , β̄j,c ) and (15), one needs to enumerate all possible subsets Λ′ ∈ {j + 1, . . . , n} with size at most k(ᾱj,c ) + k(β̄j,c ) − 1 and all ′ configurations ω ′ ∈ [q]Λ . Note that k(ᾱj,c ) + k(β̄j,c ) ≤ k(ᾱ) + k(β̄) − 1. We consider all Λ′ with size at most k(ᾱ) + k(β̄) − 2, which covers all possible subsets Λ′ in the definition of δ(j + 1, ᾱj,c , β̄j,c ). Hence, we can get an upper bound on δ(j + 1, ᾱj,c , β̄j,c ). Formally, for any Λ′ ⊆ {j + 1, . . . , n}, ′ |Λ′ | ≤ k(ᾱ) + k(β̄) − 2, and any ω ′ ∈ [q]Λ , we have k1 X (s) Y

ᾱj,c

i∈Λ′

s=1

=

k1 X

(s)

Pi (ωi′ ) −

k2 X (t) Y

β̄j,c

t=1

(t)

Qi (ωi′ )

i∈Λ′

(s)

ᾱ(s)

s=1

Pj (c) − ℓ(j, ᾱ, β̄, c) Y P̄j (c) − ℓ(j, ᾱ, β̄, c) i∈Λ′

(s)

Pi (ωi′ ) −

k2 X

(t)

β̄ (t)

t=1

k

(∗)

≤

Qj (c) − ℓ(j, ᾱ, β̄, c) Y Q̄j (c) − ℓ(j, ᾱ, β̄, c) i∈Λ′

(t)

Qi (ωi′ )

k

1 2 Y (s) X X Y (t) 1 (s) (t) ᾱ(s) Pj (c) Pi (ωi′ ) − β̄ (t) Qj (c) Qi (ωi′ ) P̄j (c) − ℓ(j, ᾱ, β̄, c) s=1 t=1 i∈Λ′ i∈Λ′

k

+

(19)

k

1 2 Y (s) X X Y (t) ℓ(j, ᾱ, β̄, c) ᾱ(s) Pi (ωi′ ) − β̄ (t) Qi (ωi′ ) P̄j (c) − ℓ(j, ᾱ, β̄, c) s=1 t=1 i∈Λ′ i∈Λ′

(20)

k

+

2   Y X 1 1 (t) (t) − β̄ (t) Qj (c) − ℓ(j, ᾱ, β̄, c) Qi (ωi′ ) P̄j (c) − ℓ(j, ᾱ, β̄, c) Q̄j (c) − ℓ(j, ᾱ, β̄, c) t=1 i∈Λ′ (21)

The inequality (∗) is obtained by first inserting two terms, −

Pk2

(t)

Q

(c)−ℓ(j,ᾱ,β̄,c) Q

(t) j t=1 β̄ P̄ (c)−ℓ(j,ᾱ,β̄,c) j

Pk2

(t)

Q

(c)−ℓ(j,ᾱ,β̄,c) Q (t) ′ i∈Λ′ Qi (ωi ) and j

(t) j t=1 β̄ P̄ (c)−ℓ(j,ᾱ,β̄,c)

(t) ′ i∈Λ′ Qi (ωi ), whose sum is 0, and then applying the triangle inequality.

Our task is now reduced to bounding the three terms in (19), (20), and (21). For (19), consider the definition of δ(j, ᾱ, β̄). We can take Λ = {j} ∪ Λ′ and ω = (ωj ← c) ∪ ω ′ in the definition. Note that Λ is allowed in the definition of δ(j, ᾱ, β̄) because |Λ′ | ≤ k(ᾱ) + k(β̄) − 2. Then 



k1 k2 δ j, ᾱ, β̄ X Y (t) Y (s) X 1 (s) (t) ᾱ(s) Pj (c) Pi (ωi′ ) − β̄ (t) Qj (c) Qi (ωi′ ) ≤ . P̄j (c) − ℓ(j, ᾱ, β̄, c) s=1 P̄j (c) − ℓ(j, ᾱ, β̄, c) t=1 i∈Λ′ i∈Λ′

29



For (20), as ℓ(j, ᾱ, β̄, c) ≤ 1, consider Λ = Λ′ and ω = ω ′ in the definition of δ(j, ᾱ, β̄); we have k1 X



k2 X

ℓ(j, ᾱ, β̄, c) (t) (s) Qi (ωi′ ) ≤ Pi (ωi′ ) − β̄ (t) ᾱ(s) P̄j (c) − ℓ(j, ᾱ, β̄, c) s=1 t=1 i∈Λ′ i∈Λ′ Y

Y

ℓ(j, ᾱ, β̄, c) · δ j, ᾱ, β̄ P̄j (c) − ℓ(j, ᾱ, β̄, c) 

δ j, ᾱ, β̄ ≤ (t)





P̄j (c) − ℓ(j, ᾱ, β̄, c)

.

(t)

For (21), as Qi (ωi′ ) ≤ 1 and Qj (c) ≥ ℓ(j, ᾱ, β̄, c) for t ∈ [k2 ] with β̄ (t) > 0, k2 X

β̄

(t)

t=1



(t) Qj (c) − ℓ(j, ᾱ, β̄, c)

k2   X (t) (t) ′ Qi (ωi ) ≤ β̄ (t) Qj (c) − ℓ(j, ᾱ, β̄, c) = Q̄j (c) − ℓ(j, ᾱ, β̄, c), t=1 i∈Λ′

 Y

thus k

2   Y X 1 1 (t) (t) − β̄ (t) Qj (c) − ℓ(j, ᾱ, β̄, c) Qi (ωi′ ) P̄j (c) − ℓ(j, ᾱ, β̄, c) Q̄j (c) − ℓ(j, ᾱ, β̄, c) t=1 i∈Λ′

≤

  1 1 Q̄j (c) − ℓ(j, ᾱ, β̄, c) − P̄j (c) − ℓ(j, ᾱ, β̄, c) Q̄j (c) − ℓ(j, ᾱ, β̄, c) 



δ j, ᾱ, β̄ P̄j (c) − Q̄j (c) = ≤ . P̄j (c) − ℓ(j, ᾱ, β̄, c) P̄j (c) − ℓ(j, ᾱ, β̄, c) Then (18) follows by combining the three inequalities above.

30