Stability analysis based parameter tuning of Social Group Optimization

Swarm-based optimization algorithms have been popularly used these days for optimization of various real world problems but sometimes it becomes hard to estimate the associated characteristics due to their stochastic nature. To ensure a steady performance of these techniques, it is essential to have knowledge about the range of variables, in which a particular algorithm always provides stable performance and performing stability analysis of an algorithm can help in providing some knowledge regarding the same. Many researchers have performed the stability analysis of several optimization algorithms and analyzed their behavior. Social Group Optimization (SGO) is a newly developed algorithm which has been proven to yield promising results when applied to many real world problems but in literature no work can be found on stability analysis of SGO. In this paper, Von Neumann stability analysis approach has been used for performing stability analysis of Social Group Optimization (SGO) to analyze the behavior of its algorithmic parameters and estimate the range in which they always give stable convergence. The results obtained have been supported by sufficient experimental analysis and simulated using eight benchmark function suite along with their shifted and rotated variations which prove that the algorithm performs better within the stable range and hence convergence is ensured.


Introduction
In many real life applications, many complex optimization problems are faced which needs effective solutions. These functions have multiple maxima, minima and saddle points, thus gradient descent/ascent based techniques fail to find the proper solution or most likely lead to a local minima or maxima. For finding the global best value efficient searching of the entire search space becomes necessary. There are several categories of algorithm that fulfill this need. The prominent of them is swarm optimization in which the search space is filled with swarm members with random initial solutions. These members interact with each other (defined by the algorithm), through a series of exploration and  [1] which emulates flocking behavior of birds, Artificial Bee Colony (ABC) [2] which emulates behavior of bees, Social Group Optimization (SGO) [3] which emulates human learning pattern of, Grey Wolf Algorithm (GWA) [5] which emulates hunting patterns of Grey Wolves etc. Social Group Optimization (SGO) is a fairly new addition to the shelve. Nonetheless it is very promising addition as it has better convergence rate as compared to many evolutionary algorithms [4]. It has been successfully applied to segmentation of skin melanoma images [6], segmentation of brain MRI images [7], task scheduling in cloud [8], brain tumor evaluation tool [9], automated detection of COVID-19 infection [10], Transformer fault analysis [16], Antenna array synthesis [17], short-term hydrothermal scheduling [18], structural health monitoring in civil engineering [19], solving Travelling Salesman problem [20], solving multiobjective problems [21] and many more.
One essential requirement to be kept in mind while working with these heuristic algorithms is to ensure that the error values remain bounded in a region and don't explode. To achieve this requirement, stability analysis of an optimiza-tion algorithm is essential as it provides the stable range of parameters that prevents the errors from growing beyond boundary making the algorithm converge faster to provide desired solution. In many of the works found in literature, it is observed that the efficacy of an algorithm is claimed by experimentation, but that experimental analysis is subject to that specific problem or situation. Stability analysis is such a method which analyzes the robustness of a system and 'Von Neumann' method is a suitable method which have been popularly used for this task for swarm optimization methodologies. Various such works can be found in the literature. Gopal et al. [11] have worked on stability analysis of PSO, Bansal et al. [12] in their work performed stability analysis of Differential Evolution, Nair et al. [13] analyzed the stability of Artificial Bee Colony, Farivar et al. [14] have reported on stability analysis of Gravitational Search Optimization, Biswas et al. [15] have analyzed the stability of Bacteria Foraging Optimization etc. To the best of authors' knowledge the stability analysis of SGO algorithm is not yet attempted. In this paper, authors have addressed the problem of stability of SGO through 'Von Neumann' stability criterion and have found stable ranges of parameters for the algorithm.
The rest of the paper is organized as follows. In the second section introduces the SGO algorithm and effects of its update equations on swarm members, in the third section provides details of the Von Neumann Stability Analysis procedure and explains stability analysis in context of swarm optimization. The mathematical analysis of SGO algorithm is carried out in the fourth section and the final section presents the results of simulation based experiments to demonstrate the simulation based find ups.

Social Group Optimization
Social Group Optimization (SGO) algorithm is a population based distributed optimization algorithm that simulates human learning pattern to iteratively search and reach the global maximum or minimum. It has two distinct phases i.e. improving phase and acquiring phase that execute sequentially in each iteration to exploit and explore the search region respectively.

Initialization
The initial position of the each member in the population is generated by Here A i denotes position of the ith member in the D dimensional space. A i, j denotes the jth dimension of the position vector. A j max A j min are the upper and the lower bounds of jth parameter. N is the population size and D is the dimensionality of the optimization problem under consideration.

Improving phase
In this phase each swarm member X i is updated as per Eq. (2.2), Here c is the self introspection factor and r is a random number. Then A new i, j is compared with A i, j and the one who provides better performance is accepted as value of A i, j . This equation has a clumping effect on the swarm members as the values of A i, j will be remapped randomly in the connected areas defined by the 4 line segments given in Eqs.
(2.3a) to (2.3d) where r ∈ [r n , r x ]. Eqs. (2.3a) and (2.3b) represent the update equation in improving phase (Eq. (2.2)) when value of r is r x and r n respectively. Equations (2.3c) and (2.3d) represent the boundary condition x min and x max defined by the optimization problem.
y = x min , and y = x max (due to boundary constraint) (2.3d) If we observe Eqs. (2.3a) and (2.3b), they intersect at the point c * gbest as shown in Fig. 1a. Thus, c * gbest is the cluster centroid. The cluster centroid is formed on the line connecting the origin and gbest. So smaller c value always forms the cluster closer to the origin. From Fig. 1a it is clear that for any value of X i in phe the range, there is always a possibility that it gets mapped to c * gbest (Can be verified by taking a vertical cross section of the grey area shown in Fig. 1a) and the farther a point is from the centroid on yaxis, the range of values of X i which gets mapped to the value decreases. As the vertical cross-section of the probable area (Grayed portion in Fig. 1a) indicates probability of y falling in the region, most points will clump around the cluster centroid. In Fig. 1, the graphical representation of transformation done by Improving phase is shown. The Fig.  1a shows the area defined by the Eqs. (2.3a)-(2.3d) and how the range of X i is reduced along each dimension(projection on x-axis shows the initial range, and projection on y-axis shows the range after improving phase). The Fig. 1b shows a 2D visualization of clustering of X i after improving phase.
The spread of the cluster roughly can be determined by the angle between the lines reddenoted by φ. It is defined by Hence the angle can be found as

Acquiring phase
In acquiring phase, each swarm member A i moves in the search space with respect to gbest and another randomly selected swarm member A s . The update equation is given by Here A s is a member of the swarm s.t. i = s and r 1 , r 2 ≥ 0. The update equations in acquiring phase can be written as where A n = r 2 * gbest ± r 1 * A s r 2 ± r 1 and r n = r 2 ± r 1 So A i updates it's position with respect to a vector which is weighted mean of A s and the gbest. Figure 2 shows the direction of movement of A i in acquiring phase. So every swarm member updates itself by moving towards or away from another swarm member while moving itself slightly closer to gbest.

Stability analysis
The stability analysis is an important requirement of numerical estimation schemes involving partial differential equations. In the pre-text context the stability reflects that the total error of a numeric scheme remains bounded. In other words, if some error is introduced it would remain under a bound and would not explode or increase beyond a limit. The Von Neumann stability analysis, which uses Fourier Decomposition, is an important technique, which is chosen to analyse the stability of a scheme.

von Neumann stability analysis
Consider the partial differential equation represents a system X , a function of χ and τ , which evolves through space (χ ) and time(τ ). The numeric estimate of the corresponding partial differential equation is given by red Equation (3.2) shows numeric estimate of Eq. (3.1) The Fourier series expansion of X (τ, χ ) can be written as where χ = i χ , τ = t τ . χ and τ represent the intervals that were used to sample values along the χ and τ axis. LetX (t, i) represent an estimate of X (τ, χ ) with an error t,i as shown in Eq. (3.4).
Here we want to learn whether the error associated with the estimate grows or shrinks with each iteration. So we assume that the error is associated with some components of the Fourier expansion of X . We replace X in the numeric estimate is replaced by its fourier component to check whether contribution of a single fourier component increases or diminishes over time. If it diminishes then the error associated with the estimate decreases with time too. If the relation between any Fourier coefficient in tth iteration with the corresponding Fourier coefficient in (t + 1)th iteration can be represented by the Eq. (3.5) and |g(α m )| ≤ 1 then contribution of each Fourier component decreases with each iteration. So any error associated with it will also diminish over time and we will conclude that the numeric scheme is stable.

Stability Analysis in context of swarm optimization
In swarm optimization algorithms the update equation is the form Here A r is another member of the swarm besides A i . A i interacts with another member of the swarm to update it's value over each iteration. To determine if it is stable we would want to know if some error gets introduced in an iteration then it diminishes or increases with each subsequent iteration.
For the stability analysis A is modelled as a continuous variable across space and time, where each iteration and each member of the swarm is taken as samples at τ = t τ and χ = i χ . i.e.   Table 1 lists eight standard objective functions to be optimized along with their shifted and rotated version. For this the steps mentioned in the CEC 2014 [22] is followed. The functions are mentioned under the column objective function of Table 1 and are expressed as a function of z. But for better control and more variety we define z as a function of x. So we can shift,scale and rotate x vectors before passing them into the functions. The mappings from x to z are defined under the column z. Different mappings from x to z are defined to generalize the properties as stated-

Optimization functions
Global Best value location: gbest = o and M is a DxD rotation matrix. In this paper we will keep R = 100 as per CEC 2014 [22] and we will use different values of M,o and D to test the algorithm.

Approach for a solution
In SGO algorithm we have two phases that operate sequentially in a single epoch which gives the updated value of A new i for the next epoch. The equations governing the improving phase is, and acquiring phase is, A acq i the value of A new i for the next iteration. The update equation can be written as is a four dimensional equation w.r.t to the parameters c, r , r 1 , r 2 . This equation is very difficult to solve as von Neumann analysis won't provide enough inequalities to perfectly define ranges of all the parameters. Moreover we will get equations where range of c, r is dependent on r 1 , r 2 and vice versa, which is not desirable as both improving and acquiring phases are independent of each other. But we can use a simple trick to simplify these issues. In Eq.
i, j and replace A t i, j and A t+1 i, j with the corresponding Fourier component. We would get an equation of the form, Doing the same for acquiring phase we would get, Fourier representation of Eq. (4.3) can be expressed as, and the corresponding condition for stability is are satisfied, then the condition of stability will be automatically satisfied. Therefore the stability analysis of improving phase and acquiring phase independently can be taken up to obtain the ranges for c, r , r 1 , r 2 .

Improving phase
The Improving Phase has 2 parameters that govern the equation, c and r . c is the self introspection parameter which is constant, i.e. its value stays the same throughout the algorithm. r is a random number. The update equation is given in Eq. (4.7) By replacing A new i, j and A old i, j with A t+1 i, j and A t i, j respectively, Eq. (4.7) can be rewritten as Eq. (4.8).
In the update equation the dimensions are independent of each other, we can drop j and assume it to be an one dimensional problem without loss of generality. So the Eq. (4.8) reduces to Eq. (4.9).
Here the superscript t denotes the value at tth iteration. For simplifying our calculations we can take gbest as a constant. This decision is justified by the fact that after some initial iterations the value of gbest updates only occasionally through the algorithm run.
As gbest is a constant we can remove the term from the calculations for stability analysis as there are no errors associated with a constant value. Mathematically, the mth Fourier component of a constant term is 0 as the Fourier series of the constant is represented by the constant itself.
Replacing A t i with it's mth Fourier component A m e −Jβ m tτ e Jα m iχ (where J is the imaginary unit) in the Eq. (4.9) an equation similar to the form Eq. (3.5) is obtained.
As per von neumann stability criterion(given in Eq. (3.8)) |e Jβ m τ | ≤ 1. Upon solving the inequality we get, We obtain the range of r as defined in Eq. (4.10). The range of r is dependent on c.
Putting these values defined in Eqs. (4.11a)-(4.11c) in Eq. (2.4), we get Therefore, if we use the complete stable range for improving phase, the angle of spread has to be π/2.

Experimental results
To verify the findings, the acquiring phase is removed from SGO algorithm and only the improving phase is considered. The stable range of r against the three schemes of unstable ranges is compared. They are Rastrigin, Rosenbrock and Alpine function are used for the testing with D = 30, number of swarm members i.e. N = 50 with 300 epochs (Or 50 * 300 = 15,000 function evaluations). Testing is performed in 3 phases.
• M,o as identity matrix and null vector respectively. Average taken over 30 iterations. • M as identity matrix and o as a random vector whose value will be changed after each 10th iteration. Average taken over 100 iterations. • M as random rotation matrix and o as a random vector.
Their values will be changed after each 12th iteration. Average and Standard deviation taken over 120 iterations.

Acquiring phase
The acquiring phase is governed by two parameters i.e. r 1 and r 2 . In acquiring phase, if A i is fitter than A s , then A i is updated as else it is updated as Here A s is a candidate solution such that i = s and r 1 , r 2 ≥ 0. We start with Eq. (4.15). The equation can be rewritten as Here the superscript t denotes the value at tth iteration. Replacing as A s is a member of the swarm and s = i, we get In the update equation the dimensions are independent of each other. So it can be assumed as an one dimensional problem without loss of generality. Then the Eq. (4.17) can be rewritten as  where θ = α m (±a) χ As per Von Neumann stability criterion, |e Jβ m τ | ≤ 1. So, Thus we obtain the Eq. (4.20). The inequality applies for all values of sin 2 θ 2 .
But as X s, j − X i, j can be negative, r 1 can be replaced with −r 1 . The region obtained by replacing with −r 1 has boundaries defined by So the complete region for r 1 and r 2 , as shown in Fig. 7a, is equal to the union of the regions defined by Eqs. (4.23a)-(4.23c) and Eqs. (4.24a)-(4.24c). For ease of calculation, this region has been approximated by the rectangular region bounded by the constraint 0 ≤ r 1 ≤ 1 and 0 ≤ r 2 ≤ 1 as shown in Fig. 7b.
The calculations for Eq. (4.14) is similar and by following the same steps which was done for Eq. (4.15), we obtain the equation as As per von neumann stability criterion |e Jβ m τ | ≤ 1. So, Thus we obtain the inequality The inequality applies for all values of sin 2 θ 2 . Putting sin 2 θ 2 = 0 in Eq. (4.26) gives us the same condition as in Eq. (4.22). Putting sin 2 θ 2 = 1 in Eq. (4.20) we get, Thus the resulting inequality is Using the condition that r 1 ≥ 0 with Eq. (4.27), the triangular region with boundaries defined by the lines Eqs. (4.24a)-(4.24c) is obtained. Again, as done before, replacing r 1 with −r 1 , the region bounded by the lines is defined by Eqs. (4.23a)-(4.23c). So for Eq. (4.14) the same rectangular region as Eq. (4.14), which is bounded by the constraint 0 ≤ r 1 ≤ 1 and 0 ≤ r 2 ≤ 1 is obtained.

Experimental results
To verify our findings, we remove the acquiring phase from SGO algorithm and only keep the improving phase. The stable range of r is compared against 3 schemes of unstable ranges of r , i.e.
• M,o as identity matrix and null vector respectively. Average and Standard deviation taken over 30 iterations. • M as identity matrix and o as a random vector whose value will be changed after each 10th iteration. Average and Standard deviation taken over 100 iterations. Their values will be changed after each 12th iteration. Average and Standard deviation taken over 120 iterations. Figure 8 shows the results of the 3 phases mentioned above. We can clearly see that the stable range performs better than the 3 unstable schemes.

Numerical experiments
All the simulations were carried out in Matlab R2016a on the system having Intel Core i7 2.67 GHz processor and 8 GB RAM.
To verify the findings, the performance of SGO algorithm (with both acquiring and improving phase) is tested with the           ranges of r , r 1 , r 2 found in the "Improving phase" section and "Acquiring phase" section. These ranges are referred to as stable range. We use r ∈ [−c − 2, −c − 1] ∪ [c + 1, c + 2] for the improving phase (defined in "Improving phase" section (b)) and r 1 ∈ [1, 1.5], r 2 ∈ [2, 3] for the acquiring phase (defined in "Acquiring phase" section (a)) to test the performance of SGO when the parameters are outside of the stable range. These ranges are referred as Unstable Range. Finally the ranges of SGO mentioned in the original paper [3] is taken to test how well the stable range performs compared to the default ranges. These ranges are referred to as Default Range. The test is performed for 10D and 30D versions of the functions mentioned in Table 1 using N = 10,epoch= 100 and N = 50, epoch= 300 respectively. Values of M and o   Table 2 for the 10D experiment and Figs. 11, 12 and Table 3 for the 30D experiment.
The results show that the default range and stable range perform nearly similar while the unstable range gives the worst results. These results are more apparent as the c value increases.   Table 4 for the 10D experiment and Figs. 15, 16 and Table  5 for the 30D experiment. The unstable range clearly gives  the worst result. Between default range and stable range,the stable range converges a bit slower than default range in some cases, but it always goes on to find a better value than the default range in all the functions. The controlled randomness may be the reason. In most of the algorithms, it is assumed that the range of random parameters should be (0,1), but the relationship and dependencies among the parameter is ignored, which should not be the case.   Table 6 could be observed that stable range performed better than the other unstable ranges for the 10D     Table 7 for the 30D experiment same observations were obtained. Here also unstable range performed worst and stable range always converged.

Conclusion
Stability analysis helps in determining the reliability factor of an algorithm. In this paper, Von Neumann stability analy-sis procedure was used to determine the appropriate range of parameters for which SGO algorithm always converges. The results were supported experimentally with suitable figures and tables. From the analysis, it is deduced that unstable range of parameters of an algorithm may degrade its performance and should be avoided. Moreover, dependencies and relationship between the parameters of an algorithm should be determined because they contribute a lot to the performance and stability analysis provides a way to solve the purpose. This theoretical analysis procedure could be combined with other experimental methodologies in determining the robustness of any algorithm and thus minimizing the chance of failure. Stability analysis of many recent algorithms such as monarch butterfly optimization (MBO) [23], earthworm optimization algorithm (EWA) [24], elephant herding optimization (EHO) [25], moth search (MS) algorithm [26], Slime mould algorithm (SMA) [27], Harris hawks optimization (HHO) [28] and Past Present Future optimization(PPF) [29] which claims to provide promising results experimentally, yet cannot be found in the literature. So, this paves a way for the researchers to work upon more such algorithms.

Conflict of interest The authors have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.