text:
On the Programmatic Generation of Reproducible Documents

Reproducible document standards, like R Markdown, facilitate the programmatic creation of documents whose content is itself programmatically generated. While these documents are generally not complete in the sense that they will not include prose content, generated by an author to provide context, a narrative, etc., programmatic generation can provide substantial efficiencies for structuring and constructing documents. This paper explores the programmatic generation of reproducible by distinguishing components than can be created by computational means from those requiring human-generated prose, providing guidelines for the generation of these documents, and identifying a use case in clinical trial reporting. These concepts and use case are illustrated through the listdown package for the R programming environment, which is is currently available on the Comprehensive R Archive Network (CRAN).

Background and concepts
R Markdown (Baumer, Cetinkaya-Rundel, Bray, Loi, and Horton 2014) facilitates the construction of computationally reproducible documents by allowing authors to insert R code for data processing, exploration, analysis, table-making, and visualization directly into structured, electronic documents. The resulting documents are made up of these chunks of R code, which we will refer to as computational components since they are generated by computational means, as well as narrative components, which (in scientific writing) is prose intended to contextualize computational components, provide background, define goals, establish themes, and convey the results of the document. The integration of narrative and computational components was identified as "Literate Programming" by Knuth (1984) and software tools, like Sweave (Leisch 2002), have supported this functionality for almost two decades.
R Markdown is particularly popular with its success likely being driven by two factors. The first is the relative ease with which these documents can be constructed. While L A T E X is more expressive, it is relatively technical and requires an investment in time to become proficient. By contrast R Markdown documents are easier to create and format and, when the document is used to create L A T E X, formatting can be passed through to the underlying .tex file. The second factor driving adoption is likely its support for creating modifiable documents, namely Microsoft Word documents. Researchers and analysts, especially those creating applied statistical analyses, often collaborate with domain experts with less technical knowledge. In these cases, the analyst focuses on creating the computational components and narrative components related to results and interpretation. After this initial document is created, the domain expert is free to develop narrative components directly in the document without needing to go through the analyst.
Since computational components are, by definition, computationally derived objects and R Markdown is a well-defined standard, it is possible to programmatically create R Markdown documents with computational components, which is the focus of this paper. Generating documents in this manner has two appealing characteristics. First, it allows us to distinguish the presentation of analytical results from other steps in a data science or data processing pipeline. These other steps including cleaning and analysis often require their own environment and configuration with requirements very different than the computational needs of creating a presentation. By separating these components each can be developed independently. At the same time, by specifying a contract for the output of those objects, we can establish a consistent means by which processed data can be passed to systems for presentating those data in a structured way. The second reason for programmatic creation of R Markdown documents is convenience. In collaborative environments, especially in the early stages, large numbers of graphs and tables are generated and discussed. By collecting these artifacts and structuring them consistently, we can quickly iterate upon and restructure the resulting documents to more clearly present the data without needing to spend time on the creation of the presentation document.
The listdown package provides functions to programmatically create R Markdown files from named lists. It is intended for data analysis pipelines where the presentation of the results is separated from their creation. For this use case, a data processing (or analysis) is performed and the results are provided in a single named list, organized heirarchically. List element names denote sections, subsections, subsubsection, etc. and the list elements contain the data structure to be presented including graphs and tables. The package has native support for workflowr (Blischak, Carbonetto, and Stephens 2019), pdf, word, or html document along with functions allowing a user to easily extend to other types of supported documents. The goal of the package is to create a documents with all tables and visualization that will appear (computational components). This serves as a starting point from which a user can organize outputs, describe a study, discuss results, and provide conclusions (narrative components).
listdown therefore provides a reproducible means for producing a document with specified computational components. It is most compatible with data analysis pipelines where the data format is fixed but the analyses are either being updated, which may affect narrative components including the result discussion and conclusion, or where the experiment is different, which affects all narrative components, but the data format and processing is consistent. An example of the former is provided later in this paper.
One area where we have found listdown is particuarly useful is in the reporting and research of clinical trial data. These collaborations tend to be between (bio)statisticians and clinicians either analyzing past trial data to formulate a new trial or in trial monitoring where trial telemetry (enrollment, responses, etc.) is reported and initial analyses are conveyed to a clinician. The associated presentations require very little context -clinicians often have a better understanding of the data collected than the statistician -often eliminating the need for narrative components. At the same time, a large number of hierarchical, heterogenous artifacts (tables and multiple types of plots) need to be conveyed thereby making the manual creation of R Markdown documents inconvenient.
In this case, data presentation can be fixed across trials. This is especially true in the initial stages, which focus on patient demographics and enrollment. This approach has made it convenient for our group to quickly generate standardize and complete reports for multiple trials concurrently. To date, we have used listdown to report on five clinical trials, with another two currently in process. Results are disseminated using the workflowr package, usually with nine tabs conveying aspects of the data from collection through several different analyses, and each tab containing approxiately five to thirty tables, plots, or other artifacts including trelliscopejs (Hafen and Schloerke 2020) environments which may hold hundreds of graphs. By generating many presentation artifacts we are able to address data-drive questions and issues during collaborative sessions and by carefully structuring these elements we allowing all members to participate in the process.
The listdown package itself is relatively simple with 6 distinct methods that can be easily incorporated into existing analysis pipelines for automatically creating documents that can be used for data exploration and reviewing analysis results as well as a starting point for a more formal write up. These methods include: The rest of this paper is structured as follow. The next section goes over basic usage and commentary. This section is meant to convey the basic approach used by the package and shows how to describe an output document using listdown, create a document, and change how the presentation of computational components can be specialized using listdown decorators. With the user accustomed to the package's basic usage, section 3 describes the design of the package. Section 4 goes over advanced usage of the package including adding initialization code to and outputted document as well as how to control chunk-level options. Section 5 provides a simplified case study of how the package is currently being used in clinical trial reporting. Section 6 concludes the paper with a few final remarks on the general types of applications where listdown has been shown effective.

Basic usage and overview
Suppose we have just completed and analysis and have collected all of the results into a list where the list elements are roughly in the order we would like to present them in a document. It may be noted that this is not always how computational components derived from data analyses are collated. Often individual components are stored in multiple locations on a single machine or across machines. However, it is important to realize that even for analyses on large-scale data, the digital artifacts to be presented are relatively small. Centralizing them makes it easier to access them, since they don't need to be found in mulitple locations. Also, storing them as a list provides a hierarchical structure that translates directly to a document as we will see below.
As a first example, we will consider the a list of visualizations from the anscombe data set. The list is composed of four ggplot2 (Wickham 2016)elements (named Linear, Non Linear, Outlier Vertical, and Outlier Horizontal) each containing a scatter plot from the famous Anscome Quartet -made available in the datasets package (R Core Team 2012). From the computational_components list, we would like to create a document with four sections with names corresponding to the list names, each containing their respective visualizations.

Creating a document with listdown
Creating a document from the computational_components requires two steps. First, we will create a listdown object that specifies how the computational_components object will be loaded into the document, which libraries and code needs to be included, and how the list elements will be presented in the output R markdown document.
R> library(listdown) R> R> ld <-listdown(load_cc_expr = readRDS("comp-comp.rds"), R+ package = "ggplot2") The ld object, along with the computational components in the comp-comp.rds file are sufficient to to create the sections, subsections, and R chunks of a document. The only other thing requires to create the document is the header. The listdown package currently supports regular R Markdown and workflowr as yml objects from the yaml package (Stephens, Simonov, Xie, Dong, Wickham, Horner, reikoch, Beasley, O'Connor, and Warnes 2020). These objects are stored as named lists in R and are easily modified to accomodate document parameters. A complete document can then be written to the console using the code shown below. It could easily be written to a file for rendering using the writeLines() function.
R> doc <-c( R+ as.character(ld_rmarkdown_header("Anscombe's Quartet", R+ author = "Francis Anscombe", R+ date = "1973")), R+ ld_make_chunks (ld) The listdown() function provides document-wide R chunk options for displaying computational components. The chunk options are exactly the same as those in the R markdown document and can be used to tailor the default presentation for a variety of needs. The complete set of options can be found in the R Markdown Reference Guide (RStudio 2020). As a concrete example, the code used to create present the plots could be hidden in the output document using the following code.

Decorators
The first example is simple in part because the ggplot objects both contains the data we want to display and, at the same time, provides the mechanism for presenting them -rendering them in a graph. However, this is not always the case. The objects being stored in the list of computational components may not translate directly to the presentation in a document. In these cases, a function is needed that takes the list component and returns an object to be displayed. For example, suppose that, along with showing graphs from the Anscombe Quartet, we would like to include the data themselves. We could add the data to the computational_components list and then create the document with: R> computational_components$Data <-anscombe R> saveRDS(computational_components, "comp-comp.rds") R> ld_make_chunks(ld)[32:36] [1] "# Data" [2] "" [3] " {r echo = FALSE}" [4] "cc_list$Data" [5] " " In this case, the listdown package will show the entire data set as is the default specified. However, suppose we do not want to show the entire data set in the document. This is common, especially when the data set is large and requires too much vertical space in the outputted document resulting in too much or irrelevant data being shown. Instead, we would like to output to an html document where the data is shown in a datatable thereby controlling the amount of real-estate needed to present the data and, at the same time, providing the user with interactivity to sort and search the data set.
In listdown, a function or method that implements the presentation of a computational component is referred to as a decorator since if follows the classic decorator pattern described in Gamma, Helm, Johnson, and Vlissides (1995). A decorator takes the element that will be presented as an argument and returns an object for presentation in the output directory. A decorator is specified using the decorator parameter of the listdown() function using a named list where the name corresponds to the type and the element correspond to the function or method that will decorate an object of that type. For example, the anscombe data set can be decorated with the DT::datatable() function (Xie, Cheng, and Tan 2020) as: R> ld <-listdown(load_cc_expr = readRDS("comp-comp.rds"), R+ package = c("ggplot2", "DT"), R+ decorator = list(data.frame = datatable)) R> R> ld_make_chunks(ld)[33:37] [1] "# Data" [2] "" [3] " {r}" [4] "datatable(cc_list$Data)" [5] " " List names in the decorator argument provide a key to which a function or method is mapped. The underlying decorator resolution is implemented for a given computational component by going through decorator names sequentially to see if the component inherits from the name using the inherits() function. The function or method is selected from the corresponding name which the element first inherits from. This means that when customizing the presentation of objects that inherit from a common class, the more abstract classes should appear at the end of the list. This will ensure that specialized classes will be encountered first in the resolution process. It should be noted that an object's type is first checked against the decorator name list and then checked to see if it is a list. This allows a user to both decorate a list and retain "list" in its class attributes.
A separate argument, default_decorator, allows the user to specify the default decorator for an object whose type does not appear in the decorator list. This allows the user to specify any class name for the decorator and avoids a potential type name collision with a default decorator whose name is determined by convention. By default, this argument is set to identity but it can be use to not display a computational component by default if the argument is set to NULL.

Design
A listdown object specifies the location of a list of computational components and options for presenting those components in an R Markdown document. The list is a hierarchical data structure that also provides the structure of the outputted document. A corresponding document has two sections "Iris" and "Sepal.Length". The latter has three subsections "Sepal.Width", "Petal.Length", and "Colored". The "Colored" subsection has two sub-subsections, "Sepal.Width" and "Petal.Length". The structure can be seen using the ld_cc_dendro() function. aes(x = Sepal.Length, y = Petal.Length, R+ color = Species)) + geom_point()))) R> R> # Create the dendrogram. R> ld_cc_dendro(comp_comp2) Both the ld_cc_dendro() and ld_make_chunks() functions work by recursively descending the computational components list depth-first. If the list containing and element has a name, it is written to the output as a section, subsection, subsubsection, etc. to a return string. If the visited list element is itself a list, then the same procedure is called on the child list through a recursive call. If the element is not a list, then it is outputted inside an R Markdown chunk in the return string using the appropropriate decorator.

Advanced usage 4.1. Initialization code
The listdown() function facilitates the insertion of initialization code through the init_expr argument. Code is inserted immediately after the libraries are loaded in the R Markdown document. In general, it is suggested that the number of initial expressions be kept small so that the R Markdown document is easy to read. If a large number of functions are required by the target R Markdown document then they can be put into a file and sourced using the initial expression. As an example, suppose we are creating an html document and presenting data using the datatable() function. However, we do not want to include the search capabilities provided. This can be easily accomplished by creating a new function, datatable_no_search(), created using the partial() function Henry and Wickham (2020) to partially apply list(dom = 't') to the options argument of datatable.

R code chunk customization
The listdown package supports also supports capabilities to further customize the presentation by specifying R code chunk options in the R Markdown document in two distinct ways. The first is used when the options we would like to specify is tied to type of the object being presented. This can be though of as a chunk-option decorator. The second is use for changing the options for an individual chunk in an ad-hoc manner.
With three different modes of chunk customization it should be noted that the increasing priority of the chunk options specification is document-wide, decorator-wide, and ad hoc. That is, decorator-wide chunk options take priority over document-wide chunk options and ad-hod options take priority over decorator-wide options. In addition, it should be noted that the use of the lowest priority scheme that accomplishes the presentation goals is preferred because it lends itself to greater code and maintenance efficiency.

Controlling options associated with a decorator
The document-wide chunk option specification provides the default chunk options for output documents generated using listdown. However, the presentation of a data object often varies by type. For example, we may want to specify the height and width of a graph, but not a table. This is accomplished in the listdown package when a listdown object is created using the decorator_chunk_opts option in the listdown() function. For example, associating all ggplot objects with R chunks having a width of 100 and a height of 200 can be accomplished with the following code and it can be seen that only chunk options associated with a plot are modified.

Controlling chunk-level options
Along with providing decorator-wide chunk options, it is also possible to control individual chunk options. The capability is distinct from the document-wide and decorator-wide specification of options in that it must be applied to the computational component list element whose associated options will be modified, rather than the listdown object. This is because the listdown only specifies how classes of objects should be presented. To modify the chunk options associated with a specific list element the list element is provided with a set of attributes that can be queried by the ld_chunk_opts() function as the output document is being generated. Because of the ad hoc nature of this capability, its use is discouraged. A better solution, that maintains the behavior is to add class information to the list element and specify decorator-wide chunk options for the new class. This maintains the separation of the computational component list, which maintains the document structure and data for presentation from the specification of how the document will be created and rendered.
R> comp_comp2$Iris <-ld_chunk_opts(comp_comp2$Iris, echo = TRUE) R> saveRDS(comp_comp2, "comp-comp2.rds") R> ld_make_chunks(ld)[12:16] [1] "# Iris" [2] "" [3] " {r echo = TRUE}" [4] "cc_list$Iris" [5] " " 5. A simple example: Reporting on the gtsummary::trial dataset As mentioned before, we have found the listdown package particularly helpful for reporting the results of clinical trials thereby creating a basis for discussion and collaboration between (bio)statisticians and the clinicians running the trials. For this use case, the context and data collection procedures are well-understood and as a result very few narrative components are needed. It is also the case that the modes of presentation (tables, scatterplots, survivalplots, consort diagrams, etc.) are standardized. The goal in presenting the trial characteristics is to identify problems in the data, monitor trial enrollment and response, quantify known relationships among the data, and test hypotheses about a therapy's efficacy.
In practice we generally separate the data cleaning, exploration, analysis, monitoring, and presentation components. The computational component list has tens of elements with thousands of visualizations. These facts, coupled with the privacy constraints make a complete example difficult for the purposes of this paper. So, below we provide a simple example of the types of documents we provide using the gtsummary::trial dataset. While it is not complete, it does convey they types of reports we are currently generating for completed and ongoing clinical trials.

Conclusion
While the programmtic generation of reproducible documents has appealing qualities, it also has fundamental limitations that should be kept in mind when deciding if tools like listdown should be employed. First and foremost, without narrative components a document has very little context. Quantitative analyses require research questions, hypotheses, reviews, interpretations, and conclusions. Computational components are necessary but generally not sufficient for constructing an analysis. This means that if narrative components must be conveyed in a document, then listdown may make the generation of their presentation more convenient. However, it does not relieve the burden of the author to create prose developing a narrative.
Second, it is difficult if not impossible to construct computational components for an arbitary analyses. Analyses themselves have context and are built with a set of assumptions and goals. Our experience shows that listdown is easiest to used for a fixed data format. This means a standard set of table and visualizations for similarly formatted data. In particular, data that is periodicaly updated, without changes to the format, is a case that is particularly amenable to document generation reuse.
Keeping these limitations in mind, listdown can be used to effectively reduce the the difficulty of generating documents in a variety of contexts and fits readily in data processing and analysis pipelines.


extfieldsofstudy:
Physics
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
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.


Improvement of genetic algorithm using artificial bee colony

Received Nov 13, 2019 Revised Jan 4, 2020 Accepted Feb 8, 2020 Genetic algorithm (GA) is a part of evolutionary computing that simulates the theory of evolution and natural selection, where this technique depends on a heuristic random search. This algorithm reflects the operation of natural selection, where the fittest individuals are chosen for reproduction so that they produce offspring of the next generation. This paper proposes a method to improve GA using artificial bee colony (GABC). This proposed algorithm was applied to random number generation (RNG), and travelling salesman problem (TSP). The proposed method used to generate initial populations for GA rather than the random generation that used in traditional GA. The results of testing on RNG show that the proposed GABC was better than traditional GA in the mean iteration and the execution time. The results of testing TSP show the superiority of GABC on the traditional GA. The superiority of the GABC is clear in terms of the percentage of error rate, the average length route, and obtaining the shortest route. The programming language Python3 was used in programming the proposed methods.

INTRODUCTION
Artificial intelligence (AI) has been shown great performance in several tasks [1][2][3][4]. Genetic algorithm (GA) is an AI algorithm that is a metaheuristic method invented in 1975 by John Holland [5]. GA has a good global search capability through its operation strategy but has a slow convergence speed [6]. GA shows high performance in complex areas because it works repetitively as the research is likely to focus on typical representative area areas found to produce good behavior [7]. Artificial bee colony (ABC) algorithm was proposed by Karaboga in 2005, which is simple in concept and requires very few initialization parameters to adjust. It has properties of fast convergence speed, good quality of solutions and good robustness, etc., but it also has the disadvantages of early convergence and ease of falling into the local optimum [8].
Important problems that test the efficiency of the algorithm are RNG and TSP, which will be adopted in this thesis as a case study, which will be explained later. Random numbers are numbers generated by an operation that unpredictable it's an outcome and cannot be sequentially reliably reproduced. Random numbers are necessary for a variety of applications. For instance, a popular cryptosystem used keys that should be produced in a random mode [9]. TSP is one of the problems known as NP-hard and is a wellknown problem in the field of computer science and mathematics. This problem has been used in many fields such as semiconductor manufacturing, logistics, and transportation [10].
The rest of the paper is organized as follows: section 2 presented some related work. In section 3, background theories were displayed. In section 4, described the design of the proposed methods. In section 5, experimental results were evaluated, and in the last section, a brief conclusion was presented.

RELATED WORK
There are many pieces of research focusing on improving GA. In Abu-Almash [9] genetic algorithm used to propose a new approach to generate keys. The Five basic tests applied as a fitness function in this experiment for each chromosome. The results proved that the competence of a genetic algorithm for key generation. The optimal representation for the genetic algorithm may be seen by the initial population, size of population, crossover, and mutation. The large population size, little mutation probability, and a high crossover probability.
Alkafaween [11] has presented a new mutation called "IRGIBNNM", that was a hybrid of two mutations which were a random-based mutation and a knowledge-based mutation. The proposed mutation was used to improve the select best mutation (SBM) strategy. The proposed mutation performance was compared with three mutations, additionally to the SBM. The performance of the proposed mutation and the SBM strategy was evaluated on several experiments on twelve TSP instances were conducted. Those instances were taken from the TSPLIB. The experiential results showed the efficiency of the "IRGIBNNM" mutation and the SBM strength where both methods have benefited from the randomization and knowledge provided by the nearest neighbor approach.
Hassanat et al. [12], studied GA was improved by dividing the large-scale TSP problem into small subsets based on regression. This was achieved by using the regression line and its perpendicular line, where cities were grouped into four sub-problems and frequently, each site locates the group in which the city belongs. This process was repeated until the size of the problem became very small or reached a size that cannot be divided. Some famous TSP instances that are derived from TSPLIB were used in experiments. The experiments showed that the proposed method outperformed the other methods using population seeding techniques such as the nearest neighbor based techniques and random techniques regarding mean convergence and error rate.
Alkharji et al. [13] describe a method for using GA to generate random keys for a fully homomorphic encryption scheme (FHE) and check its efficiency. The five statistical tests used in this experiment. This paper explains that the powerful, high-quality, non-repetitive random keys generated by GA will increase the security of FEH schemes, and making it difficult for analysts to break the data. Akter et al. [14] have presented a new crossover operator for solving the TSP. The new crossover consisted of choosing two crossover points and creating new offspring by comparing the cost. The proposed crossover was performed with two traditional crossover operators a sequential crossover operator (SCX) and a triple crossover operator (TCO), and used eight benchmark instances taken from the TSPLIB. The experimental results among the proposed crossover operator, SCX, and TCO showed that the proposed operator has outperformed other operators via providing a solution in less iteration that significantly reducing time and memory.
Gupta et al. [10] have focused on the development of a heuristic technique by combining two common optimization methods: particle swarm (PSO) and GA for TSP. The main inspiration was used to improve a hybrid GA-PSO algorithm to take advantage of PSO such as the self-improvement of individuals and high convergence rate, which compensated the weakness of GA. In the GA-PSO hybrid algorithm, new individuals were generated through the PSO mechanism as well as crossover and mutation operators. The performance of the GA-PSO hybrid algorithm against GA and PSO was used ten TSP standards in terms of finding optimal and execution time. Superior GA-PSO mixed performance between GA and PSO for TSP with concern to the standard, i.e. error rate and average computational time.

PRELIMINARIES 3.1. Genetic algorithm
GA is a non-linear, discrete random process that does not require mathematical formulation. Optima are evolving from generation to another [15]. It has grown as important in solving large complex problems that require a long time according to traditional programming techniques, compared to GA, which has a huge amount of alternative solutions, where the solution is often optimal or near to optimal, over time [16].
At first, GA selects parents from an initial population of chromosomes, and then generates offspring using crossover and mutation. It evaluates all chromosomes to determine their fitness using a fitness function. Then, fitness values were used to determine if chromosomes are retained or disposed of. The least adapted chromosomes are eliminated and the most adapted chromosomes are retained in generating new populations according to the principle of survival of the fittest. The new population replaces the old. The processes are repeated while a specific termination condition satisfied [17]. A simple GA is illustrated in Figure 1.

Artificial bee colony (ABC)
ABC is a method of improving heuristic population-based. The populations in the algorithm represent bees. Every bee searches for the best solution (food resource). It is assumed that each solution is a food resource and the amount of nectar for each resource represents the quality of each solution [18]. ABC is a population-based improvement algorithm that aims to obtain a global minimum [19].
There are two types of bees in the ABC algorithms that are employed and unemployed bees that be consisting of the onlooker bees and the scout bees. The population of the colony usually composed two parts employed that form the first half and the rest include onlooker bees. At first, employed bees exploit the sources of nectar that have been explored, and then provide information to the onlooker bees in the beehive about the location of the food sources they exploit.
After completing the search process by the employed bees, they share information about food sources and their location with the onlooker bees in the dance space. The onlooker bees evaluate information from employed bees and choose the source of food with the probability of nectar. Scouts randomly examine the environment to identify a new food source based on potential external evidence or internal motivation. Algorithm 1 represents the main algorithm of ABC.

Travelling salesman problem
TSP is considered an optimization problem [20]. It is a problem that is an easy description tricky solution and classified as a problem that has not been solved in polynomial time. This problem belongs to NP-hard [21]. The major purpose of this problem is to detect the shortest route through a group of cities (starting from city N and ending in the same city) so that each city is visited only once [22].

Random number generators
Random numbers are a wide field with a strong theoretical and empirical basis. There are many methods used to generate random numbers such as published random number tables and throwing dice. The expansion of digital computers has opened up a new area of inquiry, such as measurement of physical parameters, random numbers of real generators, etc [23]. A valid RNG will generate a sequential series of distributed numbers, at the interval [0, 1]. Using any type of computer algorithm, if we know the initial state of the algorithm which be also called the seed, the sequence of numbers can be determined [24].

Statistical tests
This section presents several tests aimed at measuring the quality of a generator called a random bit generator. The process is done by taking a sampling of the sequence and presenting it to several statistical tests. Every one of the five tests determined whether a particular attribute sequence is likely to actually display a random sequence [25]. Let s=s 0 , s 1 , s 2 …s n−1 is a binary sequence of lengths n. The statistical tests that are used in this paper are [26]: -The frequency (mono-bit) test For a random series of length n, the number 0's should be approximately equal to the number of 1's. The statistic test used in (1): (1) where: s 0 , s 1 indicates the numbers of 0's and 1's in s, respectively.
-Poker test Suppose p is a positive integer where [s/p]≥5, (2 p ), and suppose D=[s/p]. Divide the sequence s into non-overlapping segments D whose length is p, and assume that s i is the number of sequence-type appearances of length p, 1≤I≤2p. In this test, we determine whether the length p sequence shows approximately the same number of times in seconds, as expected for the random sequence. The statistic test used in (3): -Run test The purpose of this test is to determine whether the number is either zeros or numbers from a series of different lengths in the S sequence as expected in a random sequence. The expected number of gaps (or blocks) of length i in a random sequence of length s is ⁄ . Suppose L is equals to the maximum integer i was e i ≥5. Suppose B i and G i are the numbers of block and gap respectively, of length i in n for each i, 1≤i ≤L. The statistic test is calculated by (4):  (5):

THE PROPOSED METHOD
The proposed method has employed an ABC algorithm to improve the population of GA. TSP and RNG are used as a case study. For RNG problem the five statistical tests were used to compute the fitness (each time any test was failed the counter of fitness incremented), such that the random numbers with minimum fitness was token, and the termination criterion of the proposed algorithm has reached the maximum cycle number or all the population is random numbers passed the tests. In TSP, the fitness function in this work was the minimum distance between cities. Total distance for each TSP route for cities n was calculated according to the coordinates of the cities as follows: If the coordinates of any two cities i and k are (xi, yi) and (xk, yk), the distance between them was calculated according to the Euclidean distance equation: Total tour length f can be expressed as follows: where n is the total number of cities. The termination criterion of the proposed algorithm has reached a maximum cycle number.

Improving GA population using ABC for RNG
This method consists of two stages where the first stage focuses on generating random numbers using ABC while the second stage used GA on the random numbers that are generated in the first stage. Previous hybridization work was improved by generating an initial population by ABC then GA was applied to the initial population. At last, the population produced from GA was compared with the initial population and selects the best population (replaces the original parental population) as the new population to the next generation. Reproduction imitated natural selection. ABC generated numbers randomly within the boundaries of the parameters by using the following equation: (8) where; i=1. . .SN, j=1. . .D. Where SN is the solution number and D is the number of optimization parameters The created number was converted to binary. After this, the fitness function in the proposed algorithm was that the five statistical tests might be assigned to the X ij solution by (9), in Algorithm 2.
Fitnessi ∑ (9) where if it the counter value of the solution X ij . Step1: Generate random numbers using Eq. 8 to be the initial population (I). /* using ABC*/ Step2: Evaluate (I) using the five statistical tests to compute their fitness using Eq.9.

Repeat
Step3: Call algorithm (3) that takes I and produce a new solution (S) using GA.
Step4: Evaluate S using the five statistical tests to compute their fitness using Eq. 9.
Step5: Select the best populations between S and I to be a new population (K) according to their fitness. /*Reproduction operation*/ Step6: Until the termination condition is met. /* reached the maximum cycle number or all the population is random numbers passed the tests*/ Step7: Return K. End. Based on the initial population that was generated from the ABC stage GA performed genetic operators to generate new offspring. The genetic process was performed iteratively. The new population was sorted according to its fitness value; the roulette wheel operation was used in the proposed algorithm. An ordered crossover and uniform mutation were used in Algorithm 3. Step1: Sort I according to their fitness value. /* ascending order */ Step2: Select from I the best solutions depending on their fitness value (BS). /* Using roulette wheel selection method*/ Step3: Perform crossover operation between BS to produce new solutions (NS). /* ordered crossover is used, with a crossover rate of 0.8*/ Step 4: If there is a mutation rate, perform mutation operation and update NS. /* uniform mutation was used, with a mutation rate of 0.1*/ Step5: Return NS. End.

Improving GA population using ABC for TSP
This method consists of two stages in the first stage ABC generates a new route by using (8), and in the second stage, the population generated in the first stage will be using GA operators to find the best route. ABC generated random routes by using the index of cities in the default route within the boundaries of the indices by using (8). After this, (7) was used to calculate the distance for each route. The proposed method to find the shortest path is shown in Algorithm 4. Step2: Evaluate the new routes using equation (7). Repeat Step3: Call algorithm (5) that takes NR and produces NS using GA.
At the GA stage, the population produced from the ABC stage was sorted depending on their fitness value. Roulette wheel operation ordered crossover and swap mutation were used in the proposed algorithm to generate a new gene. The GA stage is shown in Algorithm 5. Step1: Sort the NR according to their fitness value. /* ascending order */ Step2: Select from NR the best solutions depending on their fitness value (BS). /* Using roulette wheel selection method*/ Step3: Perform crossover operation between BS to produce new solutions (NS). /* ordered Crossover is used, with a crossover rate of 0.85 */ Step4: If there is a mutation rate, perform the mutation operation/*swap mutation, with a mutation rate of 0.01*/ Step5: Return NS.

EXPERIMENTAL RESULTS AND DISCUSSION
The proposed method was evaluated in this section. Two experiments were carried out where the first experiment aimed to solve RNG while the second problem aimed to solve TSP. Then the solutions were compared with the traditional GA. All computational experiments were conducted on Intel Core i7-4600U 2.70 GHz machine and coded using Python 3 programming language.

RNG experiment parameters
The proposed method tested on data comprises 64, 120,128,192,256 and 512 keys' length respectively. The process was executed 5 times with various recommended parameters for the three variations of data. The population size was 10, 20, and 30, respectively, and the maximum iteration number

2131
(MIN) was 100. The rate of the crossover was 0.8 and the rate of mutation was 0.1. The goal was that all keys were random according to the five statistical tests. It is based on the number of iterations to achieve the best solution (7), and the average time execution.
Mean= ∑ Where the rate is the number of iterations to achieve stop condition, and n is the number of run times. The summary of the simulation is presented in Tables 1-3. Through the tables below, the proposed algorithm showed better results than the traditional GA, as generated random numbers in the least number of iterations and at the lowest execution time, which is what was required to get better solutions in the least time.

TSP experiment parameters
Experiments were performed to calculate the percentage of relative error, the best tour, and the average of the tour. The rate of crossover and rate of mutation was 0.85 and 0.01 respectively. The population size equaled the number of cities, and the number of generations was 5000. The simulation was performed 5 times with various recommended parameters. The percentage of relative error (%) was calculated using (11).
These algorithms were tested using five real TSP problems taken from the TSPLIB which include eil51, st70, pr76, eil76, and rd100 (the numbers attached to the problem names represent the number of cities). The default tour is taken from the library of TSPLIB for each instance. Through a Table 4, the proposed algorithm showed better results than the traditional GA in all five instances, which was what is required to get a speed convergence to an optimal solution in the least time. When comparing the proposed algorithm with the GA-PSO algorithm [6] and IRGIBNNM algorithm [7] it is showed that the results obtained from the proposed algorithm were better than the previous two algorithms. Table 5 showed the comparison among the algorithms.

CONCLUSION
In this paper, GA was improved using ABC (GABC) by using ABC to generate an initial population rather than randomly generation that used in traditional GA. Through experimental results ,   it was concluded that it using a hybrid between GA with ABC has improved GA in RNG by minimizing the number of iterations in less execution time. And using the five statistical tests has provided a good fitness function for the RNG problem. The final generated keys that were unique (not repeated), random and cryptographically strong (successfully passed the five statistical tests). The result of the proposed method in TSP was compared with traditional GA based on the percentage of relative error rate and the average of the tour's average. Results showed GABC was better than traditional GA with a lower error rate and high convergence rate.