Running-Time Analysis of Brain Storm Optimization Based on Average Gain Model

The brain storm optimization (BSO) algorithm has received increased attention in the field of evolutionary computation. While BSO has been applied in numerous industrial scenarios due to its effectiveness and accessibility, there are few theoretical analysis results about its running time. Running-time analysis can be conducted through the estimation of the upper bounds of the expected first hitting time to evaluate the efficiency of BSO. This study estimates the upper bounds of the expected first hitting time on six single individual BSO variants (BSOs with one individual) based on the average gain model. The theoretical analysis indicates the following results. (1) The time complexity of the six BSO variants is O(n) in equal coefficient linear functions regardless of the presence or absence of the disrupting operator, where n is the number of the dimensions. Moreover, the coefficient of the upper bounds on the expected first hitting time shows that the single individual BSOs with the disrupting operator require fewer iterations to obtain the target solution than the single individual BSOs without the disrupting operator. (2) The upper bounds on the expected first hitting time of single individual BSOs with the standard normally distributed mutation operator are lower than those of BSOs with the uniformly distributed mutation operator. (3) The upper bounds on the expected first hitting time of single individual BSOs with the U−12,12 mutation operator are approximately twice those of BSOs with the U(−1,1) mutation operator. The corresponding numerical results are also consistent with the theoretical analysis results.


Introduction
The swarm intelligence algorithm is one of the nature-inspired optimization algorithms that simulates the behavior of biological groups in nature [1][2][3].Over the past two decades, many different types of swarm intelligence algorithms have been proposed, such as particle swarm optimization (PSO) [4,5], ant colony optimization (ACO) [6,7], artificial bee colony (ABC) [8,9], and brain storm optimization (BSO) [10][11][12].Different from traditional methods, these algorithms solve problems by simulating the behavior of animal or human groups, with higher flexibility and adaptability.
In addition, the theoretical analysis of BSO is also very important, especially for the practical application of BSO.The theoretical analysis benefits researchers in enabling them to understand the mechanism of the algorithm in guiding its design, improvement, and application in practice.The theoretical analysis can be divided into convergence analysis and running-time analysis.Zhou et al. [31], Qiao et al. [32], and Zhang et al. [33] have performed corresponding convergence analyses on BSO.The BSO-CAPSO algorithm, proposed by Alkmini [1], effectively enhances the computational efficiency of BSO through hybridization with chaotic accelerated particle swarm optimization.However, there are few works on the running-time analysis of BSO.Some theoretical methods have been proposed as general analysis tools to investigate the running time of random heuristic algorithms, including the fitness level method [34], drift analysis method [35], switch analysis method [36], wave model [37], etc.These methods are mainly used to analyze discrete random heuristic algorithms.In contrast, fewer theoretical analysis results have been obtained for continuous random heuristic algorithms [38][39][40].However, a large number of practical application problems are continuous.Therefore, the running time of continuous random heuristic algorithms has important research significance.To analyze the running time of continuous random heuristic algorithms, Huang et al. [41] proposed the average gain model.
Huang et al. [41,42] and Zhang et al. [43] used an average gain model to evaluate the expected first hit time of the (1 + 1)-evolutionary algorithm ((1 + 1) EA), evolutionary strategy (ES), and covariance matrix adaptation evolution strategy (CMA-ES).The concept of the first hit time refers to the minimum number of iterations required before the algorithm finds an optimal solution [35].The expected first hit time represents the average number of iterations needed to find the optimal solution, which actually reflects the average time complexity of the algorithm [44].Wang [45] employed the average gain model to analyze the computational efficiency of the proposed swarm intelligence algorithm and provided theoretical evidence for its effectiveness.Therefore, the expected first hit time is a core metric in runtime analysis.Based on the average gain model, the expected first hit time of BSO is deeply analyzed in this paper.
The core of BSO consists of three key components: clustering, interruption, and update.The mutation operator plays an important role in BSO and is included in the interrupt operation and update operation.Specifically, the mutation operator helps the algorithm to jump out of the local optimal solution and further explore a broader search space by introducing random factors in the search process.For example, Zhan et al. [46] proposed an improved BSO (MBSO) in which the mutation operator employs a novel thought difference strategy (IDS).This strategy takes advantage of the thought differences among individuals in the group and increases the diversity of the group by introducing random factors, thus increasing the probability of the algorithm finding the global optimal solution.In addition, El-Abd [47] improved the step equation of the mutation operator and improved the performance of BSO by adjusting the step size and distribution.This improvement helps the algorithm to balance the local search and the global search better, so that the algorithm can find the global optimal solution more effectively when solving complex optimization problems.
The time complexity of the single individual BSO is analyzed in this paper based on the research process from simple to complex.The single individual BSO without the disrupting operation is the same as the (1+1) ES [41].However, the corresponding results can explain the influence of the mutation operator and disrupting operator on the time complexity of BSO.In this paper, we choose the three most classic and representative distributions, N(0, 1), U − 1 2 , 1 2 , and U(−1, 1), as the analysis objects for the mutation operator.Therefore, six BSO variants are obtained as the analyzed algorithms based on the combination of three mutation operators and the presence or absence of a disrupting operator.
The remainder of this paper is organized as follows.Section 2 introduces the process of BSO and the mathematical model of running-time analysis for BSO.Section 3 provides the theoretical analysis results for the running-time analysis of three different BSO variants.Section 4 presents the corresponding experimental results to evaluate the theoretical results.Finally, Section 5 concludes the paper.

Brain Storm Optimization
BSO was proposed by Shi [48,49] in 2011, and it is simple in concept and easy to implement.BSO can be simplified as follows.
In Steps 4 and 5, the new solution is generated by x = x + ∆, where x is the original individual, x is the newly generated individual, and ∆ is a vector generated according to the mutation operator.In this paper, we focus on the running time of BSO with three different mutation operators.The same mutation operators are used to generate new individuals in both Steps 4 and 5.The superposition of different mutation operators is not considered in this work.
To accurately observe the effects of the disrupting and updating operations on the running time of BSO, we select a single individual form of BSO as the analyzed object.The single individual BSO framework simplifies the effect of the population size, which helps to evaluate the effect of the disrupting and updating operations on the running time.Furthermore, following the principle from simple to difficult, the single individual BSO is a suitable starting point for the running-time analysis of BSO.Moreover, the randomness of ∆ and the design of the operations in Steps 4 and 5 are derived from the mutation operator design of evolutionary programming [50,51].Therefore, the conclusion of this analysis will have positive implications for the study of similar mutation operators in evolutionary programming algorithms.

Stochastic Process Model of BSO
BSO can be represented as a stochastic process.In this section, we introduce the terminology for the analysis of the running time of BSO.
Definition 1 (Hill-climbing task).Given a search space S ⊆ R n and a fitness function f : S → R, the hill-climbing task is to find a solution ⃗ x * ∈ S, where the fitness of ⃗ x * reaches the target value H, where f (⃗ x * ) ≥ H.
In this paper, we focus on analyzing the BSO running time in the hill-climbing task of a continuous search space.
Definition 2 (State of BSO).The state of BSO at the t-th (t = 0, 1, . ..) iteration is defined as An optimization problem is a mapping from a decision space to an objective space, and the state space of BSO represents the corresponding decision space.Definition 4 (Expected first hitting time).Let {X t } ∞ t=0 be a stochastic process, where, for any t ≥ 0, X t ≥ 0 holds.Suppose that X t is the Euclidean distance value of the t-th iteration state of BSO to the target solution, and the target threshold ε > 0, the first hitting time [35] of the ε-approximation solution, can be defined by Therefore, the expected first hitting time [44] of BSO can be denoted with E(T ε |X 0 ).
The expected first hit time refers to the average number of iterations required for the BSO algorithm to reach the target fitness value.This metric can more accurately measure the performance of an algorithm because it takes into account the probability distribution of the algorithm over different iterations.Through this metric, we can evaluate the average time complexity of the algorithm in finding the optimal solution, so as to better understand the efficiency of the algorithm.

Running-Time Analysis of BSO Based on Average Gain Model
Inspired by drift analysis [52] and the idea of Riemann integrals [53], Huang et al. [41] proposed the average gain model.Zhang et al. [43] separated this model by introducing the concepts of the supermartingale and stopping time.Based on the former research results [41,43] of the average gain model, Huang et al. [42] proposed an experimental method to estimate the running time of the continuous evolution algorithms.
The expected one-step variation is called the average gain, where {X t } ∞ t=0 is a stochastic process, and H t = σ(X 0 , X 1 , . . ., X t ).The σ-algebra H t contains all the events generated by X 0 , X 1 , . . ., X t .All the information observed from the original population to the t-th iteration is recorded in H t .
Based on Definition 2, is the state of BSO at the t-th iteration.The process of BSO in solving the hill-climbing task is considered as the gradual process of the stochastic state from the initial population to the population that contains the optimal solution.Let f (P t ) = max{ f (⃗ x) : ⃗ x ∈ P t } be the highest fitness value of individuals in P t .X t is used to measure the distance of the current population to the population of the target value.X t = f * − f (P t ), where f * is the fitness value of the optimal solution.Obviously, {X t } ∞ t=0 is a non-negative stochastic process.
The state of BSO in the (t + 1)-th iteration P t+1 only depends on P t .In other words, the stochastic process {P t } ∞ t=0 can be modeled by a Markov chain [42].Similarly, {X t } ∞ t=0 can also be regarded as a Markov chain.In this case, the average gain δ Based on Th. 2 of [43], the expectation of T ε of BSO can be estimated as follows.
Theorem 1. Suppose that {X t } ∞ t=0 is a stochastic process associated with BSO, where X t ≥ 0 for all t ≥ 0. Let h : [0, A] → R + be a monotonically increasing and integrable function. x.
Theorem 1 shows the upper bounds on the expected first hitting time of BSO based on the average gain model.The average gain δ t plays a key role in analyzing the first hitting time T ε of the ε-approximation solution for BSO.The higher average gain indicates a more efficient iteration of the optimization process.

Running-Time Analysis of BSO Instances for Equal Coefficient Linear Functions
In this section, we present the theoretical analysis results based on the average gain model to analyze the expected first hitting time of BSO for equal coefficient linear functions.The running time of BSO with three different mutation operators is analyzed from the perspective of whether the disrupting operation exists.In this paper, we refer to the BSO without a disrupting operation as BSO-I, and the BSO with a disrupting operation as BSO-II.
On this basis, the equal coefficient linear functions are selected as the research object [54][55][56].These functions are a form of basic continuous optimization problem whose function expression is as follows: where (x 1 , x 2 , . . ., x n ) ∈ S. It is assumed that the function starts from the origin and sets the target fitness value to na, where a > 0. The objective of optimizing the equal coefficient linear function is to find a solution ⃗ x * ∈ S, such that f (⃗ x * ) ≥ na.
The mutation operator that obeys the Gaussian distribution and the uniform distribution is selected for the evaluation of BSO.The Gaussian distribution and uniform distribution are common tools for the design of mutation operators [50,51,57], so it is representative to select these two distributions as research cases.

Case Study of BSO without Disrupting Operator
Since the single individual BSO is analyzed in this paper, λ is equal to 1 in the state The BSO of a single individual has only one individual, so ξ t 1 represents both the optimized individual and the random state of the algorithm.The procedure of single individual BSO can be described as follows when the disrupting operation does not exist (i.e., Step 4 in Algorithm 1 is ignored).

Algorithm 1 Brain Storm Optimization (BSO)
1: Initialization: Randomly generate λ individuals (potential solutions) to form the initial population P = {ξ 1 , ξ 2 , . . ., ξ λ } and evaluate the λ individuals; 2: while fail to achieve the predetermined maximum number of iterations do Updating: Randomly choose one or two clusters to create a new individual; Compare the newly generated individual and the original individual with the same individual index.The better one will be saved as the new individual; Update the whole population; the offspring population is recorded as Evaluate the individuals in P ′ ; 6: end while 7: Output the most optimal solution discovered.
is defined as the Euclidean distance of the t-th iteration to the optimal solution.The gain at t-th is given by 3.1.1.When z i ∼ N(0, 1) If the mutation operator obeys the standard normal distribution N(0, 1), the distribution function of η t is as presented by Lemma 1.
Lemma 1.For BSO-I, if its mutation operator ⃗ z obeys N(0, 1), the distribution function F(u) = P(η t ≤ u) of the gain η t is Proof.According to Step 4 and Step 5 of Algorithm 2, the (t + 1)-th individual is Algorithm 2 BSO-I 1: Initialization: Randomly generate an individual ⃗ x = (x 1 , x 2 , . . ., x n ) ∈ R n based on uniform distribution; 2: while stopping criterion is not satisfied do 3: Since z i ∼ N(0, 1), z 1 , . . ., z n are independent of each other.All of the z i satisfy the additivity of the normal distribution, so f (⃗ z t ) obeys the distribution of N(0, nk 2 ).
Hence, the distribution function of η t is shown as F(u) = P(η t ≤ u).
(2) If u = 0, the probability density function of N(0, n) is symmetric in the y axis, so Lemma 1 holds.
Theorem 2 is presented based on the above proof.
Theorem 2. If the mutation operator ⃗ z of BSO-I obeys N(0, 1), the upper bound on the expected first hitting time to reach the target fitness value na is derived as follows.
It is assumed that the algorithm starts from the origin at initialization, where ⃗ x 0 = (0, 0, . . ., 0), i.e., According to Theorem 1, the upper bound on the expected first hitting time is derived as Theorem 2 holds.
Theorem 2 indicates that for BSO-I, if its mutation operator ⃗ z obeys N(0, 1), the com- putational time complexity of BSO-I for the equal coefficient linear function is The uniform distribution function does not satisfy additivity like the normal distribution function.The Lindeberg-Levy center limit theorem [58] can provide an idea to find the distribution of η t .The Lindeberg-Levy center limit theorem is introduced below.
Suppose that X n is a sequence of independent and identically distributed random variables with E(X i ) = µ and then is satisfied for any real number y.The Lindeberg-Levy center limit theorem [58] shows that if n is sufficiently large, Y * n ∼ N(0, 1), it has ∑ n i=1 X i ∼ N(nµ, nσ 2 ).Generally, the case of higher dimensions requires more attention in the study of the computational time complexity of algorithms.If the mutation operator obeys U − 1 2 , 1 2 , the distribution function of η t can be represented by Lemma 2.
Theorem 3 is presented based on the above proof.
Theorem 3.For BSO-I, if its mutation operator ⃗ z obeys U − 1 2 , 1 2 , the upper bound on the expected first hitting time to reach the target fitness value na is derived as follows.
The proof of this theorem is based on the same principle as Theorem 2. The detailed derivation is given as follows.
According to Theorem 1, the upper bound on the expected first hitting time is derived as Theorem 3 indicates that if the mutation operator ⃗ z of BSO-I obeys U − 1 2 , 1 2 , its computational time complexity is E(T ε |X 0 ) = O( √ n) for the equal coefficient linear function.

When z i ∼ U(−1, 1)
If the mutation operator obeys U(−1, 1), the distribution function of η t can be represented by Lemma 3. Lemma 3.For BSO-I, if its mutation operator ⃗ z obeys U(−1, 1), the distribution function F(u) = P(η t ≤ u) of the gain η t is The proof of this lemma is based on the same principle as Lemma 2.
Theorem 4. For BSO-I, if its mutation operator⃗ z obeys U(−1, 1), the upper bound on the expected first hitting time to reach the target fitness value na is derived as The proof of this theorem is based on the same principle as Theorem 2. The proof is given as follows.
According to Theorem 1, the upper bound on the expected first hitting time is derived as Theorem 4 indicates that if the mutation operator ⃗ z of BSO-I obeys U(−1, 1), its computational time complexity for the equal coefficient linear function is The time complexity of BSO-I with three different mutation operators is O( √ n).In the next section, we will discuss the running time of BSO considering the case with a disrupting operation.

Case Study of BSO with Disrupting Operator
Based on the average gain model, this section analyzes the upper bounds of the expected first hit time in three BSO cases.When interference operations are added to the single individual BSO, the algorithm process can be simplified as follows.
In Algorithm 3, the disrupting operations, which are executed with a small probability, are shown in Steps 3 to 6.Let A = {Pb ′ |Pb ′ < Pb} indicate that replacement occurs, while Ā = {Pb ′ |Pb ′ ≥ Pb} indicates that no replacement occurs.
To highlight the effect of each mutation operator on the algorithm, we choose the same mutation operators of BSO-I in Steps 5 and 8 to generate new individuals.As a result, b i = x i + ∆x i , where mutation operator parameter ∆x i and z i follow the same distribution.
Here, ⃗ x t = (x t 1 , x t 2 , . . ., x t n ) ∈ S is still the t-th individual of the algorithm.We have , and the corresponding gain at t is given by η ) If Pb ′ ≥ Pb, it is the same as the result of the case with no disrupting operation in Section 3.1, and the average gain is (2) If Pb ′ < Pb and the mutation operator obeys N(0, 1), the distribution function of η t is represented by Lemma 4.

Algorithm 3 BSO-II
1: Initialization: Randomly generate an individual ⃗ x = (x 1 , x 2 , . . . ,x n ) ∈ R n based on uniform distribution; 2: while stopping criterion is not satisfied do 3: Randomly generate a value Pb ′ from 0 to 1 based on uniform distribution; if ⃗ y has better fitness than ⃗ x then 13: replace ⃗ x with ⃗ y 14: end if 15: end while Output: ⃗ x Lemma 4. For BSO-II, if its mutation operator ⃗ z obeys N(0, 1) and Pb ′ < Pb, the distribution function of the gain η t is F(u) = P(η t ≤ u).
Proof.According to Step 12 and Step 13 of Algorithm 3, the (t According to the definition of η t , t = 0, 1, . .., Since z i ∼ N(0, 1), ∆x i ∼ N(0, 1), and z 1 , . . ., z n are independent of each other, ∆x 1 , ∆x 2 , . . ., ∆x n are also independent of each other.All of z i and ∆x i satisfy the additivity of the normal distribution.As a result, η t obeys N(0, 2nk 2 ).
Hence, the distribution function of η t is F(u) = P(η t ≤ u). ( Lemma 4 holds. Theorem 7 indicates that if the mutation operator of BSO-II ⃗ z obeys U(−1, 1), its computational time complexity for the equal coefficient linear function is

Summary
Overall, we summarize the theoretical analysis results of the running-time analysis of the single individual BSO in solving the n-dimensional equal coefficient linear function in six different situations.The theoretical analysis results are shown in Table 1.The time complexity of the single individual BSO in these six cases is O( √ n).However, the coefficients in the display expressions are different.
Table 1 shows the correlation between the expected first hitting time, the dimension n, the slope k, and the parameter a.The upper bounds on the expected first hitting time of BSO-II are lower than those of BSO-I.This means that the performance of BSO-II is better than that of BSO-I in solving the equal coefficient linear function.The disrupting operation in BSO helps to reduce the running time of the algorithm.Moreover, the upper bounds on the expected first hitting time of the algorithm using the standard normal distribution mutation operator are lower than those of the algorithms with the uniform distribution mutation operator.In addition, the upper bounds on the expected first hitting time of the algorithm using the U − 1 2 , 1 2 mutation operator are approximately two times higher than those of the algorithm using the U(−1, 1) mutation operator.

Experimental Results
In Section 3, we obtain the theoretical analysis results of the expected first hitting time of single individual BSOs through the average gain model.To verify the correctness of the analysis results, numerical experiments are presented in this section.
As the number of samples increases, the arithmetic mean will gradually approach true mathematical expectations based on the Wiener-Khinchine theorem of large numbers [58].The Wiener-Khinchine theorem of large numbers is introduced as follows.
Suppose that X 1 , X 2 , . . . is a sequence of independent and identically distributed random variables with E(X i ) = µ; for any ε > 0, the following equation will hold.
The Wiener-Khinchine theorem of large numbers indicates that if the number of samples is sufficiently large, the mathematical expectations are approximately equal to the mean of samples X 1 , X 2 , . . ., X n .Therefore, we use the arithmetic mean of the first hitting time of multiple experiments to estimate the actual expected first hitting time.
In this section, the parameters of the proposed approach are set as follows.The fixing error is ε = 1 × 10 −8 , the initial individual is x 0 = (x 0 1 , x 0 2 , . . ., x 0 n ) = (0, 0, . . ., 0), the slope is k = 1, and the target fitness parameter is a = 10.The problem dimension n is set from 10 to 280.BSO-I and BSO-II are conducted on the n-dimensional equal coefficient linear function for 300 runs for each dimension.The termination criterion for each experiment is that the error of the optimal solution should be below a predefined threshold ε.Table 2 shows the numerical results of the practical expected first hitting time E(T ε |X 0 ) and the theoretical time upper bound, where E(T ε |X 0 ) = ∑ 300 i=1 T εi 300 , and T εi is the first hitting time of the ε-approximation solution at the i-th run.The points larger than the theoretical upper bounds are highlighted in boldface.
As shown in Table 2, the experimental results strongly fit the theoretical upper bounds, indicating that the error between the theoretical upper bounds and the actual value is within ε.The points larger than the theoretical upper bounds are highlighted in boldface.The arithmetic mean of multiple experiments is used to estimate the expected first hitting time.In the real case, only the arithmetic mean of 300 experiments is used to estimate the expected first hitting time, which allows a certain statistical error.According to the central limit theorem, the results obtained from 300 independent experiments follow a normal distribution.The null hypothesis H 0 means that the mean value of the expected first hitting time in 300 experiments is less than or equal to the corresponding theoretical upper bounds.The corresponding significance level α is 0.05 with the T testing.Moreover, as shown in Figure 1, the actual expected running time is followed with the estimated result based on our proof.All the detailed results are shown in Table 3.
Table 3 provides the numerical results, where h represents the hypothetical result, p represents the p-value of the test, and ci is the confidence interval.As shown in the T testing of Table 3, h = 0 and p > α.The null hypothesis H 0 is accepted at the significance level α = 0.05.Therefore, the analytic expression of the running time of BSO obtained based on the average gain model can characterize the actual upper bounds of the running time of BSO in these six BSO variants.

Conclusions
The running time of six BSO variants for the equal coefficient linear function is analyzed in this paper based on the average gain model.The additivity of the normal distribution and the Lindeberg-Levy center limit theorem are applied to deal with the superposition of the normal distribution mutation operator and the uniform distribution mutation operator, respectively.Furthermore, the full expectation formula is utilized to deal with the problem of individual replacement with a certain probability in the disrupting operator.This paper also concludes the upper bounds on the expected first hitting time of the single individual BSO in equal coefficient linear functions.
The analysis results show that the time complexity of BSO-I and BSO-II is O( √ n) in the equal coefficient linear function.However, their coefficients are different.In the linear function with equal coefficients, the upper bound of the expected first hit time of BSO-II is smaller than that of BSO-I.In addition, the single individual BSO using the standard normally distributed mutation operator expects a lower upper bound on the first hit time than the corresponding algorithm using the uniformly distributed mutation operator.The upper bounds on the expected first hitting time of single individual BSOs with the U(− 1 2 , 1 2 ) mutation operator are approximately twice those of BSOs with the U(−1, 1) mutation operator.
In our future work, we will analyze the running time of the population-based BSO in the equal coefficient linear function.The running time of population-based BSOs in practical optimization problems is also an important topic.Moreover, it is crucial to extend our research to practical optimization problems that are encountered in real-world applications with complex constraints, non-linear relationships, or high-dimensional spaces.

3 :
Clustering: Use clustering algorithms to divide λ individuals into m clusters; 4:Disrupting: The mutation occurs with a certain probability, and a randomly selected cluster's central individual is replaced by a randomly generated new individual; 5:

4 :if
if Pb ′ is smaller than a pre-determined probability Pb then 5: replace ⃗ x with a randomly generated individual ⃗ b = (b 1 , b 2 , . . ., b n ) based on uniform distribution; Pb ′ < Pb then 8:⃗ y = ⃗ b +⃗ z, where ⃗ z is a mutation operator;

Figure 1 .
Figure 1.The curve of the expected and the actual running time.The six figures, arranged from top to bottom and left to right, depict the three distributions corresponding to BSO-I and the three distributions corresponding to BSO-II, respectively.
where ⃗ z is the mutation operator;According to the definition of η t , where t = 0, 1, . .., (1) If f

Table 1 .
Analysis of the running time of BSO in six different situations.

Table 2 .
Comparison of the estimation of the expected first hitting time and the theoretical upper bounds.

Table 3 .
Statistical results of hypothesis testing.