{'fraction_non_alphanumeric': 0.03533042651307172, 'fraction_numerical': 0.012201467764798765, 'mean_word_length': 4.574872461738521, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Particle algorithms for optimization on binary spaces

We discuss a unified approach to stochastic optimization of pseudo-Boolean objective functions based on particle methods, including the cross-entropy method and simulated annealing as special cases. We point out the need for auxiliary sampling distributions, that is parametric families on binary spaces, which are able to reproduce complex dependency structures, and illustrate their usefulness in our numerical experiments. We provide numerical evidence that particle-driven optimization algorithms based on parametric families yield superior results on strongly multi-modal optimization problems while local search heuristics outperform them on easier problems.

using sequential Monte Carlo techniques. In the sequel, we refer to f as the objective function. For an excellent overview of applications of binary programming and equivalent problems we refer to the survey paper by Boros and Hammer (2002) and references therein. The idea to use particle filters for global optimization is not new [(Del Moral et al., 2006), Section 2.3.1.c], but novel sequential Monte Carlo methodology making use of suitable parametric families on binary spaces Schäfer and Chopin (2011) may allow to construct more efficient samplers for the special case of pseudo-Boolean optimization.
We particularly discuss how this methodology connects with the cross-entropy method Rubinstein (1997) which is another particle driven optimization algorithm based on parametric families.
The sequential Monte Carlo algorithm as developed by Schäfer and Chopin (2011) is rather complex compared to local search algorithms such as simulated annealing Kirkpatrick et al. (1983) or k-opt local search Merz and Freisleben (2002) which can be implemented in a few lines. The aim of this paper is to motivate the use of advanced particle methods and sophisticated parametric families in the context of pseudo-Boolean optimization and to provide conclusive numerical evidence that these complicated algorithms can indeed outperform simple heuristics if the objective function has poorly connected strong local maxima. This is not at all clear, since, in terms of computational time, multiple randomized restarts of fast local search heuristics might very well be more efficient than comparatively complex particle approaches.

Outline
The article is structured as follows. We first introduce some notation and review how to model an optimization problem (1) as a filtering problem on an auxiliary sequence of probability distributions. Section 2 describes a sequential Monte Carlo sampler Del  designed for global optimization on binary spaces Schäfer and Chopin (2011). Section 3 reviews three parametric families for sampling multivariate binary data which can be incorporated in the proposed class of particle algorithms. Section 4 discusses how the cross-entropy method Rubinstein (1997) and simulated annealing Kirkpatrick et al. (1983) can be interpreted as special cases of the sequential Monte Carlo sampler. In Section 5 we carry out numerical experiments on instances of the unconstrained quadratic binary optimization problem. First, we investigate the performance of the proposed parametric families in particle-driven optimization algorithms. Secondly, we compare variants of the sequential Monte Carlo algorithm, the cross-entropy method, simulated annealing and simple multiple-restart local search to analyze their respective efficiency in the presence or absence of strong local maxima.

Notation
We briefly introduce some notation that might be non-standard. We denote scalars in italic type, vectors in italic bold type and matrices in straight bold type. Given a set M , we write |M | for the number of its elements and 1 M for its indicator function. For a, b ∈ Z we denote by [[a, b]] = {a, . . . , b} the discrete interval from a to b. Given a vector x ∈ B d and an index set I ⊆ [[1, d]], we write x I ∈ B |I| for the sub-vector indexed by I and x −I ∈ B d−|I| for its complement. We occasionally use the norms x ∞ := max i x i and |x| := d i=1 |x i |.

Statistical modeling 1.2.1 Associated probability measures
For particle optimization, the common approach is defining a family of probability measures {π : ≥ 0} associated to the optimization problem max x∈B d f (x) in the sense where U S denotes the uniform distribution on the set S and M f = argmax x∈B d f (x) the set of maximizers. The idea behind this approach is to first sample from a simple distribution, potentially learn about the characteristics of the associated family and smoothly move towards distributions with more mass concentrated in the maxima. We review two well-known techniques to explicitly construct such a family π .
As increases, the modes of π become more accentuated until, in the limit, all mass is concentrated on the set of maximizers. The name reflects the physical interpretation of π (x) as the probability of a configuration x ∈ B d for an inverse temperature and energy function −f . This is the sequence used in simulated annealing Kirkpatrick et al. (1983).
Level set family We call {π : ≥ 0} a level set family, if it has probability mass functions of the form π (γ) := |L + | −1 1 L + (γ), where Indeed, L + is the super-level set of f with respect to the level c = f (x * ) − 1/ , for > 0, and π (γ) is the uniform distribution on L + . As increases, the support of π becomes restricted to the points that have an objective value sufficiently close to the maximum of the f . In the limit, the support is reduced to the set of global maximizers. The particle-driven optimization algorithms are computationally more involved than local search heuristics since we need to construct a sequence of distributions instead of a sequence of states. We shall see that this effort pays off in strongly multi-modal scenarios, where even sophisticated local search heuristics can get trapped in a subset of the state space.

Rare event simulation
While the tempered sequence is based on a physical intuition, the level set sequence has an immediate interpretation as a sequence of rare events since, as increases, the super-level set becomes a 'rare event' with respect to the uniform measure. Rare event simulation and global optimization are therefore closely related concepts and methods for rare event estimation can often be adapted to serve as optimization algorithms.
Particle algorithms for rare event simulation include the cross-entropy method Rubinstein (1997) and the sequential Monte Carlo sampler Johansen et al. (2006). The former uses the level set sequence, the latter uses a logistic potential family ) and : R → (0, 1), (x) = [1 + exp(−x)] −1 denotes the logistic function. Johansen et al. (2006) did not specifically design their algorithm for optimization but their approach to static rare event simulation is closely related to the particle optimization framework.

Sequential Monte Carlo
We discuss a static sequential Monte Carlo sampler that uses a transition kernel with independent proposals in the move step based on suitable parametric families. This methodology has been demonstrated to reliably estimate the mean of the posterior distribution in Bayesian variable selection problems for linear normal models Schäfer and Chopin (2011). In this section, we provide a self-contained description of this framework. For a more general overview of sequential Monte Carlo methods we refer to Del .

Sequential Importance Sampling
For convenience of notation, we index the sequence of distributions (π t ) t directly by t rather that by the parameter of the family t . We refer to X = (x 1 , . . . , x n ) ∈ B n×d and w ∈ [0, 1] n with |w| = 1 as a particle system with n particles. We say the particle system (w, X) targets the probability distribution π if the empirical distribution n k=1 w k δ x k converges to π for n → ∞.

Importance weights
We sample (w 0 , X 0 ) with x 1,0 , . . . , x n,0 ∼ U B d = π 0 and set w 0 = 1/n to initialize the system. Suppose we are given a particle approximation (w t , X t ) of π t and want to target the subsequent distribution π t+1 . For all k ∈ [[1, n]] and α > 0, we update the weights to andπ ∝ π denotes the unnormalized version of π . The normalizing constants ν and |L + | defined in equations (2) and (3) are unknown but the algorithm only requires ratios of unnormalized probability mass functions. We refer to α as the step length at time t. After updating, the particle system targets the distribution π t+α ≈ n k=1 w k,t,α δ x t,k .
As we choose α larger, that is π t+α further from π t , the weights become more uneven and the accuracy of the importance approximation deteriorates. If we repeat the weighting step, we just increase α and finally obtain an importance sampling estimate of π ∞ = U M f with instrumental distribution π 0 = U B d . This yields a poor estimator since the probability to hit the set M f with n uniform draws is 1 − (1 − 2 −d |M f |) n and decreases rapidly as the dimension d grows. The pivotal idea behind sequential Monte Carlo is to alternate moderate updates of the importance weights and improvements of the particle system via resampling and Markov transitions.
Effective sample size The importance weight degeneracy is often measured through the effective sample size criterion Kong et al. (1994) defined as The effective sample size is 1 if the weights are uniform, that is equal to 1/n; the effective sample size is 1/n if all mass is concentrated in a single particle.

Finding the step length
Given any increasing sequence (π t ) t , we could repeatedly reweight and monitor whether the effective sample size falls below a critical threshold. For the special case of annealing via sequential Monte Carlo, however, the effective sample size after weighting η n (w t,α ) is merely a function of α. For a particle system (X t , w t ) at time t, we pick a step length α such that η n (w t ) β = η n (w t,α ), that is we lower the effective sample with respect to the current particle approximation by some fixed ratio β ∈ (0, 1) [ ), (Del Moral et al., 2011]. This ensures a 'smooth' transition between two auxiliary distributions, in the sense that consecutive distributions are close enough to approximate each other reasonably well using importance weights; in our numerical experiments, we took β = 9/10. We obtain the associated sequence ( t ) t by setting t+1 = t +α t where α t is a unique solution of (6). Since η n (w t,α ) is continuous and monotonously decreasing in α we can use bi-sectional search to solve (6).

Conditioning the particle system
In this section, we discuss how to condition the weighted system and improve its quality before we proceed with the next weighting step.

Resampling
We replace the system (w t+1 , X t ) targeting π t+1 by a selection of particlesx 1 , . . . ,x n drawn from the current particle reservoir x 1,t , . . . , x n,t such that where n(x) denotes the number of particles identical with x. Thus, in the resampled system particles with small weights have vanished while particles with large weights have been multiplied. For the implementation of the resampling step, there exist several recipes. We could apply a multinomial resampling Gordon et al. (1993) which is straightforward. There are, however, more efficient ways like residual Liu and Chen (1998), stratified Kitagawa (1996) and systematic resampling Carpenter et al. (1999). We use the latest in our simulations, see Procedure 1.

Moving the system
If we repeated the weighting and resampling steps several times, we would rapidly reduce the number of different particles to a very few. The key to fighting the depletion of the particle reservoir is moving the particles according to a Markov transition kernel κ t+1 with invariant measure π t+1 . The particlex (0) k,t+1 is by construction approximately distributed according to π t+1 , and a draŵ is therefore again approximately distributed according to π t+1 . The last sample of the generated Markov chain (x (0) k,t+1 , . . . ,x (s) k,t+1 ) is, for sufficiently many move steps s ∈ N, almost exactly distributed according to the invariant measure π t+1 and independent of its starting point.

Stopping rule
While we could always apply a fixed number of move steps, we rather use an adaptive stopping criterion based on the number of distinct particles.
Particle diversity We define the particle diversity as Ideally, the sample diversity ζ n (X) should correspond to the expected diversity ζ n (π) := 1 ∧ n −1 γ∈B d 1 {x∈B d : cnπ(x)≥1} (γ), where c n is the smallest value that solves γ∈B d c n π(γ) ≥ n. This is the particle diversity we would expect if we had an independent sample from π(x). Therefore, if κ t+1 is fast-mixing, we want to move the system until Since the quantity on the right hand side is unknown, we stop moving the system as soon as the particle diversity reaches a steady state we cannot push it beyond Schäfer and Chopin (2011).
More precisely, we stop if the absolute diversity is above a certain threshold ζ * ≈ 0.95 or the last improvement of the diversity is below a certain threshold ζ * ∆ > 0. We always stop after a finite number of steps but the thresholds ζ * and ζ * ∆ need to be calibrated to the efficiency of the transition kernel. For slow-mixing kernels, we recommend to perform batches of consecutive move steps instead of single move steps.
If the average acceptance rate λ of the kernel (see Section 2.2.4) is smaller than ζ * ∆ , it is likely that the algorithm stops after the first iteration although further moves would have been necessary. We could adaptively adjust the threshold ζ * ∆ to be proportional to an estimate of the average acceptance rate; for our numerical experiments, however, we kept it fixed to ζ * ∆ ≈ 10 −2 .

Transition kernels
Most transition kernels in Monte Carlo simulations are some variant of the Metropolis-Hastings kernel (see e.g. Robert and Casella (2004)), where we sample from the kernel by proposing a new state γ ∼ q t+1 (γ | x) and accepting the proposal with probability or returning x otherwise. Again, we denote byπ t+1 ∝ π t+1 the unnormalized version of π t+1 since the kernel only requires the ratio of the unnormalized probability mass functions.
Symmetric kernel On binary spaces, a common choice for the proposal distribution is with weight vector p ∈ [0, 1] d normalized such that |p| = 1.
With probability p k , the kernel proposes a uniform draw from the k-neighborhood of x, We refer to this type of kernel as symmetric kernel since q(γ | x) = q(x | γ) and equation (7) simplifies. This class of kernels provide a higher mutation rate than the random-scan Gibbs kernel (see Schäfer and Chopin (2011) for adiscussion).
Locally operating transition kernels of the symmetric type are known to be slowly mixing. If we put most weight on small values of k, the kernel only changes one or a few entries in each step. If we put more weight on larger values of k, the proposals will hardly ever be accepted if the invariant distribution π is multi-modal. Ideally, we want the particles sampled from the transition kernel to be nearly independent after a few move steps which is often hard to achieve using local transition kernels.
Adaptive independent kernel For the sequential Monte Carlo algorithm, we use adaptive independent kernels which have proposal distributions of the kind which do not depend on the current state x but have a parameter θ which we adapt during the course of the algorithm.
The adaptive independent kernel is rapidly mixing if we can fit the parametric family q θ such that the proposal distribution q t+1 = q θ t+1 is sufficiently close to the target distribution π t+1 , yielding thus, on average, high acceptance rates λ q t+1 . The general idea behind this approach is to take the information gathered in the current particle approximation into account (see e.g. Chopin (2002)). The usefulness of this strategy for sampling on binary spaces has been illustrated by Schäfer and Chopin (2011).
We fit a parameter θ t+1 to the particle approximation of π t+1 according to some suitable criterion. Precisely, θ t+1 is taken to be the maximum likelihood or method of moments estimator applied to the weighted sample (w t+1 , X t ). The choice of the parametric family q θ is crucial to the implementation of a sequential Monte Carlo sampler with adaptive independent kernel. We discuss this issue in detail in Section 3.
Adaptation could, to a certain extent, also be done for local transition kernels. Nott and Kohn (2005) propose an adaptive kernel which replaces the full conditional distribution of the Gibbs sampler by an easy to compute linear approximation which is estimated from the sampled particles. This method accelerates Gibbs sampling if the target distribution π is hard to evaluate but does not provide fast mixing like the adaptive independent kernel (see Schäfer and Chopin (2011) for a comparison).
Still, the use of local kernels in the context of the proposed sequential Monte Carlo algorithm might be favorable if, for instance, the structure of the problem allows to rapidly compute the acceptance probabilities of local moves. Further, batches of local moves can be alternated with independent proposals to ensure that the algorithm explores the neighborhood of local modes sufficiently well.

Remark on discrete state spaces
Since the sample space B d is discrete, a given particle is not necessarily unique. This raises the question whether it is sensible to store multiple copies of the same weighted particle in our system. In the sequel, we discuss some more details concerning this issue which has only been touched upon briefly by Schäfer and Chopin (2011).
Let n(x) denote the number of copies of the particle x in the system (w, X). Indeed, for parsimonious reasons, we could just keep a single representative of x and aggregate the associated weights to w * (x) = n(x) w(x).

Impact on the effective sample size
Shifting weights between identical particles does not affect the nature of the particle approximation but it obviously changes the effective sample size η n (w) which is undesirable since we introduced the effective sample size as a criterion to measure the goodness of a particle approximation. From an aggregated particle system, we cannot distinguish the weight disparity induced by reweighting according to the importance function (4) and the weight disparity induced by multiple sampling of the same states which occurs if the mass of the target distribution is concentrated. More precisely, we cannot tell whether the effective sample size is actually due to the gap between π t and π t+1 or the presence of particle copies due to the mass of π t concentrating on a small proportion of the state space which occurs by construction of the auxiliary distribution in Section 1.2.

Impact on the resample-move step
Aggregating the weights means that the number of particles is not fixed at runtime. In this case, the straightforward way to implement the move step presented in Section 2.2.2 is breaking up the particles into multiple copies corresponding to their weights and moving them separately. But instead of permanently splitting and pooling the weights it seems more efficient to just keep the multiple copies.
We could, however, design a different kind of resample-move algorithm which first augments the number of particles in the move step and then resamples exactly n weighted particles from this extended system using a variant of the resampling procedure proposed by Fearnhead and Clifford (2003). A simple way to augment the number of particles is sampling and reweighting via k ) denotes the acceptance probability (7) of the Metropolis-Hastings kernel. We tested this variant but could not see any advantage over the standard sampler presented in the preceding sections. For the augment-resample type algorithm the implementation is more involved and the computational burden significantly higher. In particular, the Rao-Blackwellization effect one might achieve when replacing the accept-reject steps of the transition kernel by a single resampling step does not seem to justify the extra computational effort.
Indeed, aggregating the weights does not only prevent us from using the effective sample size criterion, but also requires extra computational time of O(n log n) in each iteration of the move step since pooling the weights is as complex as sorting. With our application in mind, however, computational time is more critical than memory, and we therefore recommend to refrain from aggregating the weights.

Parametric families on binary spaces
We review three parametric families on B d . In contrast to the similar discussion in Schäfer and Chopin (2011), we also consider a parametric family which cannot be used in sequential Monte Carlo samplers but in the context of the cross-entropy method. For more details on parametric families on binary spaces we refer to Schäfer (2012).

Suitable parametric families
We frame some properties making a parametric family suitable as proposal distribution in sequential Monte Carlo algorithms.
(a) For reasons of parsimony, we prefer a family of distributions with at most d(d + 1)/2 parameters like the multivariate normal.
(b) Given a sample X = (x 1 , . . . , x n ) from the target distribution π, we need to estimate θ * in a reasonable amount of computational time.
(c) We need to generate samples Y = (y 1 , . . . , y m ) from the family q θ . We need the rows of Y to be independent.
(d) For the sequential Monte Carlo algorithm, we need to evaluate q θ (y) point-wise. However, the cross-entropy method still works without this requirement.
(e) We want the calibrated family q θ * to reproduce e.g. the marginals and covariance structure of π to ensure that the parametric family q θ * is sufficiently close to π.

Product family
The simplest non-trivial distributions on B d are certainly those having independent components.
Product family For a vector m ∈ (0, 1) d of marginal probabilities, we define the product family

Properties
We check the requirement list from Section 3.1: (a) The product family is parsimonious with dim(θ) = d. (b) The maximum likelihood estimatorm is the weighted sample mean. (c) We can easily sample y ∼ q m . (d) We can easily evaluate the mass function q m (y). (e) However, the product family does not reproduce any dependencies we might observe in (w, X). The last point is the crucial weakness which makes the product family impractical for particle optimization algorithms on strongly multi-modal problems. Consequently, the rest of this section deals with ideas on how to sample binary vectors with a given dependence structure. There are, to our knowledge, two major strategies to this end.
(1) We construct a generalized linear model which permits to compute the conditional distributions. We apply the chain rule and write q θ as which allows to sample the entries of a random vector component-wise.
(2) We sample from an auxiliary distribution ϕ θ and map the samples into B d . We call a copula family, although we refrain from working with explicit uniform marginals.
We first present a generalized linear model and then review a copula approach.

Logistic conditionals family
Even for rather simple non-linear models we usually cannot derive closed-form expressions for the marginal probabilities required for sampling according to (11). Therefore, we might directly construct a parametric family from its conditional probabilities.
Logistic conditionals family We define, for a lower triangular matrix A ∈ R d×d , the logistic conditionals family as We readily identify the product family q m as the special case A = diag −1 (m).
The virtue of the logistic conditionals family is that, by construction, we can sample a random vector component-wise while the full probability q A (y) of the sample y is computed as a by-product of Procedure 3. We refer to the Online Supplement for instructions on how to fit the parameter A.
Procedure 3: Sampling via chain rule factorization

Properties
We check the requirement list from Section 3.1: (a) The logistic conditionals family is sufficiently parsimonious with dim(θ) = d(d + 1)/2. (b) We can fit the parameter A via likelihood maximization. The fitting is computationally intensive but feasible. (c) We can sample y ∼ q A using the chain rule factorization (11). (d) We can exactly evaluate q A (y). (e) The family q A reproduces the dependency structure of the data X although we cannot explicitly compute the marginal probabilities.

Gaussian copula family
Let ϕ θ be a family of multivariate auxiliary distributions on X and τ : X → B d a mapping into the binary space. We can sample from the copula family (12) by setting x = τ (v) for a draw v ∼ ϕ θ from the auxiliary distribution. Most multivariate parametric families with at most d(d + 1)/2 parameters appear to either have a rather limited dependency range or they do not scale to higher dimensions Joe (1996). Therefore, the natural and seemingly only viable option for ϕ θ is the multivariate normal distribution Emrich and Piedmonte (1991).
Gaussian copula family For a vector a ∈ R d and a correlation matrix Σ ∈ R d×d , we introduce the mapping and define the Gaussian copula family as , the cross-moments are equal the cumulative distribution function of the multivariate normal with respect to the entries indexed by I (see Schäfer (2012) for a more detailed discussion). In particular, the first and second moments are where Φ 1 (·) and Φ 2 (·, ·; σ ij ) denote the cumulative distribution functions of the univariate and bivariate normal distributions with zero mean, unit variance and correlation coefficient σ ij ∈ [−1, 1]. We refer to the Online Supplement for instructions on how to fit the parameters a and Σ.

Properties
We check the requirement list from Section 3.1: (a) The Gaussian copula family is sufficiently parsimonious with dim(θ) = d(d + 1)/2. (b) We can fit the parameters a and Σ via method of moments. However, the parameter Σ is not always positive definite. (c) We can sample y ∼ q gc a,Σ using y = τ a (v) with v ∼ ϕ Σ . (d) We cannot easily evaluate q gc a,Σ (y) since this requires computing high-dimensional integral expressions which is a computationally challenging problem in itself (see e.g. Genz and Bretz (2009)). The Gaussian copula family is therefore less useful for sequential Monte Carlo samplers but can be incorporated into the cross-entropy method reviewed in Section 4.2. (e) The family q gc a,Σ reproduces the exact mean and, possibly scaled, correlation structure.

Toy example
We briefly discuss a toy example to illustrate the usefulness of the parametric families.
For the quadratic function the associated probability mass function π(γ) ∝ exp(γ Fγ) has a correlation matrix which indicates that this distribution has considerable dependencies and its mass function is therefore strongly multi-modal. We generate pseudo-random data from π, adjust the parametric families to the data and plot the mass functions of the fitted parametric families. Figure 2 shows how the three parametric families cope with reproducing the true mass function. Clearly, the product family is not close enough to the true mass function to yield a suitable instrumental distribution while the logistic conditional family almost copies the characteristics of π and the Gaussian copula family allows for an intermediate goodness of fit.  (13).

Optimization algorithms
In this section, we provide a synopsis of all steps involved in the sequential Monte Carlo algorithm and connect this framework to the cross-entropy method and simulated annealing. In Table 1, we state the necessary formulas for the tempered and the level set sequence introduced in Section 1.2.

Sequential Monte Carlo
For convenience, we summarize the complete sequential Monte Carlo sampler in Algorithm 1. Note that, in practice, the sequence π t is not indexed by t but rather by t , which means that the counter t is only given implicitly.
The algorithm terminates if the particle diversity sharply drops below some threshold δ > 0 which indicates that the mass has concentrated in a single mode. If we use a kernel with proposals from a parametric family q θt , we might already stop if the family degenerates in the sense that only a few components of q θt , say less than d * = 12, are random while the others are constant ones or zeros. In this situation, additional moves using a parametric family are a pointless effort. We either return the maximizer within the particle system or we solve the subproblem of dimension d * by brute force enumeration. We might also perform some final local moves in order to further explore the regions of the state space the particles concentrated on.

Cross-entropy method
For the level set sequence, the effective sample size is the fraction of the particles which have an objective function value greater than max x∈B d f (x) − 1/( t + α); see Table 1 and equation (3). The remaining particles are discarded since their weights equal zero. Consequently, there is no need to explicitly compute α t as a solution of (6). We simply order the particles x k according to their objective values f (x k ) and only keep the n(1−β) particles with the highest objective values. Rubinstein Rubinstein (1997), who popularizes the use of level set sequences in the context of the cross-entropy method, refers to n(1 − β) as the size of the elite sample. The cross-entropy method has been applied successfully to a variety of combinatorial optimization problems, some of which are equivalent to pseudo-Boolean optimization Rubinstein and Kroese (2004), and is closely related to the proposed sequential Monte Carlo framework.
However, the central difference between the cross-entropy method and the sequential Monte Carlo algorithm outlined above is the use of the invariant transition kernel in the latter. We obtain the cross-entropy method as a special case if we replace the kernel κ t by its proposal distribution q θt . The sequential Monte Carlo approach uses a smooth family of distributions {π : ≥ 0} and explicitly schedules the evolution π t which in turn leads to the proposal distributions q θt . The cross-entropy method, in contrast, defines the subsequent proposal distribution q θ t+1 ≈ q θt 1 L + t+1 without any reference sequence π t to balance the speed of the particle evolution.
In order to decelerate the advancement of the cross-entropy method, we introduce a lag parameter τ ∈ [0, 1) and use a convex combination of the previous parameter θ t−1 and the parameterθ t fit to the current particle system, setting However, there are no guidelines on how to adjust the lag parameter during the run of the algorithm. Therefore, the sequential Monte Carlo algorithm is easier to calibrate since the reference sequence π t controls the stride and automatically prevents the system from overshooting.
On the upside, the cross-entropy method allows for a broader class of auxiliary distributions q θt since we do not need to evaluate q θt (x) point-wise which is necessary in the computation of the acceptance probability of the Hastings kernel; see Section 3.4.

Simulated annealing
A well-studied approach to pseudo-Boolean optimization is simulated annealing Kirkpatrick et al. (1983). While the name stems from the analogy to the annealing process in metallurgy, there is a pure statistical meaning to this setup. We can picture simulated annealing as approximating the mode of a tempered sequence (2) using a single particle. Since a single observation does not allow for fitting a parametric family, we have to rely on symmetric transition kernels (8) in the move step.
A crucial choice is the sequence t which in this context is often referred to as the cooling schedule. There is a vast literature advising on how to calibrate t where a typical guideline is the expected acceptance rate of the Hastings kernel. We calibrate t such that the empirical acceptance rate λ t−s:t := t r=t−s λ r , s > 0 follows approximately (t + 1) −5 for t ∈ [0, 1]. There are variants of simulated annealing which use more complex cooling schedules, tabu lists and multiple restarts, but we stick to this simple version for the sake of simplicity. Algorithm 2 describes the version we use in our numerical experiments in Section 5.3.3.

Randomized local search
We describe a greedy local search algorithm which works on any state space that allows for defining a neighborhood structure. The typical neighborhood on binary spaces is the k-neighborhood defined in (9). A greedy local search algorithm computes the objective value of all states in the current neighborhood and moves to the best state found until a local optimum is reached. The local search algorithm is called k-opt if it searches the neighborhood ∪ k i=1 N i (·) (see e.g. Merz and Freisleben (2002) for a discussion). The algorithm can be randomized by repeatedly restarting the procedure from randomly drawn starting points. There are more sophisticated versions of local search algorithms exploit the properties of the objective function but even a simple local search procedure can produce good results Alidaee et al. (2010). Algorithm 3 describes the 1-opt local search procedure we use in our numerical experiments in Section 5.3.3.
Algorithm 3: Randomized local search

Introduction
It is well-known that any pseudo-Boolean function f : B d → R can be written as a multi-linear function where a I ∈ R are real-valued coefficients. We say the function f is of order k if the coefficients a I are zero for all I ⊆ [[1, d]] with |I| > k. While optimizing a first order function is trivial, optimizing a non-convex second order function is already an NP-hard problem Garey and Johnson (1979).
In the sequel, we focus on optimization of second order pseudo-Boolean functions to exemplify the stochastic optimization schemes discussed in the preceding sections. If f is a second order function, we restate program (1) as where F ∈ R d×d is a symmetric matrix. We call (15) an unconstrained quadratic binary optimization problem (uqbo); we refer to Boros et al. (2007) for a list of applications and equivalent problems. In the literature it is also referred to as unconstrained quadratic Boolean or bivalent or zero-one programming Beasley (1998).

Particle optimization and meta-heuristics
Meta-heuristics are a class of algorithms that optimize a problem by improving a set of candidate solutions without systematically enumerating the state space; typically they deliver solutions in polynomial time while an exact solution has exponential worst case running time. The outcome is neither guaranteed to be optimal nor deterministic since most meta-heuristics are randomized algorithms. We briefly discuss the connection to particle optimization against the backdrop of the unconstrained quadratic binary optimization problem where we roughly separate them into two classes: local search algorithms and particle-driven meta-heuristics. Local search algorithms iteratively improve the current candidate solution through local search heuristics and judicious exploration of the current neighborhood; examples are local search Boros et al. (2007); Merz and Freisleben (2002), tabu search Glover et al. (1998); Palubeckis (2004), simulated annealing Katayama and Narihisa (2001). Particle driven meta-heuristics propagate a set of candidate solutions and improve it through recombination and local moves of the particles; examples are genetic algorithms Merz and Freisleben (1999), memetic algorithms Merz and Katayama (2004), scatter search Amini et al. (1999). For comparisons of these methods we refer to Hasan et al. (2000) or Beasley (1998).
The sequential Monte Carlo algorithm and the cross-entropy method are clearly in the latter class of particle-driven meta-heuristics. The idea behind sequential Monte Carlo is closely related to the intuition behind population (or swarm) optimization and genetic (or evolutionary) algorithms. However, the mathematical framework used in sequential Monte Carlo allows for a general formulation of the statistical properties of the particle evolution while genetic algorithms are often problem-specific and empirically motivated.

Particle optimization and exact solvers
If we can explicitly derive the multi-linear representation (14) of the objective function, there are techniques to turn program (1) into a linear program. For the uqbo it reads Note, however, that there are more parsimonious linearization strategies than this straightforward approach [ (Hansen and Meyer, 2009), (Gueye and Michelon, 2009)]. The transformed problem allows to access the tool box of linear integer programming which consist of branch-and-bound algorithms that are combined with rounding heuristics, various relaxations techniques and cutting plane methods [ (Pardalos and Rodgers, 1990), (Palubeckis, 1995)]. Naturally, the question arises whether particle-driven meta-heuristics can be incorporated into exact solvers to improve branch-and-bound algorithms. Indeed, stochastic meta-heuristics deliver lower bounds for maximization problems, but particle-driven algorithms are computationally somewhat expensive for this purpose unless the objective function is strongly multi-modal and other heuristics fail to provide good results; see the discussion in Section 5.2.5.
However, the sequential Monte Carlo approach in combination with the level set sequence (3) might also be useful to determine a global branching strategy, since the algorithm provides an estimator for which is the average of the super-level set L + c := {x ∈ B d : f (x) ≥ c}. These estimates given for a sequence of levels c might provide branching strategies than are superior to local heuristics or branching rules based on fractional solutions. A further discussion of this topic is beyond the scope of this paper but it certainly merits consideration.

Introduction
The meta-heuristics we want to compare do not exploit the quadratic structure of the objective function and might therefore be applied to any binary optimization program. If the objective function can be written in multi-linear form like (15) there are efficient local search algorithms Boros et al. (2007); Merz and Freisleben (2002) which exploit special properties of the target function and easily beat particle methods in terms of computational time.
Therefore, the use of particle methods is particularly interesting if the objective function is expensive to compute or even a black box. The posterior distribution in Bayesian variable selection for linear normal models is an example of such an objective function (see Schäfer and Chopin (2011) and references therein). We stick to the uqbo for our numerical comparison since problem instances of varying difficulty are easy to generate and interpret while the results carry over to general binary optimization.
In the vast literature on uqbo, authors typically compare the performance of metaheuristics on a suite of randomly generated problems with certain properties. Pardalos (1991) proposes standardized performance tests on symmetric matrices F ∈ Z d×d with entries f ij drawn from the uniform The test suites generated by Beasley (1990, OR-library) and Glover et al. (1998) follow this approach have been widely used as benchmark problems in the uqbo literature (see Boros et al. (2007) for an overview). In the sequel we discuss the impact of diagonal dominance, shifts, the density and extreme values of F on the expected difficulty of the corresponding uqbo problem.

Diagonal
Generally, stronger diagonal dominance in F corresponds to easier uqbo problems Billionnet and Sutter (1994). Consequently, the original problem generator presented by Pardalos (1991) is designed to draw the off-diagonal elements from a uniform on a different support [[−q, q]] with q ∈ N.
In this context, we point out that the impact of diagonal dominance carries over to the statistical properties of the tempered distributions (2) we defined in the introductory Section 1.2. Indeed, stronger diagonal dominance in F corresponds to exponential quadratic distributions π(γ) := exp(γ Fγ) γ∈B d exp(γ Fγ) having lower dependencies between the components of γ. We can analytically derive a parameter A ∈ R d×d for a logistic conditionals family q A that approximates π(γ) where the quality of the approximation increases as the diagonal of F becomes more dominant Schäfer (2012). We can accelerate the sequential Monte Carlo algorithm by initializing the system from q A instead of U B d . However, we did not exploit this option to keep the present work more concise.
For positive definite F 0, the optimization problem is convex and can be solved in polynomial time Kozlov et al. (1979); in exact optimization, this fact is exploited to construct upper bounds for maximization problems Poljak and Wolkowicz (1995). We observe a corresponding complexity reduction in statistical modeling. For F 0, the auxiliary distribution is a feasible mass function, and we can derive analytical expressions concerning all crossmoments and marginal distributions Schäfer (2010) which allows to largely analyze the properties of π(γ) without enumerating the state space.

Shifts
The global optimum of the uqbo problem is more difficult to detect as we shift the entries of the matrix F but the relative gap between the optimum and any heuristic value diminishes. If we sample f ij = f τ ij from a uniform on the shifted support we obtain an objective function where d = means equality in distribution. Hence, with growing |τ | the optimum depends less on F and the relative gap between the optimum and a solution provided by any meta-heuristic vanishes. Boros et al. (2007) define a related criterion and report a significant impact ofρ on the solution quality of their local search algorithms which is not surprising.

Density
The difficulty of the optimization problem is related to the number of interactions, that is the number of non-zero elements of F. We call the proportion of non-zeros the density of F. Drawing f ij from the mixture we adjust the difficulty of the problem to a given expected density ω.
Note that not all algorithms are equally sensitive to the density of F. Using the basic linearization (16), each non-zero off-diagonal element requires the introduction of an auxiliary variable and three constraints. Thus, the expected total number of variables and the expected total number of constraints, which largely determine the complexity of the optimization problem, are proportional to the density ω.
On the other hand, many randomized approaches, including the particle algorithms discussed in Section 2, are less sensitive to the density of the problem in the sense that replacing zero elements by small values has a minor impact on the performance of these algorithms. Rather than the zero/non-zero duality, we suggest that the presence of extreme values determines the difficulty of providing heuristic solutions.

Extreme values
The uniform sampling approach advocated by Pardalos (1991) is widely used in the literature for comparing meta-heuristics. Certainly, particle-driven methods are computationally too expensive to outperform local search heuristics on test problems with uniformly drawn entries; Beasley (1998) confirms this intuition with respect to genetic algorithms versus tabu search and simulated annealing. However, the uniform distribution does not produce extreme values and it is vital to keep in mind that these have an enormous impact on the performance of local search algorithms.
Extreme values in F lead to the existence of distinct local maxima x * ∈ B d of f in the sense that there is no better candidate solution than x * in the neighborhood ∪ k i=1 N i (x * ) even for relatively large k. Further, extreme local minima might completely prevent a local search heuristic from traversing the state space in certain directions. Consequently, local search algorithms, as discussed in Section 5.1.2, depend more heavily on their starting value, and their performance deteriorates with respect to particle-driven algorithms.
We propose to draw the matrix entries f ij from a discretized Cauchy distribution that has heavy tails which cause extreme values to be frequently sampled. Figure 3 shows the distribution of a Cauchy and a uniform to illustrate the difference. The resulting uqbo problems have quite distinct local maxima; in that case we also say that the function f (x) is strongly multi-modal.

Numerical comparison
In this section, we provide numerical comparisons based on instances of the uqbo problem. We generated two random test suites of dimension d = 250, each having 10 instances. For the first suite, we sampled the matrix entries from a uniform distribution U 100 on [[−100, 100]]; for the second, we sampled from a Cauchy distribution C 100 as defined in (17). For performance evaluation, we run a specified algorithm 100 times on the same problem and denote the outcome by x 1 , . . . , x 100 .

Visualization
Since the absolute values are not meaningful, we report the relative ratios where the best known solution is the highest objective value ever found for that instance and the worst solution is the lowest objective value among the 100 outcomes. We summarize the results in a histogram. The first n bins are singletons b k := { * k } for the highest values * 1 > · · · > * n ∈ { k : k ∈ [[1, 100]]}; the following n bins are equidistant intervals b < k := [ n−k n * n , n−k+1 n * n ). The graphs show the bins b 1 , . . . , b n , b < 1 , . . . , b < n in descending order from left to right on the x-axis. The interval bins are marked with a sign "<" and the lower bound. The y-axis represents the counts.
For comparison, we draw the outcome of several algorithms into the same histogram, where the worst solution found is the lowest overall objective value among the outcomes. For each algorithm, the counts are depicted in a different color and, for better readability, with diagonal stripes in a different angle. To put it plainly, an algorithm performs well if its boxes are on the left of the graph since this implies that the outcomes where often close to the best known solution.

Comparison of binary parametric families
We study how the choice of the binary parametric family affects the quality of the delivered solutions. The focus is on the cross-entropy method, since we cannot easily use the Gaussian copula family in the context of sequential Monte Carlo. We use n = 1.2 × 10 4 particles, set the speed parameter to β = 0.8 (or the elite fraction to 0.2) and the lag parameter to τ = 0.5.
The numerical comparisons, given in Figures 5(b) and 5(a), clearly suggest that using more advanced binary parametric families allows the cross-entropy method to detect local maxima that are superior to those detected using the product family. Hence, the numerical experiments confirm the intuition of our toy example in Figure 2.
On the strongly multi-modal instance 5(a) the numerical evidence for this conjecture is stunningly clear-cut; on the weakly multi-modal problem 5(b) its validity is still unquestionable. This result seems natural since reproducing the dependencies induced by the objective function is more relevant in the former case than in the latter.

Comparison of optimization algorithms
We compare a sequential Monte Carlo sampler with parametric family, a sequential Monte Carlo sampler with a single-flip symmetric kernel (8), the cross-entropy method, simulated annealing and 1-opt local search as described in Section 4.
For the cross entropy method, we use the same parameters as in the preceding section. For the sequential Monte Carlo algorithm, we use n = 0.8 × 10 4 particles and set the speed parameter to β = 0.9; we target a tempered auxiliary sequence (2). For both algorithms we use the logistic conditionals family as sampling distribution. With these configurations, the algorithms converge in roughly 25 minutes. We calibrate the sequential Monte Carlo sampler with local moves to have the same average run time by processing batches of 10 local moves before checking the particle diversity criterion. The simulated annealing and 1-opt local search algorithms run for exactly 25 minutes.
The results shown in Figures 6(b) and 6(a) assert the intuition that particle methods perform significantly better on strongly multi-modal problems. However, on the easy test problems, the particle methods tend to persistently converge to the same sub-optimal local modes. This effect is probably due to their poor local exploration properties. Since particle methods perform significantly less evaluations of the objective function, they are less likely to discover the highest peak in a region of rather flat local modes. The use of parametric families aggravates this effect, and it seems advisable to alternate global and local moves to make a particle algorithm more robust against this kind of behavior.
Further numerical results are shown in Figure 6 and Figure 7.

Discussion and conclusion
The numerical experiments carried out on different parametric families revealed that the use of the advanced families proposed in this paper significantly improves the performance of the particle algorithms, especially on the strongly multi-modal problems. The experiments demonstrate that local search algorithms, like simulated annealing and randomized 1-opt local search, indeed outperform particle methods on weakly multi-modal problems but deliver inferior results on strongly multi-modal problems. Using tabu lists, adaptive restarts and rounding heuristics, we can certainly design local search algorithms that perform better than simulated annealing and 1-opt local search. Still, the structural problem of strong multi-modality persists for path-based algorithms. On the other hand, cleverly designed local search heuristics will clearly beat sequential Monte Carlo methods on easy to moderately difficult problems.
The results encourage the use of particle methods if the objective function is known to be potentially multi-modal and hard to analyze analytically. We have to keep in mind that multiple restarts of rather simple local search heuristics can be very efficient if they make use of the structure of the objective function. For 25 minutes of randomized restarts, the heuristic proposed by Boros et al. (2007), which exploits the fact that the partial derivatives of a multi-linear function are constant, practically always returns the best known solution on all test problems treated to create Figures 6 and 7.
The numerical work was completely done in Python 2.6 using SciPy packages and run on a cluster with 1.86 GHz processors. The sources used in this work and the problems processed in this paper can be found at http://code.google.com/p/smcdss.

Acknowledgments
This work is part of the author's Ph.D. thesis at CREST under supervision of Nicolas Chopin whom I would like to thank for the numerous discussions on particle algorithms. I thank the editor and two anonymous referees for their detailed comments which helped to significantly improve this paper.

Appendix: Fitting the parameters
We briefly summarize how the parameters of the logistic conditionals family and the Gaussian copula family can be assessed for a given particle system X = (x 1 , . . . , x n ) with weights w = (w 1 , . . . , w n ). We denote bȳ the weighted first and second sample moments.
where (p

Complete separation
The data might suffer from complete or quasi-complete separation Albert and Anderson (1984) which causes the likelihood function L(a) to be monotonic. In that case there is no maximizer. We can avoid the monotonicity by assigning a suitable prior distribution to the parameter a. Firth (1993) recommends the Jeffreys prior for its bias reduction which can conveniently be implemented via data adjustment Kosmidis and Firth (2009).
For the sake of simplicity, however, we only assign a simple Gaussian prior with variance ε −1 > 0 such that, up to a constant, the log-posterior distribution is the loglikelihood function plus a quadratic penalty term and therefore always convex. The score function and its Jacobian matrix become The bias-reduced estimators are known to shrink towards 0.5 which is an undesired property when fitting a parametric family. Therefore, we attempt to keep the shrinkage parameter ε as small as possible.
If the Newton iteration at the ith component fails to converge, we can either augment the penalty term ε which leads to stronger shrinkage of the mean m i towards 0.5 or we can drop some covariates γ j for j ∈ {1, . . . , i − 1} from the iteration to improve the numerical condition of the procedure.
In practice, we drop the predictors from the regression model which are only weakly correlated with the explained variable. This step is important to speed up the algorithm and improve its numerical properties. For a proposal distribution, it is particularly important to take the strong dependencies into account but it is often sufficient to work with very sparse logistic conditionals families.
In particularly difficult cases, we might prefer to set a ii = −1 (x i ) and a i,1:i−1 = 0, wherex i is defined in (18), which guarantees that at least the mean is correct. This is an important issue since misspecification of the mean of γ i also affects the distribution of the components γ j for j ∈ {i + 1, . . . , d} which are sampled conditional on γ i .

Gaussian copula family
We adjust the Gaussian copula family q gc a,Σ by method of moments. We need to solve the non-linear equations wherex i andx ij are defined in (18) while Φ 1 (v i ) and Φ 2 (v i , v j ; σ ij ) denote the cumulative distribution functions of the univariate and bivariate normal distributions with zero mean, unit variance and correlation coefficient σ ij ∈ [−1, 1]. While the parameter a i = Φ −1 1 (x i ) is easy to assess, the challenging task is to compute the bivariate variances σ ij for all i, j ∈ [ [1, d]]. Recall the standard result [ (Johnson et al., 2002), p.255] ∂Φ 2 (x 1 , x 2 ; σ) ∂σ = ϕ(x 1 , x 2 ; σ), where ϕ(·, ·; σ) denotes the density of the bivariate normal distribution with correlation coefficient σ. We obtain the following Newton-Raphson iteration starting at some initial value σ ij ∈ (−1, 1); see Procedure 5. In the sequential Monte Carlo context, good initial values are obtained from the parameters of the previous auxiliary distributions.
We use a fast series approximation Drezner and Wesolowsky (1990) to evaluate Φ 2 (a i , a j ; σ ij ). These approximations are critical when σ ij comes very close to either boundary of [−1, 1]. The Newton iteration might repeatedly fail when restarted at the corresponding boundary σ (0) ij ∈ {−1, 1}. In any event, Φ 2 (x 1 , x 2 ; σ) is strictly monotonic in σ since its derivative (20) is positive, and we can switch to bi-sectional search if necessary.
if not Σ 0 then Σ ← (Σ + |λ| I)/(1 + |λ|) return a, Σ The locally fitted correlation matrices Σ might not be positive definite for d ≥ 3 because the Gaussian copula is not flexible enough to model the full range of crossmoments binary distributions might have (see Schäfer (2012) for an extended discussion). We obtain a feasible parameter replacing Σ by Σ * = (Σ + |λ| I)/(1 + |λ|), where λ is smaller than all eigenvalues of the locally fitted matrix Σ. This approach evenly lowers the local correlations to a feasible level and is easy to implement on standard software. In practice, we also prefer to work with a sparse version of the Gaussian copula family that concentrates on strong dependencies and sets minor correlations to zero.