Cross-entropy and rare events for maximal cut and partition problems.
We show how to solve the maximal cut and partition problems using a randomized algorithm based on the cross-entropy method. For the maximal cut problem, the proposed algorithm employs an auxiliary Bernoulli distribution, which transforms the original deterministic network into an associated stochastic one, called the associated stochastic network (ASN). Each iteration of the randomized algorithm for the ASN involves the following two phases:(1) Generation of random cuts using a multidimensional Ber(p) distribution and calculation of the associated cut lengths (objective functions) and some related quantities, such as rare-event probabilities.(2) Updating the parameter vector p on the basis of the data collected in the first phase.We show that the Ber(p) distribution converges in distribution to a degenerated one, Ber(pd&ast;), in the sense that someelements of pd&ast;, will be unities and the rest zeros. The unity elements of pd&ast; uniquely define a cut which will be taken as the estimate of the maximal cut. A similar approach is used for the partition problem. Supporting numerical results are given as well. Our numerical studies suggest that for the maximal cut and partition problems the proposed algorithm typically has polynomial complexity in the size of the network.
Introduction
Most combinatorial optimization problems are NP-hard; for example, deterministic and
stochastic (noisy) scheduling, the traveling salesman problem (TSP), the maximal cut in
a network, the longest path in a network, optimal buer allocation in a production line,
optimal routing in deterministic and stochastic networks and
ow control, optimization
of topologies and conguration of computer communication and tra-c systems. Well
established stochastic methods for combinatorial optimization problems are simulated annealing
[1], [2], [7], [37], initiated by Metropolis [31] and later generalized in [19], [22] and
[26], tabu search [16] and genetic algorithms [17]. For some additional references on both
deterministic and stochastic combinatorial optimization see [3]-[4], [23]-[25], [30], [32] and
[33]-[36].
Recent works on stochastic combinatorial optimization, which is also a subject of this
paper, include the method of Andradottir [5], [6], the nested partitioning method (NP)
[42], [43], the stochastic comparison method [18], and the ant colony optimization (ACO)
meta heuristic of Dorigo and colleagues [12], [15]. In most of the above methods a Markov
chain is constructed and almost sure convergence is proved by analyzing the stationary
distribution of the Markov chain.
We shall next review brie
y the ACO meta heuristic algorithms of [8]-[15], [20], [41],
[44]-[47], which try to mimic ant colonies behavior. It is known that ant colonies are able
to solve shortest-path problems in their natural environment by relying on a rather simple
biological mechanism: while walking, ants deposit on the ground a chemical substance,
called pheromone. Ants have a tendency to follow these pheromone trails. Within a xed
period, shorter paths between nest and food can be traversed more often than longer
paths, and so they obtain a higher amount of pheromone, which, in turn, tempts a larger
number of ants to choose them and thereby to reinforce them again. The above behavior
of real ants has inspired many researcher to use the ant system models and algorithms in
which a set of articial ants cooperate via pheromone depositing either on the edges or on
the vertices of the graph. Consider, for example, the ACS (ant colony system) approach of
Dorigo, Maniezzo and Colorni [15] for solving the TSP problem, which can be described as
follows: rst, a number of \articial ants", also called agents, are positioned randomly at
some node of the graph. Then, each agent performs a series of random moves to neighbor
nodes, controlled by suitably dened transition probabilities. Once an agent has visited all
nodes, the length of the tour is evaluated, and the pheromone values assigned to the arcs of
the path are increased by an amount depending on the length of the tour. This procedure
is repeated many times. The probability of a transition along a specic arc is computed
based on the pheromone value assigned to this arc, and the length of the arc. The higher
the pheromone value and the shorter the length of the arc, the higher the probability that
the agent will follow this arc in his rst move. Note also that while updating the transition
probabilities at each iteration of the ACS algorithm, Dorigo, Maniezzo and Colorni
[15] also introduce the so-called evaporation mechanism which discounts the pheromone
values obtained at the previous iteration. Diverse modications of ACS algorithm, which
presents a natural generalization of stochastic greedy heuristics, have been applied eciently
to many dierent types of discrete optimization problems and have produced very
results. Recently, the approach has been extended by Dorigo and Di Caro [12]
to a full discrete optimization meta heuristic, called the Ant Colony Optimization (ACO)
meta heuristic, which covers most of the well known combinatorial optimization problems.
Gutjahr [20], [21] was the rst to prove the convergence of the ACS algorithm.
This paper deals with application of the cross-entropy (CE) method to the maximal
cut and the partition problems. The CE method was rst introduced in [28] for estimating
probabilities of rare events for complex stochastic networks and then applied in [38]
for solving continuous multi-extremal and combinatorial optimization problems (COP),
namely, the shortest path (between two given nodes in a deterministic network or graph,
where each edge has a given length), the longest path and the TSP. The CE method for
combinatorial optimization employs an auxiliary random mechanism equipped with a set
of parameters, which transforms the deterministic network into a stochastic one, called
the associated stochastic network (ASN). Each iteration of the CE algorithm based on the
ASN involves the following two phases:
1. Generation of random trajectories (walks) using an auxiliary random mechanism,
like an auxiliary Markov chain with transition probability matrix
then calculation of the associated objective function.
2. Updating the parameters, like updating the elements of the transition probability
matrix on the basis of the data collected in the rst phase.
Let the original deterministic network be denoted by the graph
V is the set of nodes and E is the set of edges. Depending on a particular problem, we
introduce the randomness in the ASN by associating some form of randomness either to
(a) the edges E or to (b) the nodes V. More specically, we distinguish between the
so-called (a) stochastic edge networks (SEN) and (b) stochastic node networks (SNN).
(a) Stochastic edge networks (SEN). Here the trajectories are typically generated
using a Markov chain with transition probability matrix such that the transition
from i to j in uniquely denes the edge (ij) in the network. To SEN one can readily
reduce the TSP, the quadratic assignment problem, the deterministic and stochastic
ow
shop models, as well as some others. SEN were considered and treated earlier in [38].
(b) Stochastic node networks (SNN). Here the trajectories (walks) are generated
using an n-dimensional discrete distribution, like the n-dimensional Bernoulli (Ber (p))
distribution, such that each component of the random vector (rv)
uniquely denes the node V k network. To SNN one can
readily reduce the maximal cut problem, the partition problem, the clique problem, the
optimal buer allocation in a production line as well as some others. Application of CE
to SNN and in particular to maximal cut and maximal partition problems is the subject
of this paper.
Notice that terminology similar to SNN and SEN exists in Wagner, Lindenbaum and
Bruckstein [47] for the graph covering problem, called vertex ant walk (VAW) and edge ant
walk (EAN), respectively.
It is crucial to understand that the CE Algorithm 3.2 in [38] for SEN and Algorithm 4.1
in this paper for SNN are very similar. Their main dierence is in the sample trajectories
generation. As mentioned, in the former [38], the trajectories are typically generated using
a Markov chain, while in the latter they are generated, say, using an n-dimensional Ber(p)
distribution.
We shall show that the SNN Algorithm 4.1 for the maximal cut problem has the
following properties (similar properties apply for the maximal partitition problem):
1. The multi-dimensional Bernoulli distribution converges to a degenerated one
Ber(p
d
in the sense that some parameters of p
d , will be
unities and the rest will be zeros.
2. The unity elements of p
d uniquely dene a cut which will be taken as an estimate
of the maximal cut.
This perfectly matches with the SEN Algorithm 3.2 in [38], where
1. The probability matrix converges to a degenerated one in the sense that only
a single element at each row equals unity, while the remaining elements are
equal to zero.
2. The unity elements of the degenerated matrix uniquely dene the shortest tour, say
in a TSP problem.
We shall nally show that Algorithm 4.1, in fact, presents a simple modication of
Algorithm 2.1 for the estimation of probabilities of rare events. Algorithm 2.1 is the same
as Algorithm 1.2 in [38], and was adapted for the estimation of rare event probabilities in
SEN problems, like the stochastic TSP (Algorithm 3.1 in [38]); Algorithm 1.2 in [38] also
formed the basis of Algorithm 3.2 of that paper.
To elaborate more on this, let
represent the probability of the rare-event I fM(X)>xg ), where M(X) is the sample performance
of a stochastic system, vector with known
distribution and x is a xed number, chosen such that the probability '(x) is very small.
In order to give an example of M(X) in (1.1), consider a graph, whose edges have
random lengths given by X i 's. Then M(X) may be the length of the shortest path between
two designated nodes of the graph called the source and the sink. More formally, M(X)
can be dened as
is the j-th complete path from a source to a sink; p is the number of complete
paths.
Similar to Algorithm 3.1 in [38], Algorithm 2.1 (in this paper) is an adaptive algorithm
for the estimation of rare event probabilities using importance sampling and cross-entropy.
A distinguishing feature of both algorithms is that, when x is not xed in advance, they
automatically generate a sequence of tuples f
(see (2.5) and (2.8) below) ensuring
that
and '(
is the iteration number of Algorithm 2.1. The algorithm
stops when
Turning to COPs, note that as soon as a deterministic COP is transformed to a
stochastic one and M(X) (called the sample performance of the ASN, e.g., the cut value)
is available we can cast our ASN into the rare event framework (1.1). (Recall that in the
original formula (1.1) X is a natural random vector, while in the ASN it is an articially
constructed random vector, say a Bernoulli random vector). We shall show that in analogy
to Algorithm 2.1, Algorithm 4.1 generates a sequence of tuples f
(see (4.6) and
(4.7) below), which converges in distribution to a stationary point (
is the
true maximal cut value and p
d is the degenerated vector determining the cut corresponding
to
. In the language of rare events this also means that Algorithm 4.1 is able to identify
with very high probability a very small subset of the largest cut values. In what follows, we
shall show that COP's can be solved simultaneously with estimation of the probabilities of
rare-events for the ASN. This framework enables us to establish tight connections between
rare-events and combinatorial optimization.
As was also the case in [38], it is not our goal here to compare the e-ciency of the
proposed method with other well-established alternatives, such as simulated annealing,
tabu search and genetic algorithms. This will be done somewhere else. Our goal is
merely to establish some theoretical foundation for the proposed CE method, and to
demonstrate, both theoretically and numerically, the high speed of convergence of the
proposed algorithm and promote our approach for further applications.
In Section 2 we review the adaptive algorithm for the estimation of rare event probabilities
by citing some material from [28], [39] and Section 1.2 of [38]. In Section 3 we present
the maximal cut and partition problems, and the ASN for which we dene a probability
of rare-event exactly as in (1.1). We also present algorithms for generation of random cuts
and partitions. Section 4 presents our main CE Algorithm 4.1 for the maximal cut and
partition problems. Here we also prove a theorem stating that under certain conditions,
the sequence of tuples f
associated with the ASN converges to the stationary
point
d ). In Section 5 we give some modications to our main CE algorithm; an
important modication that we present is the fully automated CE Algorithm. In Section
6 supportive numerical results are presented and in Section 7 concluding remarks and
directions for further research are given.
2 The Cross-Entropy Method for Probability of Rare
Events Estimation
Let f(z; v) be a multivariate density with parameter vector v. Consider estimating
v). The importance sampling estimate is given by
new
I
where
is the likelihood ratio, X i  f(z; v 0 ) and v new is called the reference parameter.
To nd the optimal reference parameter v
new one can either minimize the variance of
the importance sampling estimate
new ) (see [29], [39]) or maximize the following
cross-entropy (see [38])
new )g
Below we do the latter. Given a sample new ), we can estimate
the optimal solution v
new of (2.2) by the optimal solution of the program
vnew
new
I
It is readily seen that the programs (2.2) and (2.3) are useful only in the case, where ' is
not very small, say In rare-event context (say,
and (2.3) are useless, since owing to the rarity of the events fM(X i )  xg, the random
variables I fM(X i )xg
and the associated derivatives of b
new ) at
probability, provided the sample N is small relative to the
reciprocal of the rare-event probability '(x). To overcome this di-culty we introduce an
auxiliary sequence f
We start by choosing an initial
under the original pdf f(z; v), the probability '(
g g is not too small,
say '(
specically, we set v
sequentially iterate
in both v
t and
t as outlined below.
(a) Adaptive estimation of
t . For a xed v  t derive
t from the following simple
one-dimensional root-nding program
I fM(X)
The stochastic counterpart of (2.4) is as follows: for xed
derive
t from the following
program
I
It is readily seen that
where M t;(j) is the j-th order statistics of the sequence M t;j  M(X j ).
(b) Adaptive estimation of v
t . For xed
t from the solution of the
program
I fM(X)
where v  v.
The stochastic counterpart of (2.7) is as follows: for xed
derive
t from the
following program
I
where
The resulting algorithm for estimating '(x) can be written as
1. Set v  0   v  0  v. Generate a sample X XN from the pdf f(x; v  0 ) and deliver
the solution (2.6) of the program (2.5). Denote the initial solution by
. Set t=1.
2. Use the same sample X as in (2.5) and solve the stochastic program (2.8).
Denote the solution by
3. Generate a new sample X XN from the pdf f(x;
deliver the solution
t in (2.6) of the program (2.5). Denote the solution by
t .
4. If
t  x, set
t  x and solve the stochastic program (2.8) for
x. Denote the
solution as
t+1 and stop; otherwise set t  reiterate from step 2. After
stopping:
Estimate the rare-event probability '(x) using the estimate (2.1), with v 0 replaced
by
t+1 .
The monotonicity of the sequence
crucial for convergence of Algorithm
2.1. It is proved in [27] that if X is a one dimensional random variable that is distributed
Gamma(v;) and M() is a monotonically increasing positive function of one variable,
then the sequence
generated by Algorithm 2.1 monotonically increases, provided
N !1 at each iteration t. This theorem can be readily extended for multidimensional X
and some other distributions from the exponential family, such as Normal, Beta, Poisson
and discrete. For more details see [27].
The following theorem, due to [28], states that if the sequence f
monotonically increasing, then under some mild regularity conditions, as N ! 1, the
sequence f
reaches x in a nite number of iterations.
Theorem 2.1 Let x be such that '(x) > 0. Let h be the mapping which corresponds to an
iteration of Algorithm 2.1, i.e.,
(v)
(v).
Assume rst that the following conditions hold
1. The sequence f
(v t )g; monotonically increasing.
2. The mapping v 7!
(v) is continuous.
3. The mapping v 7!
(v) is proper, i.e. if
(v) belongs to some closed interval then
v belongs to a compact set.
4. The mapping v 7!
(v) is lower semi-continuous.
Then there exist t < 1 such that
lim
Pf
Proof Given in [27].
It readily follows from the above that if x is not xed in advance, Algorithm 2.1 will
automatically generate two sequences:
such that '(
3 The Maximal Cut and Partition Problems
3.1 Cuts, Partition and the Associated Stochastic Networks
The maximal cut problem in a graph can be formulated as follows. Given a graph
E) with set of nodes set of edges E between the nodes,
partition the nodes of the graph into two arbitrary subsets V 1 and V 2 such that the sum of
the weights of the edges going from one subset to the other is maximized. Mathematically
it can be written as
f ~
where
~
and denotes the symmetric matrix of weights (distances) of the edges, which is
assumed to be known.
The partition problem can be dened similarly. The only dierence between the maximal
cut and the partition problem is that in the former the length, say  , of the vector
while in the latter it is xed.
Note that solving (3.1), for each one needs to decide whether
Since the matrix (L ij ) is symmetric, that is in order to avoid
duplication we shall assume without loss of generality that
The program (3.1) can be also written as
where
is the length (value) of the k-th cut, called the objective function, (V is the k-th
is the set of all possible cuts in the graph and jX j is the cardinality
of the set fXg. We denote the maximal cut and the maximal cut value (the optimal value
of the objective function)by (V
, respectively.
It is readily seen that the total number of cuts is
Similarly, the total number of partitions with  being xed and equal to n=2 (for simplicity
assume n to be even) is

Figure

3.1: A 6-node network
As an example, consider Figure 3.1 and the associated 6x6 distance matrix
with cardinalities for maximal cut and partition jX
Consider, for instance, the following two cuts
and
The function value (partition cost)
in the rst and the second cases are
and
respectively.
As mentioned before, in order to generate the stationary tuple (
d ) (see Algorithm
4.1 below), we need to transform the original (deterministic) network into an associated
stochastic one. To do so for the maximal cut problem, we associate an n dimensional
random vector (rv) with the n dimensional vector
Each component of X is independent and Bernoulli distributed, i.e., X k  Ber(p k ), and
has the interpretation that if X . If not stated
otherwise, we set p 1  1
Then for the maximal cut problem, each iteration of our main Algorithm 4.1 comprises
the following two phases:
(a) Generation of random cuts (see Algorithms 3.1, below) from the ASN using the
calculating the associated sample performance M(X).
(b) Updating the sequence of tuples f
t+1 g at each iteration of Algorithm 4.1, where
is the parameter vector of Ber( p  t+1 ) having independent components. This is
the same as updating the sequence of tuples f
t+1 g at each iteration of Algorithm
2.1. Note that as soon as the auxiliary discrete distributions is dened, the sequence
f
t+1 g can be viewed as a particular case of the sequence f
t+1 g with
t+1 .
For the case of the partition problem, the generation of random partitions is performed
by a dierent algorithm (Algorithm 3.2 below), but Algorithm 4.1 for updating
the sequence of tuples f
t+1 g is the same for both problems.
We consider separately both phases (a) and (b). More specically, the rest of this
section deals with phase (a), while Section 4 deals with phase (b), where Algorithm 4.1 is
presented.
3.2 Random Cut Generation
The algorithm for generating random cuts in the ASN is based on Ber(p) with independent
components and can be written as follows:
Algorithm 3.1 Random cut generation :
1. Generate an n-dimensional random vector
independent components.
2. From X construct two vectors, V 1 and V 2 , such that the V 1 contains the set of
is a -dimensional vector containing the set of indices corresponding to unities and
the V 2 is a (n ) dimensional vector containing the set of indices corresponding
to zeros. Note that  (0    n) is a random variable.
3. Calculate the sample function M(X) (see also (3.1)), associated with the random
cut
Consider the example in Fig. 3.1. Assume that the cut
is the maximal one. In this case, starting from an arbitrary 6-dimensional vector p with
the goal of Algorithm
4.1 is to converge to the degenerated Bernoulli distribution with the parameter vector
after a nite number of iterations.
Algorithm 3.1 (along with the Algorithm 4.1 below) can be readily extended to randomly
partitioning the nodes V of the graph E) into r  2 subsets and such that
the sum of the total weights of all edges going from one subset to another is maximized.
In this case one can follow the basic steps of Algorithm 3.1 using n r-point distributions
r
instead of n 2-point Ber(p j distributions.
3.3 Random Partition Generation
Here, unlike the independent Bernoulli case, the sample will be generated using a sequence
of m (recall that m is the number of nodes we want in V 1 ) dependent discrete distributions
denoted as
The goal of Fm (p) is to generate an associated random walk of length m, i.e.,
through the nodes of the network (i.e., each i k takes values from the set
such that the nodes visited in the random walk are not repeated. These nodes will then
constitute the set V 1 .
Let
(1)
Clearly,
(1)
is the discrete distribution where the probability of
selecting node i is  (1)
. The node thus selected is denoted by i 1 . The sequence of distributions
will be derived recursively starting from M (1) and will be
used for generating
Algorithm 3.2 :
1. Generate i 1 from the discrete pdf M (1) and set X i 1
1.
2. Derive M (2) from M (1) as follows. First eliminate the element  (1)
from  (1) , and
then normalize the remaining (n 1)-dimensional vector. Call the resulting vector
(2) . Then M
3. Generate i 2 from M (2) and set X i 2
1.
4. Proceed with steps 2 and 3 recursively m 2 times and for each
1. Here i k is generated from M (k) , which is derived from M (k 1) as follows.
Eliminate the element  (k 1)
corresponding to the node i k 1 from  (k 1) , and then
normalize the remaining (n k 1)-dimensional vector. Call the resulting vector
(k) . Then M
5. Set the remaining n m elements of X equal to 0.
6. Calculate the objective function M(X) as in (3.8).
Note that the algorithm does not assume X 1  1. Its modication with X 1  1 is
straightforward: one needs only to replace (3.11) with
(1)
while the rest is similar.
4 The Main Algorithm
We assume below that we are given algorithms for generating random cuts and random
partitions and we are able to calculate the sample function M(X). As mentioned before,
we shall cast the ASN into the rare-events context.
4.1 The Rare-Event Framework
Consider (1.1) for the ASN. Assume for a moment that x is \close" to the unknown true
maximal cut
, which represents the unknown optimal solution of the programs (3.1)
and (3.3). With this in mind we shall adapt below the basic single-iteration and multiple
iteration (see (2.2)-(2.3) and (2.4)-(2.8), respectively) used in Algorithm 2.1 for the ASN
and in particular for the maximal cut, bearing in mind that X  Ber(p). Similar to (2.2)
and (2.3) (the single-iteration program), we have
new
new
and
new
fi:X i;k =1g
respectively. Here it is assumed that the expectation is taken with respect to Ber(p) and
a sample of X.
The optimal solutions of (4.1) and (4.2) can be derived by straightforward application
of the Lagrange multipliers technique. They are
I fM(X)xg X r
I fM(X)xg
and
I fM(X k )xg
I fM(X k )xg
respectively, where expectation in the numerator of (4.3) is the expectation
over all possible cuts for which node V r belongs to V 1
in (4.4) (the sum is over all generated cuts). Also, in (4.4), we set the
new;r to some
arbitrary value, say 1/2, if
I fM(X k )xg
realizing that the chance of the
latter shrinks to 0, as N !1. Let  p
new;n ).
Assume that the optimal solution (V
2 ) of the program (3.1) is unique and consider
the probability '(x) in (1.1). It is obvious that if x >
, then irrespective of the
choice of the parameter vector p in the Bernoulli distribution. We shall present rst an
important observation for the case when
. Let p
d denote a degenerated probability
vector, i.e., it contains a combination of unities and zeros. Moreover, let the components of
dene the unique maximal cut, in the sense that the unity components of p
d correspond
to the components of V
1 and the zero components of p
d correspond to the components of
2 . We shall call p
d the optimal degenerated vector (ODV).
Proposition 4.1 Assume that the maximal cut (V
2 ) is unique. Let X be a random
vector with independent components distributed Ber(p). Then for
the optimal
vector p
new in (4.3) reduces to the ODV p
d irrespective of p, provided
Proof The proof follows immediately from (4.3). To clarify, let X  be the vector from
a degenerated Bernoulli distribution uniquely dening
and let (V
2 ) be the corresponding
cut. Then for any random vector X and the corresponding cut
must have that
I fM(X)
and therefore the rst part of the proposition is proved. For the second part of the
proposition note that 0 <
the second part follows from the fact that (due to similar reasoning as for the rst
new;r in (4.4) is either p
d;r or 1/2 (using our convention), and the fact that
I
It is not di-cult to verify that the variance of the estimate  ' N (x; new ) in (2.1) (note
that v is now replaced by p and similarly for v new ) with
d , equals zero.
As we already mentioned, we shall approximate the unknown true solution (
d ) by the
sequence of tuples f
generated by Algorithm 4.1 below.
4.2 Main Algorithm
To proceed, note again that the single-iteration program (4.2) and its optimal solution
new;r in (4.4) are of little practical use, since for an arbitrary vector p in Ber(p), all
indicators I fM(X j )xg
very high probability, provided x is
close to the optimal value
and the sample size N is small relative to the reciprocal of
the rare-event probability PfM(X)  xg.
To overcome this di-culty we use Algorithm 2.1 where we use
t instead of
t .
(a) Adaptive estimation of
t .
For a xed  p
derive
t as the solution of
I
(b) Adaptive estimation of p
t .
For xed
derive
t;n ) from the solution of
fi:X i;k =1g
Note again that both programs (4.6) and (4.7) can be solved analytically. The solution
of (4.6) is given by (2.6). The solution of (4.7) (see (4.4)) can be written as
I
I
for
I
> 0; as before we set
otherwise, but this case never
happened in the examples we tried.
The resulting algorithm for estimating
and the vector p  is as follows:
Algorithm 4.1 :
1. Choose p in the ASN such that
Generate N random vectors
corresponding cuts using Algorithm 3.1. Deliver the solution (2.6) of the program
(4.6). Denote the initial solution by
2. Use the same N vectors deliver the solution (4.8) of the
program (4.7). Denote the solution by
t .
3. Generate N new random vectors
using Algorithm
3.1 and deliver the solution
t of (2.6) of the program (4.6).
4. If for some t  k and some k, say
stop and deliver
t as an estimate of
. Otherwise, set go to Step 2.
As an alternative to the estimate
t of
and to the stopping rule in (4.9) one can consider
the following:
4  . If for some t  k and some k, say
0st
as an estimate of
. Otherwise, set go to Step 2.
Remark 4.1 Smoothed Probability Vectors Instead of  p
t (see (4.8)) we typically use
its following smoothed version
~
where 0:5 <   1. Clearly, for we have that ~
The reason for using ~
instead of
t;i is twofold: (a) to smooth out the values of
t;i , (b) to reduce the probability
that some values of
t;i will be zeros or unities, specially in the beginning iterations. It
can be readily seen that starting with, say, p
that 0 <
for some indices i  2. We also found empirically that for 0:7    0:9, Algorithm 4.1
is typically more accurate than for in particular for noisy COP. In our numerical
studies we used 0:9. Note that according to the ant-based terminology [12], [15], we
can call  and ~ p
t the evaporation parameter and the pheromone vector, respectively.
Remark 4.2 Relation to Root Finding As mentioned, Algorithm 4.1 might be viewed
as a simple modication of Algorithm 2.1. More precisely, it is similar to Algorithm 2.1 in
the sense of nding the root x (which is associated with the optimal solution
than the rare event probability '(x) by itself. This in turn implies that Algorithm 4.1
involves neither likelihood ratio calculations nor estimation of probabilities of rare events.
For that reason, both Algorithm 2.1 and Algorithm 4.1 have dierent stopping rules.
At this point it may be worth mentioning the following theorem.
Theorem 4.1 Assume that the maximal cut (V
2 ) is unique. If the conditions of
Theorem 2.1 hold, then there exists a t < 1 such that the sequence of tuples f
generated from Algorithm 4.1 converges in distribution to the constant tuple (
d ) as
!1, irrespective of the choice of p, provided
Proof According to Theorem 2.1, setting
, there exists a t < 1 such that
lim
Using (4.8) and a reasoning similar to that of Proposition 4.1, we get that
Combining the two previous facts, we have that there exists t < 1 such that
thus proving the statement of the theorem.
Remark 4.3 To apply this theorem to the maximal cut and partition problems one further
needs to prove that the conditions of Theorem 2.1 hold. For such results in settings
other than the maximal cut and partition problems one is referred to [27] (see, e.g., Proposition
3.1 in [27]).
Remark 4.4 It follows from Theorem 4.1 and Proposition 4.1 that (irrespective of the
initial choice of p), the multiple-iteration procedure involving the sequence f p
converges
to the same degenerated parameter p
d as does the single-iteration procedure.
5 Modications and Enhancements to the Main Algorithm
5.1 Alternative Sample Functions
A natural modication of the two-stage iterative procedure of Algorithm 4.1 would be
to update  p
t (in the second stage; see (4.7)) using some alternative sample functions
rather than the indicators I fM i >
is the abbreviated
notation for M(X i ) that is used, for example, in (4.7).
Consider rst a maximization problem. As for alternatives to I fM i >
one could use
(in (4.7))
1.
where  > 0.
2. Boltzmann type function
exp
where  > 0.
3. Linear loss function with insensitive zone
4. Huber loss function
We found that the above modications typically lead to an increase in the convergence
speed of Algorithm 4.1 two to four times. The reason is that using indicators we put an
equal weight associated with each of the top dNe values of M
(4.8)), while in the modied versions we put a weight proportional to the respective value
of
Consider now a minimization problem. Recall that in this case we use in Algorithm
4.1 the bottom dNe values of M instead of the top ones, since now
we use the indicator function I fM<
. As for a simple modication of the sample function
I fM<
g we could use the top dNe values of 1
MN , (or say the top dNe values
We found, however, that this policy (modication) results typically in a
worse performance of the CE Algorithm 4.1 regardless of . The intuitive explanation is
that the function 1
M  is quite nonlinear and \non symmetric" relative to M  (used in the
maximization). To nd a \symmetric" to the M  function we rst nd a \symmetric"
to the M function as follows: instead of the sample M (N) we can use (in a
minimization problem), say the following sequence [M (N) +M (1)
We can then use the top dNe samples of this new sequence (i.e., raise each element of
this sequence to the power , etc.
One can also use a Boltzmann type function that can be written (in analogy to (5.2))
as
exp
where again  > 0. Similarly for the functions (5.3) and (5.4).
5.2 Single-Stage CE Algorithms versus Two-Stage CE Algorithm
Here the term, single-stage means that at each iteration, the CE algorithm updates  p
alone, i.e., it does not involve the program (4.6) and, thus the sequence
t . In particular,
say in a maximization problem, this would imply updating the vector
t (similar to (4.8))
but taking the entire (rather than the truncated) sample M (or the entire
sample M
N ). For example, using the entire sample M
obtain instead of (4.7) the following stochastic program
fi:X i;k =1g
Such single-stage version would simplify substantially Algorithm 4.1. The disadvantage
of using single-stage sample functions is that, typically, it takes too long for Algorithm
4.1 to converge, since the large number of \not important" (untruncated) trajectories
slow down dramatically the convergence of f p
t g to p
d . We found numerically that the
single-stage CE algorithm is much worse than its two-stage counterpart in the sense that
it is both, less accurate and more time consuming compared to the original two-stage
Algorithm 4.1 (with Practically, we found that it does not work for
maximal cut problems of size n > 30.
Hence it is important for the CE method to use both stages, as in Algorithm 4.1.
This is also one of the major dierences between CE and ant-based methods of Dorigo,
Maniezzo and Colorni [15] and others, where a single-stage sample functions (for updating
t alone) is used.
5.3 Fully Automated CE Algorithm
Here we present a modication of Algorithm 4.1 in which both  and N in (4.6) and (4.7)
are updated adaptively in t. We call this modication of CE, the fully automated CE
Algorithm. In addition, the FACE Algorithm is able to identify some \di-cult"
and pathological problems where it fails. For such problems, the FACE Algorithm stops
if the sample size N t becomes prohibitively large, say there is no improvement
in the objective function.
Let M be the order statistics corresponding to
(that is used while updating the tuple f
t+1 g), but arranged
in decreasing order (note that this is a departure from the convention we used earlier; we
have done this to simplify the notation used below). Notice from (4.6) that
The main assumption in the FACE Algorithm is
for is a xed positive constant in the interval 0:01
c  1. Thus cn corresponds to the number of samples in the set M t;j
lie in the upper 100 t % of the samples. We will refer to the latter as elite samples. Our
second parameter is  as used in (4.11). Note that for c close to 1,  may be chosen close
to unity, say 0.99. In contrast, if we choose should be less, say
We found numerically that combinations, like
and are good choices. If not stated otherwise, we shall use the
combination bear in mind a maximization problem, in particular
the maximal cut problem.
As compared to Algorithm 4.1 we introduce into the FACE Algorithm the following
two enhancements:
1. Prior to the t-th iteration, we keep a portion ; 0    1 of the elite samples
and we incorporate M t (with N and
being replaced by N t and  t , respectively). More precisely, while implementing
(4.6) and (4.7), in addition to the current sample M(X
t , of size
generated from Ber( p
we also include the elite
samples M t from the (t 1)-th iteration. It is readily
seen that for the sequence f
t+1 g in (4.6) and (4.7) is based on the elite
sampling (of size  t N obtained during all t iterations. We use 0:2 in the
numerical experiments reported in this paper, even though other values of  in the
2. Let ~
r be a constant such that ~
r  1. For each iteration t of the FACE Algorithm we
design a sampling plan which ensures that
Note that for ~
implies improvement of the maximal order statistics,
t;(1) , at each iteration. If not stated otherwise, we shall take ~ r = 1.
The rst ~
r iterations (i.e., iteration number 0 to iteration number ~
r 1) of the FACE
Algorithm coincide with that of Algorithm 4.1 where as before we take  p
Also, we typically take  We then proceed
as follows:
Algorithm 5.1 : FACE Algorithm
1. At iteration t;
a sample of size ~
Combine these samples with the elite samples M t
obtained until iteration number t 1, thus giving a total sample
size of N these samples by M t;j ,
M t;j .
2. If (5.8) holds,  proceed with (4.6)-(4.7) using the N t samples mentioned
in Step 1.
3. If (5.8) is violated, check whether or not
holds, where, say
stop and deliver M t;(1) as an estimate of the optimal solution. We call such M t;(1)
a reliable estimate of the optimal solution. Otherwise
4. Increase the sample size ~
holds or ~
very
large, say N In the former case, set N
and proceed again with (4.6) and (4.7). In the latter case, stop and deliver M t 1;(1)
as an estimate of the optimal solution. We call such M t 1;(1) , an unreliable estimate
of the optimal solution.
Remark 5.1 In Step 3, we use c in the numerical experiments reported in this
paper, even though other values of c 1 in the given range gave similar results.
Remark 5.2 To save sampling eort, we terminate Step 4 in a slightly dierent manner
than as stated in Algorithm 5.1. Let N  be such that N 0 < N  < N  , say N
(implicitly, we are assuming that 20n < 10 6 which is usually the case in practice; also,
for SEN networks we take N In any iteration t, while increasing ~
N t as given
in Step 4, if we obtain that ~
and (5.8) is still violated, then we interrupt Step
4, and directly proceed with updating (
do (4.6) and (4.7) (before proceeding with the next iteration). However,
if with this slight change in Step 4, Algorithm 5.1 keeps generating samples of size N
for several iterations in turn, say for then we complete Step 4, as given
in Algorithm 5.1.
Remark 5.3 Zig-zag policy Let ~ k be a given constant such that ~ k  1. If for some t
then we can decrease N 0 by some factor starting from iteration t + 1. Similarly, if (5.8)
does not hold for ~ k consecutive iterations, then we can increase N 0 by some factor. Note
that we generate at least ~ k after the latest increase or decrease in N 0 , before
we start to check (5.11) or start to check the (5.8)'s (for the previous ~ k iterations). In the
numerical experiments reported in this paper we use ~ we increase or decrease
by a factor of 2.
Note that Algorithm 5.1 uses a maximal order-statistics based stopping criterion given
by (5.9). This seems more natural compared to the quantile based stopping criterion used
in Algorithm 4.1. Note also that (5.9) is based on the fact that while approaching the
degenerating solution
d , more and more trajectories will follow the path associated with
d and thus, the number of dierent trajectories will be less than ~
Note nally that if (5.8) holds for all t  ~ r we automatically obtain that N
8t. In that case Algorithm 5.1 reduces to the original Algorithm 4.1, provided
6 Numerical Results
Below we present the performance of Algorithm 4.1 and Algorithm 5.1 for both the maximal
cut and partition problems. By the performance of the algorithms, we mean the
convergence of estimators
t and
t  (see, e.g., (4.9), (4.10)) to the true unknown optimal
value
. As mentioned before, we choose in the respective ASNs.
Since the maximal cut and partition are NP hard problems, no exact method is known
for verifying the accuracy of our method except for the naive total enumeration routine,
which is feasible for small graphs, say for those with n  nodes. To overcome this
di-culty, we construct an articial graph such that the solution is available in advance
and then verify the accuracy of our method. As an example, we consider the following
symmetric distance matrix
| {z }
| {z }
where all components of the upper left-hand and lower right-hand quadrants (of sizes
m)  (n m), respectively, where 0 < m < n) are equal and generated
from a given distribution, such as U(a; b) (uniformly distributed on the interval (a; b)),
etc., and the remaining components equal We choose
C 2 such that the partition V
will be the optimal one. More precisely, assume that the random variable Z has a pdf
with bounded support, say Z  C 1 . Let, for example, clearly for
2 ) to be the optimal solution, it su-ces that C 2 > C 1 (similarly for m < n=2).
In all our examples we set  = 0:01 in (4.6), took stopped
Algorithm 4.1 according to the stopping rule (4.9) with the parameter 5. We found
that for r  10, the relative error, dened as
equalled zero in all our experiments. Here T corresponds to the stopping time of Algorithm
4.1. The running time in seconds of Algorithm 4.1 on a Sun Enterprise 4000 workstation
(12 CPU, 248 MHz) is reported as well.

Table

6.1 presents the relative errors  (and the associated stopping times
as function of the sample size for the maximal cut problem
with for the following 6 cases of Z in (6.1): Z
Here Beta(; ; a; b) is the Beta distribution whose probability density function is given
by
()() (b a

Table

6.1 The relative errors  (and the associated stopping times
as functions of the sample size for the maximal cut problems
with for the above 6 cases of Z given in (6.1).

Table

6.2 The CPU times T as functions of the sample size
the same data as in Table 6.1
The results of Tables 6.1 - 6.2 are self-explanatory. Similar results were obtained with
Algorithm 4.1 for the partition problems.
The rest of our numerical results are for the partition problems. Tables 6.3 - 6.4
presents the performance of
t along with as functions of t, where
are dened as
t;s  0:5; ng
and
t;s < 0:5;
In particular Table 6.3 corresponds to

Table

6.4 presents data similar to Table 6.3 for In both
cases we obtained a relative error of 0 with CPU times
respectively. The results of Tables 6.3 - 6.4 are self-explanatory.

Table

6.3 Performance of Algorithm 4.1 for
Z  U(4:5; 5).
6 1.220356e+06 0.510000 0.490000
9 1.224679e+06 0.510000 0.490000
13 1.238418e+06 0.510000 0.470000
14 1.240769e+06 0.590000 0.420000
19 1.249318e+06 0.890000 0.090000

Table

6.4 Performance of Algorithm 4.1 for
Z  U(4:5; 5).
6 1.308071e+06 0.516667 0.483333
9 1.410264e+06 0.833333 0.466667
14 1.456611e+06 1.000000 0.233333
19 1.458000e+06 1.000000 0.000000

Table

6.5 presents the dynamics of  p
t;14 ) for another smaller example.

Table

6.5 Dynamics of
t for for
3 1:00 0:97 0:95 0:99 0:87 0:99 0:98 0:01 0:03 0:00 0:11 0:00 0:02 0:00
4 1:00 1:00 1:00 0:99 1:00 0:99 0:98 0:01 0:00 0:00 0:00 0:00 0:01 0:00
5 1:00 1:00 1:00 0:99 1:00 0:99 0:99 0:01 0:00 0:00 0:00 0:00 0:00 0:00

Table

6.6 presents the performance of Algorithm 5.1 for the same input data as Table
6.4 with N that M t;(N) denotes the best (elite) sample value
of M(X) obtained at iteration t. We found that Algorithm 5.1 is at least as accurate as
Algorithm 4.1 and is typically 2-3 times faster than Algorithm 4.1.

Table

6.6 Performance of Algorithm 5.1 for the same input data as Table 6.4 with
9 1.42e+06 1.39e+06 3000 0.01 0.87 0.43
At this point we would like to note that we have performed extensive simulations case
studies with Algorithm 4.1 for dierent SNN models. We found that in approximately
99% of the cases, Algorithm 4.1 performs well in the sense that the relative error  dened
in (6.3) does not exceed 1%. These results will be reported somewhere else.
6.1 Empirical Computational Complexity
Let us nally discuss the computational complexity of Algorithm 4.1 for the maximal cut
and the partition problems, which can be dened as
Here T n is the total number of iterations needed before Algorithm 4.1 stops; N n is the
sample size, that is the total number of maximal cuts and partitions generated at each
iteration; G n is the cost of generating from n independent Bernoulli distributions (needed
in Algorithm 3.1) or from the sequence of distributions Fm (p) (needed in Algorithm 3.2);
is the cost of updating the tuple (
). The latter follows from the
fact that computing M(X) in (3.8) is a O(n 2 ) operation.
For the model in (6.1) we found empirically that T
1000. For the maximal cut problem, considering that we take n  N n  10n and that G n
is O(n) , we obtain  In our experiments, the complexity we observed was
more like
The partitition problem has similar computational characteristics.
It is important to note that these empirical complexity results are solely for the model
with the distance matrix (6.1).
7 Concluding Remarks and Directions for Further Research
This paper presents an application of the cross-entropy (CE) method [38] to the maximal
cut and the partition problems. The proposed algorithm employs an auxiliary discrete dis-
tribution, which transforms the original deterministic network into an associated stochastic
one, called the associated stochastic network (ASN). Each iteration of the CE method
involves two major steps: (a) generation of trajectories (cuts and partitions) using an
auxiliary discrete distribution with a parameter vector p and calculation of the associated
objective function M(X) and some related quantities, such as indicator functions associated
with rare-event probabilities, and (b) updating the parameter vector p on the basis
of the data collected in the rst step. Our numerical studies with the maximal cut and the
partition problem, as well as some other COPs (not reported in this paper) suggest that
the proposed algorithms typically perform well in the sense that in approximately 99% of
the cases the relative error  dened in (6.3) does not exceed 1%. In addition, experiments
suggest that Algorithm 4.1 and Algorithm 5.1 have polynomial complexity in the size of
the network. Some further topics for investigation are listed below.
1. Establish convergence of Algorithm 4.1 for nite sampling (i.e., N < 1) with emphasis
on the complexity and the speed of convergence under the suggested stopping
rules.
2. Establish condence intervals (regions) for the optimal solution.
3. Apply the proposed methodology to a wide variety of combinatorial optimization
problems, such as the quadratic assignment problem, minimal cut, vehicle routing,
graph coloring and communication networks optimization problems.
4. Apply the proposed methodology to noisy (simulation based) networks. Our preliminary
studies suggest that Algorithm 4.1 performs well in the case of noisy networks.
5. Apply parallel optimization techniques to the proposed methodology.

Acknowledgments

I would like to thank Alexander Podgaetsky for performing the numerical
experiments. I would also like to thank three anonymous referees, Pieter-Tjerk
de Boer and Victor Nicola (Guest Editor) at University of Twente, The Netherlands, and
Perwez Shahabuddin (Guest Editor) at Columbia University, U.S.A., for many valuable
suggestions.



--R

John Wiley
Local search in combinatorial optimization
Network ows theory













Genetic algorithms in search





Introduction to global optimization










Global optimization in action

Modern heuristic search methods



Discrete event systems: Sensitivity analysis and stochastic optimization via the score function method







--TR
Discrete optimization
Simulated annealing and Boltzmann machines: a stochastic approach to combinatorial optimization and neural computing
Approximating the permanent
Network flows
Polynomial-time approximation algorithms for the Ising model
search
A method for discrete stochastic optimization
Ant-based load balancing in telecommunications networks
New parallel randomized algorithms for the traveling salesman problem
The ant colony optimization meta-heuristic
ACO algorithms for the quadratic assignment problem
Ant algorithms for discrete optimization
A Graph-based Ant system and its convergence
Genetic Algorithms in Search, Optimization and Machine Learning
Local Search in Combinatorial Optimization
Simulated Annealing
Efficiently searching a graph by a smell-oriented vertex process
Nested Partitions Method for Global Optimization

--CTR
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient importance sampling heuristics for the simulation of population overflow in Jackson networks, Proceedings of the 37th conference on Winter simulation, December 04-07, 2005, Orlando, Florida
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient simulation of population overflow in parallel queues, Proceedings of the 37th conference on Winter simulation, December 03-06, 2006, Monterey, California
P. T. de Boer , D. P. Kroese , R. Y. Rubinstein, Rare event simulation and combinatorial optimization using cross entropy: estimating buffer overflows in three stages using cross-entropy, Proceedings of the 34th conference on Winter simulation: exploring new frontiers, December 08-11, 2002, San Diego, California
Zdravko Botev , Dirk P. Kroese, Global likelihood optimization via the cross-entropy method with an application to mixture models, Proceedings of the 36th conference on Winter simulation, December 05-08, 2004, Washington, D.C.
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient heuristics for the simulation of population overflow in series and parallel queues, Proceedings of the 1st international conference on Performance evaluation methodolgies and tools, October 11-13, 2006, Pisa, Italy
Victor F. Nicola , Tatiana S. Zaburnenko, Efficient importance sampling heuristics for the simulation of population overflow in Jackson networks, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.17 n.2, p.10-es, April 2007
cross-entropy;importance sampling;rare event simulation;combinatorial optimization