Replica Exchange Using q-Gaussian Swarm Quantum Particle Intelligence Method
Faculty of Engineering, Department of Computer Engineering, International Balkan University, Skopje, R. of Macedonia
To cite this article:
Hiqmet Kamberaj. Replica Exchange Using q-Gaussian Swarm Quantum Particle Intelligence Method. Engineering and Applied Sciences. Vol. 1, No. 2, 2016, pp. 20-25. doi: 10.11648/j.eas.20160102.12
Received: June 15, 2015; Accepted: July 1, 2015; Published: July 16, 2016
Abstract: We present a newly developed Replica Exchange algorithm usingq -Gaussian Swarm Quantum Particle Optimization (REX@q-GSQPO) method for solving the problem of finding the global optimum. The basis of the algorithm is to run multiple copies of independent swarms at different values of q parameter. Based on an energy criterion, chosen to satisfy the detailed balance, we are swapping the particle coordinates of neighboring swarms at regular iteration intervals. The swarm replicas with high q values are characterized by high diversity of particles allowing escaping local minima faster, while the low q replicas, characterized by low diversity of particles, are used to sample more efficiently the local basins. We compared the new algorithm with the standard Gaussian Swarm Quantum Particle Optimization (GSQPO) and q-Gaussian Swarm Quantum Particle Optimization (q-GSQPO) algorithms, and found that the new algorithm is more robust in terms of the number of fitness function calls, and more efficient in terms of ability to convergence faster to the global minimum. In additional, we also provide a method for optimally allocating the swarm replicas among different q values. Our algorithm is tested for three benchmark functions, which are known to be multimodal problems, at different dimensionalities.
Keywords: Swarm Quantum Particle, q-Gaussian Distribution, Global Optimization, Replica Exchange
The problem of finding the global optimum in a multimodal and multidimensional space can be extremely difficult since the number of stable optima increases as the search space increases, for instance, the search for the global minimum energy in a surface energy landscape of the atomic structures. [1,2] Swarm Particle Optimization (SPO) is a population-based optimization technique, similar to evolutionary algorithms.  Kennedy & Eberhart introduced the method to solve the problem of finding the global optimum of a q dimensional function.  The SPO method is based on the swarm intelligence algorithms, which concern with the design of intelligent multi-agent systems based on the collective behavior of insects (ants, termites, bees, and wasps) or other animal societies (flocks of birds and schools of fish). 
In SPO method, the swarm particles, representing possible solutions, search the phase space, defined by their velocities and coordinates, which are updated based on the particle’s own experience and experience of the particle’s neighbors or the experience of the whole swarm. The method has already been used to solve many optimization problems,  with some interests also in other fields, such as statistical mechanics. 
Since the standard SPO algorithm has a low convergence rate, [7,8] several improvements and variants of the SPO algorithm have been proposed. [9,10,11,12,13] The new variant of the SPO method, the so-called Swarm Quantum Particle Optimization (SQPO), has been considered as an improvement versus the classical SPO method, since there is a nonzero probability to escape the local minima even for very high barriers.  Efforts have been made to improve the SQPO method. [15,16,17,18,19,20,21,22] These improvements focus primarily on parameter selection criteria, [18,22] and maintaining diversity of the swarm. [19,20,21] A detailed review of all these methods is described in Ref.  Use of different forms of attractive potential-energy surfaces for SQPO algorithm is also considered for improvement of the algorithm.  Different potentials yield different probability distributions, which describe the probability of finding the swarm quantum-like particle at a certain position in the phase space. 
In our previous study , we showed that the use of q-Gaussian probability distribution of swarm particles (q-GSQPO algorithm) improved significantly the efficiency of searching for the global minimum when compared with Gaussian probability distribution of swarm particles (GSQPO algorithm). The q-Gaussian distribution is characterized by long tails, which allow reaching long distant regions in phase space by increasing the diversity of the swarm particles.  The application of the probability distributions with heavy tails is found to be useful in allowing the system escaping from local optima in multimodal problems also in other areas [25,26,27,28,29], since the nonzero occupations taken from long tails of the distribution implies jumps of scale-free sizes. The q - Gaussian distribution is also known as Tsallis distribution  in statistical mechanics, and it is known to produce a smoothed potential energy surface. 
The problem facing GSQPO algorithm is the premature convergence to a local minima due to low diversity of the swarm particles.  This problem of GSQPO algorithm can efficiently be solved by using the q-GSQPO algorithm , which increases the diversity between the swarm particles, and hence allowing the swarm to explore long distant regions of the searching space.
In this paper we will explore the use of the popular replica exchange method [33,34], which is known for overcoming the problem of the sampling convergence in dynamical systems. In this method several copies of the system (so-called replicas) are run independently at different values of some internal parameter of the system (e.g., temperature [34,35,36], strength of interaction , and q parameter when combined with Tsallis statistics [27,38]). At regular time intervals the coordinates of the replicas are swapped based on an energy criterion that satisfies the detailed balance. The higher parameter value replicas are used to enhance the barrier crossings and the low parameter value replicas are used to sample the local basins.
Similarly, in this study, we will create M copies of swarm particles (replicas) and run them independently at different values of q, where the lowest level replica is running at q=1. Then, in analogy with standard replica exchange methods, at regular iteration intervals, we swap the particles of swarm replicas between two neighboring values of q based on an energy criterion, which is explained in details in the next section. To increase the efficiency of the algorithm, we also discuss the optimization of the round-trip time of the swarm replicas from the lowest to the highest value of q and vice-versa. To test our algorithm, we studied three benchmark functions, namely the Ackley , Griewank  and Rastrigin  functions at different dimensions, d =5, 10, 20 and 50.
2. Materials and Methods
2.1. Theoretical Basics of Generalized q-GSQPO Algorithm
The swarm quantum particle optimization algorithm determines the probability of finding the swarm particle at the position at any time t.  Details of the derivation of GSQPO and q-GSQPO algorithms can be found in our previous work.  Here, we re-write the equations, proposed to describe the q-GSQPO algorithm in a generalized form as :
where has to be less than 1.7 in order to guarantee the convergence of the particles,  and it is chosen here to have a simplified sinusoidal expression as
with g being a scaling constant fixed in this study to 0.5, , and A was varying from 0.01 to 1. , also called Mainstream Thought or Mean Best of the population , defines the mean of the best positions of all particles (i=1, 2, …, N) in each dimension and it is given by 
Here, Fq(u) is a random number following q- Gaussian distribution if and Gaussian distribution if q=1.
2.2. Replica Exchange Method
In the Replica Exchange q-GSQPO (REX@q-GSQPO) algorithm, we replicate M copies of the swarm particles among different values of q: q1, q2,…, qM. A geometrical distribution in the interval [1, qmax) was used to select the values of q, where qmax is a maximum value of q, which is chosen here to be between 2 and 3. In order to control the acceptance probability of swapping swarm particles between two neighboring q, we introduce a "temperature" parameter at each level:
where i=1,2,…,M. k is considered an adjustable parameter, which will be optimized to minimize the time of a round trip from the lowest to the highest value of q and vice-versa. At regular iteration steps we swap the configurations of swarm particles between two randomly chosen neighboring q, let say i and j, such that the detailed balance is satisfied:
where and are the transition probabilities from the swarm replica i to j and j to i, respectively. We will assume that the probability of finding the swarm i at a value q is proportional to the Boltzmann factor as:
where Ei is the minimum of the best local scoring value of the swarm particles in replica i at a given iteration step. Choosing Ei as the minimum of the best local scoring value is optional. Other choices may also be allowed, for instance, the mean best local scoring value or the best global scoring value of the replica. It is important that the choice should yield a high value of the diversity of swarm particles on each of replica in order to allow the swarm particles exploring more efficiently the searching space. This is in particular important for low q replicas, which are used to explore more efficiently the local basins.
Then, the probability of accepting an attempting swap between swarm replicas i and j is given by:
The algorithm stops for searching the phase space when at least one of the swarm replicas converges to the optimal solution. In general, the global minimum of the problem is not known; therefore the searching does not have to stop when one of the replicas converges. For instance, one can restart the replica that has converged randomly from a new position. However, we found, for the problems under investigation, stopping the phase space search after the observation of the first replica converging to some optimal solution is one efficient way for determining criteria for stopping the algorithm.
3. Results and Discussions
We considered the motion in three d-dimensional potential functions determined mathematically as in Table 1.
In order to estimate the efficiency of the algorithm, we calculated the frequencies of visiting each q from each swarm replica and compared the observed frequencies for each swarm with expected frequencies, which depend only on the number of replicas. We used the test with a confidence level of 95 % to investigate the goodness of the comparison.  In this study we chose a geometrical distribution of q among 5 or 6 replicas. The replicated exchange simulations were stopped when the first global best score of the swarms was less than some minimum value, chosen here to be 10-5. As an estimation of the computational robustness of the method we used the number of iterations needed until the convergence was reached.
In addition, we also measured the highest and the lowest values of the diversity, which is given by 
3.1. Adjusting Temperature Parameter k
Value of k will influence on the average allocation of the swarm replicas over the q space ladder, since the overlapping of probability distributions between two neighboring q values, and hence the acceptance probability of swapping the swarm replicas, depend on k. Due to the barriers that exist in the q space, the swarms at the low q values may not be able to visit the upper values and vice-versa. In order to adjust the value of k we did several test runs, and compared the computed frequencies of observing the swarms at a value of q with the expected frequencies assuming equal probability of finding a swarm at each value of q.
In Figure 1 we show the ratio of the computed with critical value of at confidence level of 95 % for the three benchmark functions at different values of parameter k. The dimension the searching space was fixed to d=10. For, the probability distribution of swarm replicas in the q space is considered to be uniform, at a confidence level of 95 %, otherwise there is not enough evidence to say that this distribution is uniform. We found that for values of , depending also on the benchmark function, there is enough evidence to say at the confidence level of 95 %, that the distribution of the swarm replicas in the q space is uniform (see Figure 1).
Figure 1. The ratio between the calculated and critical value of at the confidence level of 95 % for different values of the parameter k. (A) For Ackley function; (B) Griewank function and (C) Rastrigin function. Dimensionality of the space was fixed to d=10, and q was chosen from a geometrical distribution in the range between one and two.
In order to examine the robustness of the REX@q-GSQPO algorithm in comparison with other algorithms such as GSQPO and q-GSQPO in terms of function calls, we studied the number of iterations needed for each algorithm until the convergence was reached for a fixed value of . For the REX@q-GSQPO algorithm we replicated five swarms among five values of q in the range from 1 to 2 chosen from a geometrical distribution. The swapping of the swarm particle configurations between two neighboring values of q was performed at each iteration step. For the q-GSQPO algorithm we chose these values . The results of the number of iterations for the GSQPO (q=1) and q-GSQPO algorithms were averaged over different values of .
Our results are presented graphically in Figure 2 for the three benchmark functions versus the search space dimension d (Ackley (A), Griewank (B) and Rastrigin (C)): in gray the average values for GSQPO and q-GSQPO algorithms, and in black for REX@q-GSQPO algorithm. Our data indicate that REX@q-GSQPO algorithm is more robust than GSQPO and q-GSQPO algorithms since the number of iterations needed for reaching the convergence is much smaller compare to GSQPO and q-GSQPO algorithms, in particular for high dimensionality which is of interest for many applications.
To examine the efficiency of the method on finding the global minimum, we investigated the average best score and compared the three algorithms presented here. We used the same setup as described above. In Figure 3 we are plotting the average best score for the REX@q-GSQPO algorithm (in black) and an average value for the GSQPO (q=1) and q-GSQPO () algorithms (in gray). Results are presented for the three benchmark functions and different dimensions of the searching space d. It can be seen that REX@q-GSQPO algorithm provides much smaller best scoring values compare to standard GSQPO and q-GSQPO algorithms. In addition, for the REX@q-GSQPO algorithm, in contrast to GSQPO and q-GSQPO algorithms, the convergence to the global minimum was achieved for all benchmark functions and for all dimensions considered here (see Figure 3).
In this paper we presented a new algorithm for determining the global minimum of a multimodal problem using the replica exchange approach of the swarm quantum particles over the q searching space. The algorithm was compared with standard Gaussian swarm quantum particles and q-Gaussian swarm quantum particles algorithms in terms of the efficiency of convergence to the global minimum and computational robustness.
It was found the new algorithm, REX@q-GSQPO, outperforms the standard algorithms GSQPO and q-GSQPO in both, the efficiency of the convergence and the computational robustness. In addition, we also examine and showed how to optimize the probability distribution among q values in order to minimize the round-trip of the swarm replicas between the two extreme q values. This was an important step to ensure the ergodicity in the q space.
In the REX@q-GSQPO algorithm, the high q values swarm replicas, which are characterized by high diversity between the particles, can be used to explore distant regions in the sampling space, while the low q values, which are characterized by low diversity between the particles, can be used to sample efficiently the local basins.
We envision that the approach proposed in this study can be applied to other examples, as well, of swarm formation during collective migration including collective behaviour of biomolecular dynamics  by improving searching mechanism as proposed very recently  using the principle component analysis.
The author would like to thank the International Balkan University for the support.