SPADE: Sequential-clustering Particle Annihilation via Discrepancy Estimation

For an empirical signed measure $\mu = \frac{1}{N} \left(\sum_{i=1}^P \delta_{x_i} - \sum_{i=1}^M \delta_{y_i}\right)$, particle annihilation (PA) removes $N_A$ particles from both $\{x_i\}_{i=1}^P$ and $\{y_i\}_{i=1}^M$ simultaneously, yielding another empirical signed measure $\nu$ such that $\int f d \nu$ approximates to $\int f d \mu$ within an acceptable accuracy for suitable test functions $f$. Such annihilation of particles carrying opposite importance weights has been extensively utilized for alleviating the numerical sign problem in particle simulations. In this paper, we propose an algorithm for PA in high-dimensional Euclidean space based on hybrid of clustering and matching, dubbed the Sequential-clustering Particle Annihilation via Discrepancy Estimation (SPADE). It consists of two steps: Adaptive clustering of particles via controlling their number-theoretic discrepancies, and independent random matching among positive and negative particles in each cluster. Both deterministic error bounds by the Koksma-Hlawka inequality and non-asymptotic random error bounds by concentration inequalities are proved to be affected by two factors. One factor measures the irregularity of point distributions and reflects their discrete nature. The other relies on the variation of test function and is influenced by the continuity. Only the latter implicitly depends on dimensionality $d$, implying that SPADE can be immune to the curse of dimensionality for a wide class of test functions. Numerical experiments up to $d=1080$ validate our theoretical discoveries.


Introduction
Given two sequences X = (x 1 , . . . , x P ) ⊂ Q P and Y = (y 1 , . . . , y M ) ⊂ Q M in a finite rectangular domain Q ⊂ R d , X and Y are sets of positive and negative particles, carrying opposite weights 1 and −1, respectively. An empirical signed measure µ is given by where N > 0 is a prescribed normalizing constant and δ x is Dirac measure concentrated at x. Typical choices of N include N = N tot with N tot = P + M the total particle number or N = P −M for P > M in order to make µ(Q) = 1. Removing N A positive particles from X and N A negative ones from Y simultaneously yields two subsets X A = (x 1 , . . . ,x P −N A ) ⊂ X and Y A = (ỹ 1 , . . . ,ỹ M −N A ) ⊂ Y and another empirical signed measure Of direct interest in various applications is seeking a strategy to diminish the error function for a suitable class of test functions f . The aforementioned problem is called particle annihilation (PA) since positive and negative particles are paired and annihilated, so that the opposite importance weights are cancelled out. It is also known as particle cancelation or particle resampling on some occasions. PA is originated from a vast number of Monte Carlo applications ranging from classical particle transport to quantum Monte Carlo. In classical regime, the signed measure describes fluctuations over the equilibrium. 1,33 While in quantum mechanics, negative weights emerge in Monte Carlo integration of oscillatory and determinantal functions, including the path-integral representation of quantum observables, 4 the pseudo-differential operator in the Wigner dynamics 13,25 and the Rayleigh quotients with Slater determinant ansatz of fermionic wave-functions. 22 Unfortunately, stochastic estimators of the form in Eq. (1.1) usually suffer from the notorious numerical sign problem due to the near-cancelation of contributions from positive and negative weights, which stem from oscillating functions or even and odd permutations in determinants. 9 As a consequence, it leads to an exponential increase of particle number and stochastic variance, as well as the computational complexity, especially when simulating long-time particle dynamics. 4,22,25,33 The sign problem is believed to be NP-hard since there might not exist an algorithm to achieve a relative statistical error scaling polynomially with the dimensionality d and the system size in general. 12, 30 A well-known example is given in 30 by mapping a 3-D Ising spin glass model, which belongs to the complexity class NP, to a quantum system with the sign problem.
Many efforts have been made to potentially overcome the numerical sign problem in certain situations, albeit usually in the settings that rely heavily on concrete problems. Several related approaches include the particle resampling in kinetic theory by filtering out high-frequency components, 33 the grid-based annihilation in quantum Boltzmann simulations, 13, 25 the fixed-node approximation in diffusion Monte Carlo, 22 the annihilation of determinants in Fermion Monte Carlo, 4 the resummation in an inchworm algorithm 7 and the stationary phase approximation in the Wigner branching random walk. 26 In spite of their distinct appearances, the core of all the above approaches is essentially the same, that is, to fully utilize the cancelation of positive and negative weights, thereby reducing the growth of particle number as well as stochastic variance. This creed will be faithfully inherited in subsequent trials for a new approach.
Our approach, termed the Sequential-clustering Particle Annihilation via Discrepancy Estimation (SPADE), hybridizes sequential clustering of particles and independent random matching in each cluster, as depicted by the diagram: The SPADE algorithm utilizes both number-theoretic and statistical properties of point distributions. Detailedly speaking, SPADE consists of two steps. The first step is to perform adaptive clustering of positive and negative particles through a sequential binary partition of the domain Q. Each time we pick up a cluster of particles and make a further split until their star discrepancies, say, the number-theoretic irregularities, reach a prescribed threshold, partially borrowing the pioneering idea from the non-parametric density estimation method based on discrepancy sequential partitioning. 16 Then it succeeds at the second step by seeking independent random matchings in each cluster, which is essentially equivalent to sampling from a finite population without replacement, and removing paired particles. In contrast to existing approaches for PA, such as 'fixed-node approximation', 22 SPADE requires no a priori information of the exact nodal hypersurface of underlying integrands. Instead, it will 'learn' the nodes based on the point distribution in an adaptive manner. In a sense, SPADE can be intuitively treated as a 'flexible-node approximation'.
From the perspective of numerical analysis, it requires to bound the error function E(f ) and to ensure its convergence as N → ∞. The interplay between the continuous integrals and discrete (combinatorial) nature of point distributions will be extensively utilized in analysis of SPADE, and the discrepancy theory directly links both sides. First, the deterministic error bounds under arbitrary matching are given with the help of the famous Koksma-Hlawka inequality 8,15 and a combinatorial property of the star discrepancy. The obtained bound can be separated into two parts: one comes from irregularities of point distributions (measured by the star discrepancy) and the other from the variation of test function f in the sense of Hardy and Krause, without any a priori knowledge of continuity. Since only the latter implicitly depends on the dimensionality d, the error bounds may be immune to the curse of dimensionality for some good test functions. This explains why SPADE is able to achieve a satisfactory numerical accuracy for certain problems up to d = 1080 and work well for slow-varying test functions without continuous derivates. Second, non-asymptotic stochastic error bounds under random matching further sharpens the deterministic ones with the help of concentration inequalities for the sum of independent random variables 3, 5, 11 and sampling from a finite population without replacement. 2,24,29 As a comparison, the optimal assignment 14 is also adopted to pair the positive and negative particles in each cluster. Numerical results show that the random matching usually outperforms the optimal assignment.
The rest is organized as follows. Notations and definitions are given in Section 2, proceeding with a statement of main results in Section 3. The sequential clustering is issued in Section 4 with a proof of deterministic bounds. The random matching is discussed in Section 5 and the bounds of error function will be sharpened thanks to concentration of random variables. It follows by numerical experiments in Section 6. The paper is concluded in Section 7 with a few remarks.

Notations and definitions
We use the short hand notation Q = [a, b] = d j=1 [a j , b j ] to denote a hyperrectangle in R d , where a = (a 1 , . . . , a d ) and b = (b 1 , . . . , b d ) with a j and b j specifying the lower and upper bound of Q along dimension j. A partition of Q at level K is denoted by by sorting the i-th coordinates of all the lower and upper bounds of Q k in an increasing order. The collection of all partitions is denoted by P.
The subsequences of X and Y in Q k are denoted by M k ), respectively. P k = |X k | and M k = |Y k | are the cardinal numbers, and we let P k ∧ M k = min{P k , M k }.
which is equivalent to sampling M times from a finite population {1, . . . , P } without replacement. By contrast, an optimal assignment is a matching σ to attain min σ M i=1 x σ(i) − y i for a given norm. To quantitively measure the difference between the continuous integral and the arithmetic mean of the sequences {f (x 1 ), . . . , f (x P )}, i.e., we need both concepts of the variation of the function f and the star discrepancy of sequence (x 1 , . . . , x P ). Here λ(Q) gives the Lesbesgue measure of Q. Given a partition of Q, we define an operator (2. 3) The

4)
where the supremum is extended over all possible partitions P. Then the variation of f on [a, b] in the sense of Hardy and Krause is . . , i l }, and the second sum is extended over all l-dimensional faces F l . In particular, for continuous functions f with continuous mixed partial derivates, we have For functions with vanishing mixed derivates, e.g., b] (f ) may depend on d linearly (see Eqs. (6.7) and (6.8)). By contrast, for f ( 20 Another characterization of local variation of f on [a, b] is the oscillator norm: and the local L 2 -norm of f on Q is defined by The star discrepancy, 8,15 stemming from number theory, is adopted to measure the irregularity of a sequence (x 1 , . . . , Definition 3. Let X = (x 1 , . . . , x P ) be a finite sequence of points in the d-dimensional space [0, 1] d . Then the number is called the star discrepancy of sequence X . Here A(I, P, X ) = P i=1 χ I (x i ) gives the number of the sequence X in I.
For a general sequence (x 1 , . . . , x P ) ⊂ [a k , b k ], its irregularity can be measured by the star discrepancy after mapping it from [a k , b k ] to [0, 1] d via a linear scaling: In this work, we always deal with test functions in the class with bounded variations in the sense of Hardy and Krause, say, V HK [a,b] (f ) < ∞. Under the partition P K , their oscillator norms osc [a k ,b k ] (f ) are also bounded for k = 1, . . . , K.

Statement of main theorems
The SPADE algorithm consists of two steps: Clustering and matching. Clustering: Establish a partition P K of Q of level K associated with clustering of particles X and Y such that the star discrepancies of both clustered subsequences X k and Y k in Q k are all bounded for k = 1, . . . , K: where ϑ ∈ (0, 1) is a prescribed parameter to control the error function E(f ).
Matching: Within each cluster, seek a matching σ from {1, . . . , M k } to {1, . . . , P k } independently when M k ≤ P k , or from {1, . . . , P k } to {1, . . . , M k } when M k > P k , and then remove the matched pairs from X k and Y k .
Under the partition P K , the error function between two empirical signed measures µ (before annihilation) and ν (after annihilation) is The condition (3.1) is used to bound the second and third terms in Eq. (3.3) owing to the Koksma-Hlawka inequality. As for the first term, the bound can be obtained by either the combinatorial property of the star discrepancy or the bounded difference condition osc [a k ,b k ] (f ) < ∞. It deserves to mention that K, the partition level, depends on d implicitly as the star discrepancy does. In general, K has to be smaller than total particle number N tot to avoid the 'overfitting' phenomenon.
The first bound (3.5), obtained using the Koksma-Hlawka inequality and the combinatorial property of star discrepancies, characterizes the basic convergence of SPADE. A key observation follows that the sole factor that may implicitly depend on dimensionality d is the variation function V HK [a,b] (f ). In other words, the bound (3.5) might not be bothered by the curse of dimensionality for slow-varying test functions.
The worse-case bound of K k=1 ξ k appears in the second bound (3.6), which turns out to be a gross over-estimate in general. A useful improvement is to adopt independent random matchings in each cluster, so that ξ 1 , . . . , ξ K become independent random variables with Eξ k = 0 and the sum K k=1 ξ k will be concentrated near the mean value zero. Meanwhile, the second and third terms in Eq. (3.3) remain non-stochastic and are bounded by the Koksma-Hlawka inequality. In consequence, stochastic upper bounds can be obtained owing to the concentration bounds of Hoeffding's type and Bernstein's type for sampling with or without replacement. 2,5 The stochastic bounds can dramatically sharpen the deterministic worse-case bound with a high probability, albeit not tight. Furthermore, estimation of the stochastic variance Var(E(f )) = E(E(f ) − EE(f )) 2 is also of the great importance for the random matching. That is equivalent to estimate Theorem 2 (Stochastic bounds under random matching). Under the assumptions in Theorem 1 and independent random matching between X k and Y k in each cluster, for any function f with bounded variation V HK [a,b] (f ) < ∞ in the sense of Hardy and Krause, it holds that

8)
and Var where H 1 (N, γ) = γ 4 + 1 2N 1/2 + 1 4γN . Furthermore, for any ε ∈ (0, 1), with probability higher than 1 − ε, it has that Eqs. (3.8) and (3.9) characterize the bounds for the bias and variance of random matching, respectively. In the sense of expectation, the bias of PA is totally determined by the irregularity of point distributions as Eξ k vanishes under sampling without replacement. It notes that the constant factor is tighter compared to Eq. (3.5) as γ/4 + ϑ 2 /4γ ≥ ϑ/2. The variance of random matching is of the order N −3/2 and the constant factor relies on the oscillator norm of f . The bound (3.10) may dramatically improve (3.6). For the first term, the constant factor is sharpened as Actually, it gains an extra rate of N −1/4 . For the second term, the constant factor is diminished provided The third term keeps unchanged due to the bias of matching. Therefore, with a high probability, the random matching is advantageous due to the concentration property near the mean value K k=1 Eξ k = 0. We shall provide an alternative estimation of Var(E(f )) by taking the irregularity of point distribution into account under an additional assumption that f 2 has a bounded variation. The discrepancy theory will be again used for bounding the local variance of random matching in each cluster, which in turn give an estimate of the second moment K k=1 Eξ 2 k .

Clustering and discrepancy theory
As the first and crucial step of the SPADE algorithm, the clustering of particles is to divide the PA problem into several local parts. The prototype is the girdbased annihilation, 1,13,25 in which the domain is partitioned into bins of equal size, so that the positive and negative particles located in the same bin are annihilated. Apparently, the advantage of such method is its simplicity and the numerical errors are controlled by the bin size. 31 It is also closely related to the histogram statistics in non-parametric density estimation since the annihilation can be realized based on the difference of two probability densities postulated from positive and negative particles, respectively.
However, the setting of uniform grid is no longer applicable for higher dimensional PA problem arising in particle simulations of many-body problems. The reasons are as twofold. First, the number of bins grows exponentially in d, known as the curse of dimensionality. Second, the required sample size for given accuracy also grows with the increasing number of bins. 31 A potential way to resolve these challenges is to replace the uniform grid mesh by an adaptive one driven by the distribution of particle populations, borrowing the basic idea from high-dimensional statistical learning theory. 16,17 However, the negative weights emerged in the PA problem make it essentially different from the density estimation.
In practice, an adaptive partition is established through the sequential binary splitting, which in turn clusters positive and negative particles in a hierarchical manner. Two key ingredients must be specified for the adaptive partition. One is the criterion for stopping the binary splitting. SPADE tries to utilize the numbertheoretical properties (discrepancy) of two kinds of particles to make the decision like in discrepancy sequential partitioning, 16 instead of relying on any a priori probability distribution as adopted in the Bayesian interference. 17 The other is the choice of nodes to be split, which is relevant in monitoring the nodal hyper-surfaces of underlying integrand function. SPADE chooses a node to maximize a gap function between the distributions of positive and negative particles.

Setting of sequential clustering
An adaptive partition is built up through the binary splitting. The binary partition can be defined in the following recursively way. It starts with P 0 = Q. Suppose at level K we have P K = {Q 1 , Q 2 , . . . , Q K } with Q k mutually disjoint and Q = ∪ K k=1 Q k , then the binary partition at level K + 1 is obtained by choosing one sub-rectangle Q k and dividing it into two parts along one of its coordinate, parallel to one of dimensions.
For example, for Q k = [a k , b k ], we pick up j-th dimension and choose a node We can establish the partition P K+1 at level K + 1 by

2)
Algorithm 1 Clustering via adaptive partitioning Choose a node c k,j to attain the maximum of Eq. (4.4) Divide the sub-region Q k into Q (1) k and Q (2) k as given in Eq. (4.1) Divide the pointsets: Although there are various choices of c k,j , a rule of thumb for PA is to dig out the possible nodal points in the hyper-surfaces of the underlying integrand function. To this end, we can adopt the following strategy in practice. We can pick up the j-th dimension and m equidistant points in the interval [a k,j , b k,j ]: then select the node c k,j from (m − 1) × d choices (typically, m = 2, 4, 8) to attain the maximal value of the following difference gap: which measures the gap between the distributions of positive and negative particles after splitting. Heuristically, a large difference gap means positive and negative particles tend to concentrate in different bins. That is, concentrations of positive or negative particles should be observed in opposite sides of the nodal hyper-surfaces.

Deterministic bounds under arbitrary matching
The clustering allows us to analyze the error function E(f ) by bounding the local pieces E k (f ), which depend on both variations of test function f and the bounds for star discrepancies of X k and Y k due to the celebrated Koksma-Hlawka inequality.
In fact, the discrepancy theory plays a crucial role in high-dimensional numerical analysis and the most outstanding application is the low-discrepancy sequence in numerical integration, known as the quasi-Monte Carlo method. 6,19 Owing to the deterministic sequences with the star discrepancy of the order (log n) d /n, where n is the required number of samples, the quasi-Monte Carlo achieves a theoretical convergence order of n −1 , but may still lose its effectiveness in higher dimension due to the factor (log n) d . 6,28 Nonetheless, in certain applications it has been reported to obtain a satisfactory accuracy in evaluating the integrations up to d = 360. 21 Here we would like to focus on the potential application of discrepancy theory in clustering, without any specified construction of low-discrepancy sequences. It deserves to mention that the bounds of star discrepancies given in Eq. (3.1) are still of order n −1 with n the count of sequences, but their numerators are free from dimensionality d. As a consequence, in the deterministic bounds for SPADE as stated in Theorem 1, the dependance on d is only embodied in the variation of test function, instead of the irregularity of points.
We begin to prove Theorem 1. Before that, we need Lemma 1 on the summation property and scaling-invariant property of variation in the sense of Hardy and Krause, 20 as well as Lemma 2 on the combinatorial property of the star discrepancy.

Lemma 1. Let f be defined on the hyper-rectangle
The proof of Lemma 1 can be found in. 20 Combining the Koksma-Hlawka inequality and Lemma 1 yields the generalized Koksma-Hlawka inequality.

Proposition 1 (Generalized Koksma-Hlawka inequality). Let f be of bounded variation on Q = [a, b] in the sense of Hardy and Krause and φ is linear scaling from
Proof. Since the scaling φ in strictly monotone increasing, we have that where the second inequality is the standard Koksma-Hlawka inequality, and the last one uses Lemma 1.
In order to characterize the star discrepancy of a subsequence after removing k points, we need a combinatorial lemma.
Lemma 2 (Combinatorial property of star discrepancy). For any sequence X = (x 1 , . . . , x P ) ⊂ [0, 1] d×P and arbitrary k points X k ⊂ X in it, we have that Proof. Suppose I is the critical box 18 that attains the star discrepancy of the se- Since 0 ≤ A(I, P, X ) − A(I, P − k, X \ X k ) ≤ k, it suffices to take A(I, P, X ) − A(I, P − k, X \ X k ) = m. Then The last inequality utilizes the fact that and m ≤ k. The proof is completed. Now we give a deterministic bound for ξ k in Eq. (3.3), which also bounds the random variables in Section 5.

Lemma 3 (A deterministic bound for random variables). Suppose
Proof. It suffices to consider P k ≥ M k . From Eq. (3.4), a direct calculation yields With these preparations, we are able to finish the proof of Theorem 1.
For the first term, when P k ≥ ϑ √ N, using Lemma 2 yields In this case, we use the fact that

The second term is bounded by
where the second inequality uses Proposition 1 and Lemma 3.
A remarkable observation in the proof of Theorem 1 is that the worse-case upper bound of error function is attained when M k = (P k + ϑ √ N )/2. Actually it accords with our intuition. If M k is sufficiently small, only a few particles are removed and errors should be small. When M k approaches to P k , the bound of discrepancy for Y k also controls the errors. A bad situation may occur when M k is close to P k /2.

Matching and concentration inequalities
The matching of two kinds of particles is the second issue we would like to address. Intuitively, we can pair the particles that are 'sufficiently close to each other' and resort to combinatorial algorithms for solving optimal assignment problems. 14 However, the choice of metrics is rather subtle since the Euclidian distance is highly influenced by dimensionality d and an embarrassing power factor 1/d seems to inevitably occur in error bounds (for instance, see 27 ). Besides, we would like to extend our discussion to a general class of functions, such as that of bounded variation, and drop the assumption of continuity. For these purposes, it requires us to fully utilize the discrete nature of point distributions.
In this section, we would like to show that the random matching, in spite of its simple setting, proves to be surprisedly efficient and accurate. When P k ≥ M k , a random matching from Y k to X k is equivalent to sampling M k times from a finite population X k without replacement, and several powerful concentration inequalities will characterize the non-asymptotic bounds for the random variables ξ k in Eq. (3.3) and their independent sum K k=1 ξ k . Roughly speaking, the concentration inequalities used in our analysis are divided into two categories. One is for sampling with replacement, initialized by the pioneering works of Hoeffding, Bennett and Bernstein, 3,11 which characterizes the concentration behaviors of K k=1 ξ k . The other category is for sampling without replacement, realized first by Serfling, 24 which bounds the random variables ξ k with a martingale structure. The latter has later been improved by Bardenet and Maillard. 2 The Talagrand-type inequality of sampling without replacement was given by Tolstikhin. 29 There are still a vast of researches on concentration inequalities and a good survey can be found in. 5 An exhaustion of them is definitely beyond the scope of this article. Here we only use several typical concentration inequalities, such as those of Hoeffding's type and Bernstein's type. The first part mainly discusses Bernstein's inequality and the Serfling-Hoeffding inequality. The second part concentrates on the Serfling-Bernstein inequality. Again, the discrepancy theory will be used for bounding the local variance in sampling without replacement.

Stochastic bounds under random matching
We first investigate the Serfling-Hoeffding concentration inequalities for sampling without replacement. Without loss of generality, we consider the situation that the number of positive particles is no less than that of negative particles.
Consider the partial summation where X t are randomly sampled from a finite population X = (x 1 , . . . , x P ) without replacement as defined in Eq. (2.1). It is easy to verify that and the two-side tail probabilities U M (δ) and V M (δ) are bounded by the Serfling-Hoeffding inequality, which characterizes the deviation between S k /k and the arithmetic mean of X . Here we only list the concentration bounds in Lemma 4 and the detailed proof can be found in.
The factor M/P in Eq. (5.7) for ρ M reflects the influence of sampling fraction 24 and leads to stronger concentration compared to independent sampling. 29 We would like to use Lemma 4 to bound the random variable 2 . This coincides with our intuition since M = 0 means no particle is annihilated, while M = P corresponds to the case in which all positive particles are selected from the population, so that errors induced by random matching vanish.
The variance estimation of the random matching is also obtained by Lemma 4. In contrast to a trivial bound obtained from Eq. (4.8): we would like to show that the second moment Eξ 2 gains an extra order of N −1/2 with the concentration bound in Eq. (5.5).
Proposition 3 (The first bound for the variance of random matching). Suppose P ≥ M > 0, osc [a,b] (f ) < ∞ and there exists a constant γ < 1 such that P < γ √ N, then Proof. It starts from Owing to the concentration bound (5.5), it further has that Since Eξ = 0, it completes the proof.
In order to characterize the concentration of the independent sum K k=1 ξ k , we need a basic version of Bernstein's inequality.
Proof of Theorem 2. Since Eξ k = 0, it has that the proof of which is the same as that in Eq. (4.9). Combining the above estimate Var(E k (f )).

Variance estimation via bounding discrepancies
We first illustrate concentration bounds of Bernstein's type for the tail probabilities U M (δ) and V M (δ), see Lemma 6, which uses the local variance Its detailed proof is omitted and the interested readers can refer to. 2 Afterwards we will complete the proof of Theorem 3, using the bounds of discrepancies to estimate the local variances and the concentration bounds to estimate the second moment Eξ 2 k .
Lemma 6 (The Bernstein-Serfling inequality). For X = (x 1 , . . . , x P ) ⊂ Q P , Q = [a, b], let (X 1 , . . . , X M ) be a list of a flexible size M < P sampled without replacement from X . Then for any δ > 0 and any ε ∈ (0, 1), we have Moreover, it holds with probability higher than 1 − ε that where ρ M is given in Eq. (5.7) and When osc [a,b] (f ) < ∞, a trivial bound of the local variance is For the points (x 1 , . . . , x P ) with bounded star discrepancy, we could provide another deterministic bound for σ 2 .

Theorem 4 (A deterministic bound for local variance). Suppose
Proof. Since P i=1 f (x i ) = P S P , a direct calculation yields that and then The bounded star discrepancy in the first term gives the first part of the bound in Eq. (5.21). For the second term, we have that which recovers the last two parts of the bound in Eq. (5.21).
Combining Theorem 4 with Lemma 6 produces the second stochastic bound for the random variable ξ in Eq. (5.8).

Proof of Theorem 3. According to Eqs. (3.7), (5.15) and (5.26), the variance Var(E(f )) is bounded by
Var For the first term, we have As a result, it recovers Eq. (3.11).

Numerical validation
To illustrate the negative particle weights, we consider the numerical integration of a determinantal function the prototypes of which are the Slater-determinant-type wave functions in quantum physics and the determinantal point process, where G(v) = (G ij (v)) m×m is a m × m matrix-valued function with G ij (v) being a function of d variables.
Since g(v) is usually not positive semi-definite, the standard Monte Carlo approach is to adopt |g(v)| as the unnormalized instrumental probability density and w(v) = g(v)/|g(v)| ∈ {−1, 1} as the importance weight as follows It yields two kinds of estimators I + (f ) and I − (f ): 23,30 I and where w(x i ) = 1 and w(y i ) = −1. The latter estimator I − (f ) utilizes the strong law of large number Both estimators have the form of Eq. (1.1). In particular, I − (f ) avoids calculating the normalizing constant Z.

Experiment setup
We adopt the convention that (1) A linear function with vanishing mixed derivates: (2) A quadratic function with vanishing mixed derivates: (3) A continuous function with discontinuous derivates: (4) A continuous function with discontinuous derivates: (5) A continuous function with its variation depending on d exponentially:

Implementation
Sampling: The point distributions of (x 1 , . . . , x P ) and (y 1 , . . . , y M ) are obtained by the Markov Chain Monte Carlo (MCMC), where a simple Gaussian random walk Metropolis sampling is adopted to draw from |g(v)| with standard deviation 0.1. 23 Under different settings of dimensionality d and rank m, we need to adjust the parameter ǫ to ensure a high acceptance ratio in the Metropolis sampling and avoid the over-repetitions of samples.
Discrepancy: The calculation of the star discrepancy is a NP-hard problem. 10 Nonetheless, it can be realized by an efficient algorithm based on threshold accepting and probabilistic sampling. 10 Here we set the total number of iterations I = 128 and the trial time T = 5 for each test in order to ensure a satisfactory approximation.
Matching: Two kinds of matching strategies are tested. One is the random matching and the other the optimal assignment. The random matching is realized by sorting a randomly generated sequence in an increasing order. The optimal assignment is realized by the Hungarian algorithm 14 where the distance is chosen as the Manhattan distance. Although the complexity of the Hungarian algorithm scales polynomially in the number of points, its cost still grows dramatically when P k and M k become large. Thus for the sake of comparison, we only use the optimal assignment for the experiments of relatively small sample size (N tot ≤ 10 6 ).
Parallellization: For each group of experiments, we use 64 threads to generate samples from mutually independent Markov chains and collect all data through Message Passing Interference (MPI) standard. In order to realize the sequential clustering in a distributed platform, we decompose the domain into mutually disjoint 64 bins via balancing the particle numbers in each processor, which in turn strikes a balance in overload. All the numerical examples are obtained with our own Fortran implementations and run on the High-Performance Computing Platform of Peking University with the platform: 2×Intel Xeon E5-2697A v4 (2.60GHz, 40MB Cache, 9.6GT/s QPI Speed, 16 Cores, 32 Threads) and 256GB Memory.

Numerical results
To clarify the efficiency of SPADE, we calculate the ratio N A tot /N tot and relative errors for all experiments with ϑ = 0.08 in Table 1, where N tot and N A tot denote the total particle number before and after annihilation, respectively. The partition levels K and the total wall time in sequential clustering, under different N tot , N, ϑ and d, are collected in Table 2 and Table 3, respectively.
For the readers' convenience, all the raw data, ranging from d = 12 to d = 1080, are collected in Tables 4-9 in the appendix as supplementary materials. There Rand(1) corresponds to the results obtained by one random matching in each cluster, Rand(100) the averaged results of 100 random matchings, and Hungarian the results produced by the optimal assignment. The parameter ǫ and the acceptance ratio in sampling, the parameter ϑ in adaptive clustering, the positive particle number P , the negative particle number M and the estimators I − (f ) for all test functions are recorded. The results with relative error exceeding 5% are marked in bold font.
Clustering: We first investigate the performance of adaptive clustering. According to Table 2, the partition level K depends on the parameter ϑ, the total sample size N tot and the dimensionality d. By increasing ϑ, it will dramatically reduce the partition level K so that more particles are annihilated. Too small ϑ is not suggested as it will lead to a very large K, resulting in an inefficient annihilation and a severe over-fitting problem (sometimes K even exceeds the total sample size N tot ). In addition, we observe that for a fixed ϑ, there is a positive correlation between partition level K and total sample size N tot . In particular, when ϑ = 0.08, K is almost linearly proportional to √ N tot . Accuracy and variation: It is readily observed that the accuracy of SPADE is strongly influenced by the variation of test functions. As presented in Table 1, for f 1 , f 2 , f 3 and f 4 , regardless of the continuity of their derivates, SPADE can achieve a satisfactory accuracy even after removing more than 70% of total particles (see the groups d = 36, N tot = 10 7 , ϑ = 0.08 and d = 60, N tot = 10 7 , ϑ = 0.08). And accurate results are still obtained after removing nearly half of particles in very high dimensional problem (see the group d = 360, N tot = 1 × 10 7 , ϑ = 0.08). Therefore, for an appropriate class of test functions, SPADE seems to be immune to the curse of dimensionality. By contrast, a severe fluctuation of errors is observed for f 5 , because its total variation depends exponentially on d. This coincides with the predictions of the deterministic bound (3.5) in Theorem 1 and the random bound (3.8) in Theorem 2, as the performances of E(f ) are distinct for different f even under the same sequences before and after annihilation. Nevertheless, the bound provided by (3.5) and (3.8) sometimes turns out be a gross over-estimate. With an appropriate choice of parameters, SPADE may still be able to preserve the accuracy for test function f 5 in high-dimensional cases after removing more than 10% of total particles (see Table 1 for the groups: d = 360, N tot = 10 5 , ϑ = 0.08, and d = 1080, N tot = 10 6 , ϑ = 0.08).
Matching: We compare the performance of random matching and optimal assignment. With a good partitioning, the accuracy of optimal assignment can be ensured by the deterministic bounds (3.5) in Theorem 1. In general, the optimal assignment is less accurate than the random matching, especially when d becomes large, because the distance in Euclidian space and its geometric structure are heavily influenced by dimensionality. In contrast, the random matching, regardless of once and many times, can improve the accuracy. This validates the improvement in Eq. (3.8) compared to Eq. (3.5) and the concentration bound (3.10) in Theorem 2. In practice, it suggests to do random matching in each cluster for just once to save the computational time, without loss of accuracy. In fact, the result produced by random matching once is sometimes superior to the average of 100 independent trials.
Computational time: We record the total wall time in second for sequential clustering in Table 3, which occupies more than 90% of total computational time. Three observations are listed as follows. (1) For fixed ϑ and N tot , the total wall time almost increases as d increase. The main reason is that more directions are searched in calculating the star discrepancies. 10 (2) For fixed ϑ and d, the total wall time is almost linearly dependent on the sample size N tot , as well as on the partition level K. This observation is rather important in estimating the computational cost of sequential clustering. (3) For fixed N tot and d, the total wall time in general increases as ϑ is chosen smaller, along with a deeper level of the decision tree. But there are some exceptions. One can see that in the group d = 1080, N tot = 1 × 10 7 , the computational cost is the highest when ϑ = 0.08, even though the partition level K is the smallest.

Conclusion and discussion
We proposed an algorithm, dubbed Sequential-clustering Particle Annihilation via Discrepancy Estimation (SPADE), for efficiently removing particles with opposite weights in an empirical signed measure within an acceptable accuracy. SPADE first seeks adaptive clustering of particles via controlling their number-theoretic discrepancies, then pairs the positive and negative particles via random matching and finally removes the paired ones. It does not require any a priori knowledge of nodal hyper-surface of underlying integrand as adopted in fixed-node approximation, and alleviates the restriction of mesh size in grid-based annihilation, thereby providing a potential approach to overcoming the numerical sign problem in general settings. We have proved that both the deterministic and the random error bounds of SPADE are affected by two factors. One factor measures the irregularity of the point distribution via bounding the star discrepancy in each cluster, and is proportional to n −1 (n is the particle number in one cluster). That is, it has the same denominator as the bound of low-discrepancy sequences in quasi-Monte Carlo. Moreover, a unified numerator is set to be ϑ √ N (ϑ ∈ (0, 1) is the error-controlling parameter and N is the normalizing constant), thereby avoiding the embarrassing scaling of (log n) d in quasi-Monte Carlo (d is the dimensionality). The other factor is the variation of the test function, and implicitly depends on d. Therefore SPADE can be immune to the curse of dimensionality for suitable test functions, which has been validated by numerical experiments up to d = 1080. A direct application of SPADE is to alleviate the exponential growth of both particle number and stochastic variances in simulating classical transport and quantum many-body dynamics, especially those where long-time behaviors of several non-oscillating averaged quantities, or physical observables, are of the most interest. In fact, we have recently employed SPADE to alleviate the numerical sign problem in stochastic Wigner simulations. Numerical results in 6-D phase space are very promising 32 and the generalization to 12-D problems is still ongoing.  8 Supplementary: Raw data in numerical experiments Table 4: Numerical results for d = 12, m = 4, ǫ = 0.6. Acceptance ratio in MCMC is 85%. The results with relative error exceeding 5% are marked in bold font.