https://doi.org/10.20965/jaciii.2021.p0177
Self-Organized Subpopulation in GP on GPU

Paper:

Self-Organized Subpopulation Based on Multiple Featuresin Genetic Programming on GPU
Keiko Ono∗ and Yoshiko Hanada∗∗∗ Doshisha University

1-3 Tatara Miyakodani, Kyotanabe, Kyoto 610-0394, Japan
E-mail: kono@mail.doshisha.ac.jp∗∗ Kansai University
3-3-35 Yamate-cho, Suita, Osaka 564-8680, Japan[Received August 9, 2019; accepted December 3, 2020]

Genetic Programming (GP) is an Evolutionary Computation (EC) algorithm. Controlling genetic diversity
in GP is a fundamental requirement to obtain various types of local minima effectively; however, thiscontrol is difficult compared to other EC algorithms
because of difficulties in measuring the similarity between solutions. In general, common subtrees and theedit distance between solutions is used to evaluate thesimilarity between solutions. However, there are noclear guidelines regarding the best features to evaluate
it. We hypothesized that the combination of multiple
features helps to express the specific genetic similarity of each solution. In this study, we propose a selforganized subpopulation model based on similarity interms of multiple features. To reconstruct subpopulations, we introduce a novel weighted network basedon each normalized feature and utilize network clustering techniques. Although we can regard similarity as a correlation network between solutions, the useof multiple features incurs high computational costs,
however, calculating the similarity is very suitable forparallelization on GPUs. Therefore, in the proposed
method, we use CUDA to reconstruct subpopulations.
Using three benchmark problems widely adopted instudies in the literature, we demonstrate that performance improvement can be achieved by reconstructing subpopulations based on a correlation network ofsolutions, and that the proposed method significantly
outperforms typical methods.

Keywords: genetic programming, subpopulation model,
genetic diversity, multiple features

1. Introduction
In Genetic Programming (GP), it is important to maintain genetic diversity, which can be achieved by improving the replacement strategies, such as re-selection andcrowding [1], and the sharing function [2]. A subpopulation model easily encourages genetic diversity and avoids
premature convergence toward local optima; it frequently
Vol.25 No.2, 2021

shows a better search performance than single population
models. Studies have been conducted on subpopulation
models in GP to further improve its performance. For
example, Hu et al. [3, 4] proposed the hierarchical faircompetition model (HFX), in which individuals are separated into subpopulations according to their fitness values
and evolve in a fair competition environment. Hornby [5,6] proposed the age-layered population structure model(ALPS), in which individuals are evaluated by their age,that is, a measure of how long the genetic material has
been in the population, and separated into subpopulations
according to this age measure. Luna et al. [7] proposeda meaning-based model for data mining to extract mining
rules. Moreover, other researchers have proposed migration models [8–11], in which some individuals migrate to
the next subpopulation to maintain genetic diversity. HFXand ALPS require appropriate boundary parameters foreach subpopulation, which are difficult to determine: a fitness value in HFX and an age value in ALPS. By contrast,in migration models, it is not necessary to set boundary
parameters, and these models can be easily implemented
in various applications. Therefore, we focus on migration
models in this study and consider enhancing the subpopulation models in the GP framework.
To enhance the subpopulation model, it is important toconsider the amount of genetic diversity that should bemaintained during evolution. The higher the genetic diversity of a subpopulation, the wider the range of characteristics of each solution, which hinders the convergenceto local minima. Previous subpopulation models focusedon increasing genetic diversity by dividing solutions intosome populations; the level of genetic diversity was notcontrolled. The genetic diversity can be easily controlled
if the distance between solutions can be measured as Euclidian distance. However, this distance is difficult tomeasure because solutions in GP are typically expressedby a tree. In general, the similarity based on common
subtrees and edit distance between solutions is used toevaluate the distance between solutions. Nevertheless, thebest manner to express distance is uncertain because eachproblem domain has its own most effective feature. We
hypothesize that multiple features extracted from each solution can enable an improved evaluation of the distance

Journal of Advanced Computational Intelligenceand Intelligent Informatics© Fuji Technology Press Ltd. Creative Commons CC BY-ND: This is an Open Access article distributed under the terms of
the Creative Commons Attribution-NoDerivatives 4.0 International License (http://creativecommons.org/licenses/by-nd/4.0/).

177

Ono, K. and Hanada, Y.

between solutions compared to a single feature, and adaptive reconstruction of subpopulations consisting of similar
individuals in terms of multiple features can enhance local
search. Moreover, both mutations in each subpopulationand independent evolution using the subpopulation model
can enhance genetic diversity.
In this study, we propose a novel self-adaptive subpopulation model that, on the basis of multiple features, creates a subpopulation in which similar individuals exist andenhances mutation in this subpopulation. More specifically, to create subpopulations with similar features, we
introduce a weighted network of individuals based on
similarity obtained from features, and separate individuals into subpopulations using network clustering. Using
these strategies, we can control the balance between local
search and genetic diversity to avoid redundant searching. Obtaining similarity based on multiple features is atime-consuming task, whereas calculating the similarity
between individuals is very suitable for parallelization.
Therefore, we use CUDA to construct a correlation network. It is well known that a considerable computational
effort is required to evaluate fitness in GP. The first study
using CUDA (GPUs) in GP was reported in [12]; otherresearchers have also proposed using CUDA models in
GP to enhance acceleration [13–16]. To the best of ourknowledge, our proposed method is the first subpopulation model in GP, executed on CUDA, which reconstructs
subpopulations according to the individuals’ features.
The remindar of this paper is organized as follows. In
Section 2, we formulate the optimization problem discussed in this study and provide a detailed explanationof the subpopulation model in GP. In Section 3, we describe the CUDA-based proposed adaptive subpopulation
model in GP in detail. We first present an overview ofthe proposed method, and then explain the CUDA-based
approach of the proposed method with the objective of
reducing computational effort. Moreover, we explaina technique for adaptively reconstructing subpopulations
using the proposed method. In Section 4, we examine
the effectiveness of the proposed method through experiments using three benchmark problems widely adopted instudies in the literature. We first explain the implementation of the proposed method, and then introduce two comparison methods to analyze the effects of the proposed
method.1

2. Problem Definition
Model in GP

andapply subpopulation model methods. In the GP subpopulation model, the total population V is partitioned into
M subpopulations, and the subpopulations execute GP
searches for minimizing f (x) in parallel, although theyexchange information by migration. There are several migration topologies such as ring, random, and grid, but thering topology is the most basic. Thus, we focus on thering topology and refer to it as the subpopulation modelin GP in this paper. It should be noted that the number ofindividuals in each subpopulation does not change during
searching.
We begin by recalling the details of the subpopulation
model in GP with ring topology. Let Pm (m = 1, . . . , M)
denote the m-th subpopulation, where the number of individuals in Pm is |Pm | = |V |/M. The individuals in Pm aredenoted by V m = {vmi ; i = 1, . . . , I}, and the total population is obtained by
V=

M


V m,

m=1

where I = |V |/M, and |V | is the number of all individuals. The total population is V = {vi ; i = 1, . . . , I × m}.
It should be noted that the set of subpopulations P ={Pm ; m = 1, . . . , M}, and the maximum generation is G.
The algorithm of the subpopulation model in GP is shownas Algorithm 1. Here, the topology type of the subpopulation model in GP is a ring; therefore, the subpopulation P0 receives individuals from PM , and PM sends themto P0 .

Subpopulation

In the framework of GP, we consider the problem ofminimizing a function f (x), where its input variable x isrepresented as a labeled ordered tree. Given the total population size N and a crossover-mutation-selection strategy, a common approach to improving a GP search is to
1. In our conference paper [17], we proposed a basic framework for sucha migration strategy. This paper extends our preliminary study [17] by
improving the technique of measuring the similarity between subpopulations incorporating mutation.

178

Algorithm 1 Algorithm for subpopulation model
R0: Set g ← 1.
R1: Initialize individuals in each subpopulation Pm .
R2: If g ≤ G, then perform Steps R3 to R8; otherwise,stop.
R3: Apply genetic operations, such as selection andcrossover.
R4: If g is a migration generation, then perform Steps R5
to R7; otherwise, perform Step R8.
R5: Choose the best n individuals Bm and the worst n individuals W m in Pm .
R6: Send Bm to the next subpopulation Pm+1 .
R7: Replace W m with Bm−1 , that is, the best individuals
in the next subpopulation Pm−1 .
R8: Set g ← g + 1.

3. Self-Organized Subpopulation Model in GP
on GPU: SoS-GP
3.1. Overview of SoS-GP
In sequence, we consider enhancing the subpopulation
model in GP for the minimization problem of the objective function f (x).
In GP, the structure of an individual can be expressedby a tree. There are multiple features of an individual,

Journal of Advanced Computational Intelligenceand Intelligent Informatics

Vol.25 No.2, 2021

Self-Organized Subpopulation in GP on GPU

als

idu

div

Ini

3

In

id

ivi

du

als

Fig. 1. Proposed model using CUDA.

such as fitness, tree size, and subtree. However, fitness isused only as a feature for migration in the common subpopulation model in GP; elite individuals with lower fitness values than others in the same subpopulation migrate
to other subpopulations. According to the building block
hypothesis, partial solutions are assembled into the entire solution. Here, we hypothesize that incorporating thestructure information of individuals in the migration strategy is important for enhancing the combination of partial
solutions. We also assume that the balance between local
and genetic searches is important. If abrupt changes occur in the search area, individuals randomly evolve, whichleads to disruption of genetic operations. To effectively
explore partial solutions, we incorporate a local search
into the subpopulation model in GP. Thus, from the viewpoint of encouraging local search to combine partial solutions, we suppose that automatic reconstruction of subpopulations, which consist of similar individuals in termsof their multiple features, enhances local search and improves the combination of partial solutions. Therefore,
we propose a self-adaptive subpopulation model in GP.
In the proposed method, we generate a weighted networkof individuals based on fitness and node size and create
subpopulations by network clustering.
However, fitness evaluation is the most timeconsuming part of reconstructing subpopulations. Moreover, the proposed method evaluates the similarity between individuals based on their features, which also increases the computational effort. To solve this issue, we
use GPGPU to create a similarity network using multiple features. Fig. 1 shows the model of the proposed
method. The proposed method forms subpopulations inwhich individuals are evolved in parallel. In spite of migration, subpopulations are automatically reconstructedby the similarity network of individuals V . In Section 3.2,
we describe the proposed model using GPGPU in detail.
The algorithm of the proposed method is shown as Algorithm 2.

3.2. Feature Evaluation Using GPGPU
In the proposed method, self-organized subpopulation
model in GP on GPU (SoS-GP), an automatic reconstrucVol.25 No.2, 2021

Algorithm 2 Algorithm for the proposed method
P0: Set g ← 1.
P1: Initialize individuals in each subpopulation Pm , where
the number of individuals I is set to |V |/M.
P2: If g ≤ G, then perform Steps P3 to P13; otherwise,stop.
P3: Apply selection and crossover.
P4: Apply mutation γ times.
P5: If g is a reconstruction generation, then perform
Steps P5 to P11; otherwise, perform Step P13.
P6: Send V m in each Pm to the subpopulation P0 .
P7: If the subpopulation is P0 , then perform Steps P7
to P11; otherwise, perform Step P12.
P8: Evaluate the similarity matrix H between vi and v j(i = j).2
P9: Create a network based on H.
P10: Divide the population V into M subpopulations V̂ mby network clustering to reconstruct V m .3
P11: If |V m | is not 1, then perform Steps P11 to P12; otherwise, perform Step P13.
P12: Send V̂ m to Pm .
P13: Replace V m with V̂ m .
P14: Set g ← g + 1.
tion of subpopulations that consist of similar individuals
in terms of their features, is conducted to enhance local
search. In GP, as an individual can be represented by atree, there are several features that can be used to measure their similarity. We assume that the similarity canbe accurately measured using multiple features. However, there are |V |2 /2 combinations for one type of featurewhen the number of individuals is |V |, and the required
computational effort is thus increased. To decrease thecomputational effort when multiple features are applied,
we use GPGPU to create a similarity network.
Figure 1 shows the model of the proposed methodthat uses both CUDA and MPI. The proposed method
evolves subpopulations in parallel and applies the same
type of genetic operators in all subpopulations. MPI is
2. Here, the similarity H is a V ×V matrix, and the element of H is huv . For
more detail, see Section 3.2.
3. For more detail, see Section 3.3.

Journal of Advanced Computational Intelligenceand Intelligent Informatics

179

Ono, K. and Hanada, Y.

mostly known for message-passing multiprocessing programming and allows a subpopulation model to be easily
implemented using multiple processors on a single systemor a cluster of systems. However, the proposed method
evaluates the similarity between individuals in terms ofeach feature in a different CUDA kernel in parallel toachieve acceleration. Therefore, the proposed method
uses both CUDA and MPI. The three steps for evaluating
the similarity are as follows:
1. Gather individuals in each subpopulation to subpopulation 1; the information of individuals is sent to a
global memory;
2. Evaluate the similarity in terms of each feature using
a different kernel: node size or fitness. Each similarity result is sent to subpopulation 1;
3. Combine the similarity results using cuBLAST.
Here, to combine the similarities evaluated by different kernels, we use cuBLAST, which is the CUDA basic
linear algebra subroutine library on NVIDIA CUDA runtime. It automatically allocates the required matrices andvectors in the GPU memory space. cuBLAST has many
calculators of matrix vectors divided into three levels interms of parallelization. In the proposed method, eachsimilarity consists of |V |2 /2 elements; therefore, combining these elements requires a matrix and matrix multiplication operation. This operation is a BLAST Level 3 operation, which has higher scalability than other levels.
Therefore, the proposed method uses the cublasGetMatrix
operator in cuBLAST to process multiple feature similarities.

Algorithm 3 Algorithm for Newman clustering
S1: Evaluate the initial values of ΔQi j and ai according toEqs. (1) and (2).
S2: Select the largest element of each row of the matrix ΔQ, and save it to the max-heap.
S3: Select the largest ΔQi j from the max-heap, and join
the corresponding communities.
S4: Update ΔQ, the max-heap and ai according to Eq. (3),and increment Q by ΔQi j .
S5: Repeat Steps 2 and 3 until only one community remains.
S6: Extract communities where Q is the maximum value.
the similarity matrix H. Thus4⎧⎪⎨huv if the weight between vertices uAuv =and v is not 0,⎪⎩0otherwise.
ai is obtained by
a i = ∑ ei j ,j

ei j =

1∑ Auv δ (cu, i)δ (cv, j),
2W ∑
u v

where W is the amount of weight and δ function δ (i, j)
is 1 if i = j and 0 otherwise. If community i belongs to
the community cu , then δ (cu , i) = 1. As the first step, leteach vertex be the sole member of a community of one;then, ei j = 1/2W if i and j are connected and 0 otherwise,
and ai = ki /2W , whereku = ∑ Auv .
v

3.3. Newman Clustering
Newman clustering [18, 19] is the most well-known
algorithm for extracting communities. It uses a quality function known as “modularity Q” over the possible divisions of a network. By optimizing the modularity, density-connected groups of vertices with only sparseconnections between groups are detected. As in the proposed method we reconstruct subpopulations based on thesimilarity between individuals, Newman clustering canbe applied to detect subpopulations in which individuals
have a similar characteristic easily if we refer the similarity matrix H to an adjacency matrix in a weighted network. In this study, we used the Newman clustering algorithm to reconstruct subpopulations. We refer to an individual as a vertex, and an element of similarity H is theweight between vertices. The modularity Q is defined as
Q = ∑(eii − a2i ),iwhere eii is a portion of the weight of the joint vertices incommunity i in the total weight, and ai is a portion of theweight of nodes that connect vertices in community i toother groups in the total weight. Let Auv be an element of
180

Thus, the initial ΔQi j is obtained by⎧ki k j1⎪⎨−
if the weight between
2W
(2W
)2
ΔQi j =
vertices i and j is not 0,⎪⎩0otherwise,. . . . . . . . . . . . . . . . . . . . (1)andki. . . . . . . . . . . . . . . . (2)2W
The Newman clustering algorithm with the similarity matrix H is shown as Algorithm 3.
⎧
ΔQik + ΔQ jk if community k is connected
⎪
⎪
⎪
⎪
⎪
to both i and j,⎪⎪⎨ΔQik − 2a j ak if community k is connected
ΔQ jk =⎪
to i but not to j,
⎪
⎪
⎪
⎪
−
2a
a
ΔQ
⎪
j k if community k is connectedjk⎪⎩
to j but not to i.
. . . . . . . . . . . . . . . . . . . . (3)ai =4. Auv is an element of the adjacency matrix of the network.

Journal of Advanced Computational Intelligenceand Intelligent Informatics

Vol.25 No.2, 2021

Self-Organized Subpopulation in GP on GPU

4. Experimental Evaluation
Through experiments using three benchmark problems
widely adopted in studies in the literature, we evaluatedthe effectiveness of the proposed method. Moreover, weanalyzed the effects of mutation and genetic diversity, andverified that the proposed framework based on similarity
between individuals can control genetic diversity.

4.1. Implementation
Here, we explain the implementation of Step P7 of the
proposed method in the experiments. Although there areseveral features for measuring the similarity between individuals i and j, we used the difference in the fitness value
and in the node size between i and j for simplicity. The
similarity matrix H f in terms of the difference in the fitness value is obtained by
hifj f (vi ) = f (v j ),f
Hi j =
0otherwise,
where hifj = η /| f (vi ) − f (v j )| and η is a real-valued paf
rameter. We set it to 1.0 in our experiments. Hi j becomeslarge when f (vi ) and f (v j ) are scalar. Moreover, the similarity matrix H n in terms of the difference in node size isobtained byhnij |vi | = |v j |,Hinj =
0otherwise,

where hnij = η /||vi | − |v j ||. We constructed the similarity
of individuals H as follows. Let Hˆ f and Hˆn be the normalized matrices of H f and H n , respectively. Let H denotethe total similarity with respect to the fitness values andnode size. Then, we define H by
H = α Hˆ f + (1 − α )Hˆn .5 . . . . . . . . . (4)A matrix and matrix multiplication operation were used
to normalize the matrices H f and H n and to obtain H in
Eq. (4); therefore, we used the cublasGetMatrix operationin cuBLAST.

4.2. Comparison Methods
The proposed method evaluates the similarity between
individuals using a combination of H f and H n . H f measures the difference between individuals in fitness values,and H n measures the difference in node size. To analyze the effect of H f and H n separately, we introduce two
methods as extreme cases of the proposed method andcompare them with the proposed method, as describedin Section 4.5. First, we investigate the effectiveness ofsimilarity in terms of fitness values. The similarity H isdefined by
H = Hˆ f .
5. α is a parameter. In our experiments, we adopted α = 0.5.

Vol.25 No.2, 2021

We refer to this method as the F-method. Then, we investigate the effectiveness of the similarity in terms of nodesize. Similarity H is then defined by
H = Hˆn .
We refer to this method as the N-method.

4.3. Benchmark Problems
We evaluated the effectiveness of the proposed method
on well-known benchmark problems, the symbolic regression problems. Here, we briefly explain these problems (see [20–22] for detailed descriptions).
We investigated the symbolic regression problem forthe function space X constructed by the labeled ordered
trees of functional nodes {+, −, ×, /, sin, cos, log} andterminal nodes {s, 0.05, 0.10, 0.15, 0.20, . . ., 0.95, 1.00},
where s denotes a variable. Our training set was composed of 30 data points {(s j , x∗ (s j )) ∈ R2 ; j = 1, . . . , 30},
where s j = 0.2( j − 1), and x∗ (s) ∈ X is the true function to be identified. For any x(s) ∈ X , we define thefitness f (x) by
30

f (x) = 50 ∑ |x(s j ) − x∗ (s j )|,j=1

and consider the maximization problem of f (x). In our
experiments, we employed three functions x∗ (s),Function A :x∗ (s) = (2 − 0.3s) sin(2s) cos(3s) + 0.11s2,
Function B :x∗ (s) = s cos(s) sin(s) sin2 (s) cos(s) − 1 ,Function C :x∗ (s) = s3 cos(s) sin(s)e−s sin2 (s) cos(s) − 1 .

4.4. Experimental Setting
To evaluate the benchmark problems, we used the following standard parameter settings [20–22]: the recombination rate was 1.0, the mutation rate was 0.0, random
selection was used (non-elitist), the total population size
was N = 250, the number of subpopulations was M = 5,
the maximum depth was 16, the number of generations
was I = 100, the migration interval was 10, and the initialindividuals were created by using “ramped half-and-half”
with maximum depth. Fig. 2 shows the history of fitness
in terms of Ring and the proposed method. From this figure, we observe that the Ring method converged aroundthe 50th generation; therefore, we set the number of generations to I = 100.
All our experiments were conducted on a single PCwith 6 Intel Xeon E5-1660 3.3 GHz processors, with
24 GB of memory and NVIDIA Quadro K6000 runningon Linux. The NVIDIA Quadro K6000 is capable of approximately 5.2 TFlops and 2880 CUDA cores.
4.5. Performance Evaluation
We first compared the proposed method with the N,
F, Ring, and Simple methods in terms of solution qual-

Journal of Advanced Computational Intelligenceand Intelligent Informatics

181

Ono, K. and Hanada, Y.

600

800

700

Fit ness

650

600

650

600

550

550

500

500

450

0

20

40
60
Generat ion(a) Ring

80100Island0Island1Island2Island3Island4

750

450

550

Fitnes700

Fit ness

800

Island0Island1Island2Island3Island4

750

500

450
0

20

40
60
Generat ion

80

100

400

(b) Proposed

ProposedN-method

F-method

Ring

SimpleRing

SimpleRing

SimpleMethods

Fig. 2. History of fitness (Function A).

(a) Function A
1000
950
900

182

Fitnes

850
800
750
700
650
600
550
500
ProposedN-method

F-method

Methods(b) Function B
400

350

Fitnesity. In the experiments, we conducted 20 trials and evaluated the best individuals xmf in Pm at the final generation, that is, f (xmf ) ≥ f (x) (∀x ∈ V m ). Fig. 3 graphicallypresents f (xmf ), where m = 1, . . . , M and M = 5.
Figure 3 shows that the Simple method was not superior to the other methods in any functions. This result implies the necessity of migration in a subpopulation model.
Fig. 3(a) shows the results of Function A. The N-method
was not superior to the proposed, F and Ring methods
in performance, and the F-method was comparable to the
Ring method. Moreover, the performance of the proposed
method was superior to that of these methods.
Figure 3(b) shows the results for Function B. It is observed that the performance of the Ring method compared to that of the F-method was superior for Function B,
and similar for Functions A and C. By contrast, the performance of the F-method was superior to that of the
N-method. These results imply that similarity based on
fitness is effective for solution quality. In Section 4.6, weanalyze the reason why the F-method was superior to the
N-method.
From Fig. 3(c), it is observed that the performanceof the proposed method was comparable to that of theF-method for Function C, whereas for other functions
the proposed method outperformed the N and F methods.
These results imply that the proposed method that reconstructs subpopulations based on multiple features of individuals, that is, incorporating both the similarity of node
size Ĥ n and the similarity of fitness Ĥ f , leads to performance improvement.
Next, we evaluated the effects of mutation in Step P4.
We hypothesized that the proposed SoS model can enhance local search by constructing a subpopulation; it cansimultaneously maintain genetic diversity through mutation. Here, individuals in the same subpopulation haverelatively the same features. We evaluated the fitness ofeach individual every five generations while changing thenumber of mutations γ . Figs. 4–6 show the distributionsof individuals every five generations using the proposed
method. We found that fitness became lower than that of
the 5th generation as the number of generations increased.
In particular, the results showed better performance in all
problems when γ was set to 10. From these results, we
confirmed that the balance in terms of mutation rate isimportant for performance improvement.

300

250

200

150

100
ProposedN-method

F-method

Methods(c) Function C

Fig. 3. Performance evaluation.

To verify why the proposed method showed better
performance when γ was set to 10, we analyzed thedifferences between the elite and the other individuals.
Figs. 7–9 show the differences between the best individual and each individual in terms of fitness and edit distance, in which individuals evolve without mutation when
γ is set to 0. Compared to γ = 5, 10, and 15, the plots in
Figs. 4(a), 5(a), and 6(a) were gathered around the upperleft corner. These results without mutation show less diversity, as expected from the high similarities in each subpopulation. By contrast, diversity increases as the numberof mutations increases. From these results, we verify thatthe proposed method can control local search and genetic
diversity by the proposed SoS and mutation, and an appropriate balance exists. These results demonstrate theeffectiveness of the proposed method.
Moreover, we evaluated the average of nodes in termsof the best individuals xmf in Pm at the final generation toinvestigate the influence of the bloat phenomenon in GP.

Journal of Advanced Computational Intelligenceand Intelligent Informatics

Vol.25 No.2, 2021

Self-Organized Subpopulation in GP on GPU

1600

Fitness

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

1400
1200
1000
800

600

600

600

600

400

400

400

400

200

200

200

0

0
50

100

150

200

250

200

0
50

100

150

200

250

gene-5gene-10gene-50
gene-95

1800

Fitness

1600

Fitness

2000

gene-5gene-10gene-50
gene-95

1800

Fitness

2000

0
50

100

150

200

250

50

100

150

200

#Indivisuals#Indivisuals#Indivisuals#Indivisuals(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

250

Fig. 4. Distribution of individuals (Function A).

1600

Fitness

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

1400
1200
1000
800

600

600

600

600

400

400

400

400

200

200

200

0

0
50

100

150

200

250

200

0
50

100

150

200

250

gene-5gene-10gene-50
gene-95

1800

Fitness

1600

Fitness

2000

gene-5gene-10gene-50
gene-95

1800

Fitness

2000

0
50

100

150

200

250

50

100

150

200

#Indivisuals#Indivisuals#Indivisuals#Indivisuals(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

250

Fig. 5. Distribution of individuals (Function B).

1600

Fitness

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

2000

gene-5gene-10gene-50
gene-95

1800

1600

1400
1200
1000
800

1400
1200
1000
800

600

600

600

600

400

400

400

400

200

200

200

0

0
50

100

150

200

250

200

0
50

100

150

200

250

gene-5gene-10gene-50
gene-95

1800

Fitness

1600

Fitness

2000

gene-5gene-10gene-50
gene-95

1800

Fitness

2000

0
50

100

150

200

250

50

100

150

#Indivisuals#Indivisuals#Indivisuals#Indivisuals(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

200

250

Fig. 6. Distribution of individuals (Function C).

A common method to prevent bloat is to limit the maximum number of nodes; therefore, we introduced a new
function shown in Eq. (5).
30g(x) = 50 ∑ |x(s j ) − x∗ (s j )| + β · ||x||,. . . (5)j=1where || || indicates the number of nodes in x and β isa weight parameter that controls the number of nodes.
Fig. 10 shows the fitness values in terms of fitness g(x)and f (x) = g(x) − β · ||x|| with respect to β using the proposed method with γ = 0. From Figs. 3 and 10, Functions A and B with β = 1 showed better performance
than the proposed method. As the first step, the proposed
method does not adopt an approach to reduce bloat; however, these results indicate that we should consider the
bloat in GP in the proposed method for further performance improvement by incorporating a bloat resistancemethod.

Vol.25 No.2, 2021

4.6. Behavior Analysis
In the proposed method, multiple features are used tocreate a correlation network of individuals. Therefore,
considerable computational effort is required by the proposed method compared to the Ring and Simple methods
for evaluating the similarity between individuals. First,
we investigated the scalability of the proposed method interms of the number of individuals and used this method
over 20 trials to analyze the best and worst trials based on
the average execution time, as shown in Fig. 11.
We set the number of individuals as 50, 150, 250, 350,
and 450. The average execution time is different in eachproblem because although the maximum depth of an individual is limited, the number of nodes of individuals
is not constant. In Fig. 11, we see that the execution
time increased linearly as the number of individuals increased from 50 to 150. However, linear computational
effort with respect to the number of individuals was not
required when the number of individuals exceeded 150.

Journal of Advanced Computational Intelligenceand Intelligent Informatics

183

-1000

-500

-1000

-1500

-2000
0

50

100

200

0

50

100

200

-500

-1500

gene-5gene-10gene-50
gene-95

-2000

150

0

-1000

-1500

gene-5gene-10gene-50
gene-95

-2000

150

-500

-1000

-1500

gene-5gene-10gene-50
gene-950Difference (Fitness)

-500

0

Difference (Fitness)

0

Difference (Fitness)

Difference (Fitness)Ono, K. and Hanada, Y.

0

50

100

gene-5gene-10gene-50
gene-95

-2000

150

200

0

50

100

150

Edit distance

Edit distance

Edit distance

Edit distance(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

200

-1000

-500

-1000

-1500

-2000
0

50

100

200

0

50

100

200

-500

-1500

gene-5gene-10gene-50
gene-95

-2000

150

0

-1000

-1500

gene-5gene-10gene-50
gene-95

-2000

150

-500

-1000

-1500

gene-5gene-10gene-50
gene-950Difference (Fitness)

-500

0

Difference (Fitness)

0

Difference (Fitness)

Difference (Fitness)

Fig. 7. Differences between the elite and the others (Function A).

0

50

100

gene-5gene-10gene-50
gene-95

-2000

150

200

0

50

100

150

Edit distance

Edit distance

Edit distance

Edit distance(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

200

-1000

-500

-1000

-1500

-2000
0

50

100

200

0

50

100

200

-500

-1500

gene-5gene-10gene-50
gene-95

-2000

150

0

-1000

-1500

gene-5gene-10gene-50
gene-95

-2000

150

-500

-1000

-1500

gene-5gene-10gene-50
gene-950Difference (Fitness)

-500

0

Difference (Fitness)

0

Difference (Fitness)

Difference (Fitness)

Fig. 8. Differences between the elite and the others (Function B).

0

50

100

gene-5gene-10gene-50
gene-95

-2000

150

200

0

50

100

150

Edit distance

Edit distance

Edit distance

Edit distance(a) γ = 0(b) γ = 5(c) γ = 10(d) γ = 15

200

Fig. 9. Differences between the elite and the others (Function C).

1600

1200

Fitness
Fitness-β ⋅ #Nodes

1000

800

Fitness
Fitness-β ⋅ #Nodes

1400

600

600

Fitness

Fitness

Fitness

1200
800

1000
800

400

200

3

5

10

20

30

40

50

400

500
400
300

600

1

Fitness
Fitness-β ⋅ #Nodes

700

200
1

3

5

10

20

30

40

50

100

1

3

5

10

20

30

β

β

β

(a) Function A(b) Function B(c) Function C

40

50

Fig. 10. Influence of β on performance.

The proposed method uses CUDA to evaluate the similarity between individuals to provide acceleration; we set
the number of threads at |V |2 /2, where |V | is the numberof individuals. Therefore, the proposed method can adaptively reconstruct subpopulations without a substantial increase in computational effort. These results demonstrate
the effectiveness of the proposed method.
184

Finally, we investigated the difference in the solution
quality of the proposed, N, and F methods, as describedin Section 4.5. Thus, we evaluated the success rate of
subpopulation reconstructions of the three methods. In
Step P11 in Section 2, we describe an error procedurein the case where the subpopulation reconstruction fails.
These methods use Newman clustering to adaptively di-

Journal of Advanced Computational Intelligenceand Intelligent Informatics

Vol.25 No.2, 2021

Self-Organized Subpopulation in GP on GPU

Fig. 11. Scalability in terms of the number of individuals.

Table 1. Success rate and average of |Pm |.
MethodProposedN-method
F-method

Success rate [%]
100
30
100

Average of |Pm |Best trial Worst trial
44.5
65.7
–
–
56.1
65.2

vide the population P into subpopulations Pm based on
the features of individuals. If the value of each Hi j is similar, the number of communities extracted by Newman
clustering may be 1. Table 1 presents the success rateof subpopulation reconstruction over 20 trials. From thetable, it can be seen that the proposed and F-methods canreconstruct subpopulations with a 100% success rate, but
the N-method failed several times to reconstruct subpopulations. This result shows that the weighted network with
only the difference in the node size cannot create subpopulations, whereas the proposed method using multiple
features succeeded in creating subpopulations and exhibited a better performance than the F-method. We hypothesized that incorporating the multiple structure informationof individuals in the migration strategy is important forenhancing the combination of partial solutions and leads
to performance improvement. These results support our
hypothesis.

5. Conclusion

Acknowledgements
This work was supported by JSPS KAKENHI Grant Number
JP19H02032 and JP20K11982.

References:

To improve the subpopulation models in GP, it is essential to obtain a search strategy that has a good balance between local search and genetic diversity. We attempted to enhance subpopulation models by adaptively
changing subpopulations according to the features of individuals to improve local search, where individuals canbe represented by a tree structure in GP. In this paper, we
introduced a novel similarity H to reconstruct subpopulations based on the similarity of fitness values and alsonode size between individuals, and proposed a novel subpopulation model (SoS-GP). The proposed method generates a weighted network of individuals according to the
proposed measure H, and creates subpopulations using a
network clustering technique. In this case, considerable
computational effort is required to evaluate the similarity H. We proposed a CUDA-based model to shorten thisVol.25 No.2, 2021

time-consuming task. Moreover, we adopted a mutationin the proposed method. We hypothesized that mutation
can enhance the generation of various individuals withdifferent structures, and can maintain the balance between
local search and genetic diversity by cooperating with the
SoS model.
In well-known benchmark problems widely adopted instudies in the literature, we confirmed that the proposed
subpopulation model offers a significant performance advantage over the comparison method. We also analyzed
the behavior of the proposed method and demonstrated
its validity.
The experimental results indicated that the performance of the proposed method and F-method is high, andconfirmed that performance improvement can be achievedby incorporating multiple features and reconstructing subpopulations. Although in this paper we showed only the
results for three problems, we achieved similar results inother problems. Our immediate future work is to evaluate the proposed method using various other benchmark
problems. Moreover, other subpopulation models such as
ALPS require appropriate boundary parameters; however,in the proposed method, some parameters like the number of mutations γ and crossover rate are fixed as the firststep. We will investigate methods for adaptively controlling these parameters, and compare the proposed methodwith them. The proposed method uses fitness and node
size to evaluate the similarity between individuals for simplicity. However, various features exist for estimating thesimilarity between individuals. Therefore, future work includes combining these features to create a weighted network with the objective of further performance improvement. Additionally, we will attempt to develop a more
effective parallel method using CUDA to reduce computational effort.

[1] K. A. De Jong, “An Analysis of the Behavior of a Class of Genetic
Adaptive Systems,” Ph.D. thesis, University of Michigan, 1975.
[2] B. L. Miller and M. J. Shaw, “Genetic Algorithms with Dynamic
Niche Sharing for Multimodal Function Optimization,” Proc. of
1996 IEEE Int. Conf. on Evolutionary Computation (ICEC’96),
pp. 786-791, 1996.
[3] J. Hu, E. D. Goodman, and K. Seo, “Continuous Hierarchical Fair
Competition Model for Sustainable Innovation in Genetic Programming,” R. Riolo and B. Worzel (Eds.), “Genetic Programming Theory and Practice,” pp. 81-98, Kluwer Academic Publishers, 2003.
[4] J. Hu, E. D. Goodman, K. Seo, and M. Pei, “Adaptive Hierarchical
Fair Competition (AHFC) Model for Parallel Evolutionary Algorithms,” Proc. of the Genetic and Evolutionary Computation Conf.
2002 (GECCO 2002), pp. 772-779, 2002.
[5] G. S. Hornby, “ALPS: the age-layered population structure for
reducing the problem of premature convergence,” Proc. of the
8th Annual Conf. on Genetic and Evolutionary Computation(GECCO’06), pp. 815-822, 2006.
[6] G. S. Hornby, “Steady-state ALPS for real-valued problems,” Proc.
of the 11th Annual Conf. on Genetic and Evolutionary Computation(GECCO’09), pp. 795-802, 2009.
[7] J. M. Luna, J. R. Romero, C. Romero, and S. Ventura, “Discovering Subgroups by Means of Genetic Programming,” Proc. of the

Journal of Advanced Computational Intelligenceand Intelligent Informatics

185

Ono, K. and Hanada, Y.

16th European Conf. on Genetic Programming (EuroGP 2013),
pp. 121-132, 2013.
[8] Y. Jin, “Surrogate-assisted evolutionary computation: Recent advances and future challenges,” Swarm and Evolutionary Computation, Vol.1, Issue 2, pp. 61-70, 2011.
[9] D. Andre and J. R. Koza, “A Parallel Implementation of Genetic
Programming That Achieves Super-Linear Performance,” Proc. of
the Int. Conf. on Parallel and Distributed Processing Techniques andApplications (PDPTA’96), Vol.III, pp. 1163-1174, 1996.
[10] W. F. Punch, “How Effective Are Multiple Populations in Genetic
Programming,” Genetic Programming 1998: Proc. of the 3rd Annual Conf., pp. 308-313, 1998.
[11] M. Harwerth, “Experiments on Islands,” Proc. of the 14th European
Conf. on Genetic Programming (EuroGP 2011), pp. 239-249, 2011.
[12] S. Harding and W. Banzhaf, “Fast Genetic Programming and Artificial Developmental Systems on GPUs,” Proc. of the 21st Int.
Symp. on High Performance Computing Systems and Applications(HPCS’07), doi: 10.1109/HPCS.2007.17, 2007.
[13] S. Harding and W. Banzhaf, “Genetic programming on GPUs forimage processing,” Int. J. of High Performance System Architecture, Vol.1, Issue 4, pp. 231-240, 2008.
[14] I. Arnaldo, K. Veeramachaneni, and U.-M. O’Reilly, “Flash:
A GP-GPU Ensemble Learning System for Handling Large
Datasets,” Proc. of the 17th European Conf. on Genetic Programming (EuroGP 2014), pp. 13-24, 2014.
[15] W. B. Langdon and M. Harman, “Genetically Improved CUDA C++
Software,” Proc. of the 17th European Conf. on Genetic Programming (EuroGP 2014), pp. 87-99, 2014.
[16] W. B. Langdon, “A Many Threaded CUDA Interpreter for Genetic
Programming,” Proc. of the 13th European Conf. on Genetic Programming (EuroGP 2010), pp. 146-158, 2010.
[17] K. Ono and Y. Hanada, “A CUDA-based self-adaptive subpopulation model in genetic programming: CuSASGP,” Proc. of the 2015
IEEE Congress on Evolutionary Computation (CEC), pp. 15431550, 2015.
[18] M. Girvan and M. E. J. Newman, “Community structure in social
and biological networks,” Proc. of the National Academy of Scienceof the United States of America, Vol.99, No.12, pp. 7821-7826,
2002.
[19] A. Clauset, M. E. J. Newman, and C. Moore, “Finding community structure in very large networks,” Physical Review E, Vol.70,
Issue 6, doi: 10.1103/PhysRevE.70.066111, 2004.
[20] K. Yanai and H. Iba, “Estimation of distribution programming basedon Bayesian network,” Proc. of the 2003 Congress on Evolutionary
Computation (CEC’03), Vol.3, pp. 1618-1625, 2003.
[21] W. B. Langdon and R. Poli, ”Foundations of Genetic Programming,” Springer, 2002.
[22] W. B. Langdon and R. Poli, “Why Ants are Hard,” Genetic Programming 1998: Proc. of the 3rd Annual Conf., pp. 193-201, 1998.

Name:
Keiko Ono

Affiliation:
Doshisha University

Address:
1-3 Tatara Miyakodani, Kyotanbe, Kyoto 610-0394, Japan

Brief Biographical History:
2010- Assistant Professor, Ryukoku University
2014- Associate Professor, Ryukoku University
2020- Associate Professor, Doshisha University

Main Works:• K. Ono, Y. Hanada, M. Kumano, and M. Kimura, “Enhancing Island
Model Genetic Programming by Controlling Frequent Trees,” J. of
Artificial Intelligence and Soft Computing Research, Vol.9, No.1,
pp. 51-65, 2019.

Membership in Academic Societies:• The Institute of Electrical and Electronics Engineers (IEEE)• Information Processing Society of Japan (IPSJ)• The Japanese Society for Evolutionary Computation (JPNSEC)

Name:
Yoshiko Hanada

Affiliation:
Kansai University

Address:
3-3-35 Yamate-cho, Suita, Osaka 564-8680, Japan

Brief Biographical History:
2008- Assistant Professor, Kansai University
2016- Associate Professor, Kansai University

Main Works:• K. Matsumura, Y. Hanada, and K. Ono, “Probabilistic Model-BasedMultistep Crossover Considering Dependency Between Nodes in Tree
Optimization,” R. Lee (Ed.), “Software Engineering, Artificial
Intelligence, Networking and Parallel/Distributed Computing,”
pp. 187-200, Springer, 2017.

Membership in Academic Societies:• The Institute of Electrical and Electronics Engineers (IEEE)• Information Processing Society of Japan (IPSJ)• The Institute of Energy Economics, Japan (IEEJ)

186

Powered by TCPDF (www.tcpdf.org)

Journal of Advanced Computational Intelligenceand Intelligent Informatics

Vol.25 No.2, 2021