Improved Differential Evolution based Feature Selection through Quantum, Chaos, and Lasso

Abstract

Modern deep learning continues to achieve outstanding performance on an astounding variety of high-dimensional tasks. In practice, this is obtained by fitting deep neural models to all the input data with minimal feature engineering, thus sacrificing interpretability in many cases. However, in applications such as medicine, where interpretability is crucial, feature subset selection becomes an important problem. Metaheuristics such as Binary Differential Evolution are a popular approach to feature selection, and the research literature continues to introduce novel ideas, drawn from quantum computing and chaos theory, for instance, to improve them. In this paper, we demonstrate that introducing chaos-generated variables, generated from considerations of the Lyapunov time, in place of random variables in quantum-inspired metaheuristics significantly improves their performance on high-dimensional medical classification tasks and outperforms other approaches. We show that this chaos-induced improvement is a general phenomenon by demonstrating it for multiple varieties of underlying quantum-inspired metaheuristics. Performance is further enhanced through Lasso-assisted feature pruning. At the implementation level, we vastly speed up our algorithms through a scalable island-based computing cluster parallelization technique.

1 Introduction

The fundamental objective of feature selection is to identify the most important and discriminative features from a given set of features. Its prominence has garnered the attention of researchers and practitioners in domains replete with big, high-dimensional datasets. Feature selection [4, 46] entails a great reduction in computational complexity, improvement in human comprehensibility, and easy deployment of models in production. Wrapper methods based on metaheuristics solve the feature subset selection (FSS) problem by posing it as a combinatorial optimization problem. These methods attempt to identify an efficient feature subset with the least cardinality and associated high accuracy among the $2^{n}-1$ possible number of combinations, where $n$ is the total number of features in the dataset. Among the metaheuristics employed for this purpose, evolutionary algorithms have been proven to be most efficient for determining optimal feature subsets owing to the inherent parallelism present in the population-based search [27, 45].

Table 1: Notation used in the current study

$\phi$	Empty set
$N$	Population size
$X$	Population comprising N number of solutions
$n$	Number of dimensions / features
$X_{i}$	ith solution of population X and having $n$ dimensions
$M^{t}$	Mutated vectors at generation $t$ comprising $N$ solutions
$M_{i}^{t}$	$i\textsuperscript{th}$ solution of mutation population of $n$ dimensions
F	Mutation factor
$U^{t}$	trial vectors at generation t comprising ps solutions
$U_{i}^{t}$	$i\textsuperscript{th}$ solution of mutation population of $n$ dimensions
$CR$	Crossover rate
$MAXITR$	Maximum number of iterations
$randi$	Randomly chosen index
$d_{t}$	Chaotic number at tth time step
$rand(0,1)$	Random number generated between $[0,1]$
$\lambda$	Logistic map control parameter
AUCi	AUC score of an ith solution
cardinalityi	Cardinality score of an $i\textsuperscript{th}$ solution
$cvalue_{t}$	Chaotic random number at $t\textsuperscript{th}$ time step
$P$	RDD of population X
$Xtrain$	Train dataset
$Xtest$	Test dataset
$localN$	Local population size
$mMig$	Maximum number of migrations
$mGen$	Maximum number of generations
$localN$	Local population / sub-population pertained to a data island comprising lps solutions

Among the evolutionary algorithms, Differential Evolution (DE) proved to be robust while solving many combinatorial and continuous optimization problems [9, 8, 32]. Despite its supremacy over other algorithms, DE suffers from a tendency to get stuck prematurely in local optimal solutions, which affects its exploration and exploitation capabilities [9]. To alleviate these issues, researchers proposed several quantum-inspired algorithms (QIEA) to solve a myriad of combinatorial optimization problems such as knapsack problem, truck trail problem, and portfolio optimization. Often, QIEAs accelerate the evolution process owing to their quantum parallelization [38, 19] and entanglement of the quantum state. Further, the diversity and convergence rate are improved well enough to increase the probability of getting better global optimal solutions.

In today’s data-rich environment, the humongous growth of high-dimensional datasets warrants the critical need for developing scalable algorithms [17, 36]. Despite the popularity of Hadoop and its ecosystem as a big data framework, Spark’s unique features, including in-memory computing and seamless integration, have made it a viable and competitive alternative to Hadoop [47]. However, the extant Quantum-Inspired Evolutionary Algorithms (QIEAs) are not scalable to large, high-dimensional datasets.

Our contributions in this paper include:

• proposing chaotic, quantum-inspired evolutionary algorithms for FSS in high-dimensional data and demonstrating their superiority over the extant methods on four problems.

• utilizing the Lyapunov exponent to ensure that we work in a truly chaotic regime.

• introducing LASSO LR in place of LR as a classifier for the FSS wrapper

• proposing parallel versions of the above-mentioned algorithms operating in a single-objective environment under an island-based approach in Apache Spark framework.

The paper is structured as follows: Section 2 reviews the literature, Section 3 presents our new algorithms, Section 4 describes the datasets analyzed in the study, Section 5 analyses the results obtained, and section 6 concludes the paper.

Table 2: Full form of the acronyms of the Algorithms employed in the current study

No.	Algorithm	Description
1	LR	Logistic Regression
2	LLR	Least Absolute Shrinkage and Selection Operator LR
3	BDE	non-quantum counterpart of QBDE
4	QBDE-I	a non-gate quantum variant of BDE with threshold trick and random numbers generated from the uniform distribution
5	QBDE-II	a gate quantum variant of BDE with threshold trick and random numbers generated from uniform distribution
6	CQBDE-I	a chaotic maps guided variant of QBDE-I but without Lyapunov exponent guidance
7	CQBDE-II	a chaotic variant of QBDE-II without Lyapunov exponent guidance
8	CLQBDE-I	a chaotic variant with Lyapunov exponent of QBDE-I
9	CLQBDE-II	a chaotic variant with Lyapunov exponent of QBDE-II
10	CQIEA	an algorithm from [35]
11	CTQIEA	a variant of CQIEA [35]which integrates threshold trick and chaos without Lyapunov exponent guidance
12	CLTQIEA	a variant of CQIEA [35] which incorporates both threshold trick and Lyapunov exponent guidance.

(a) Schematic of the BDE

(b) Schematic of the proposed QBDE variants

(c) Schematic of the proposed CQBDE variants

(d) Block diagram of the proposed wrapper

Figure 1: Generic framework of the proposed wrappers

2 Background and Literature Review

2.1 Differential Evolution

Binary Differential Evolution (BDE), a stochastic population-based global optimization algorithm, starts by initializing the random population, consisting of $N$ candidate solution vectors ($X_{i}$), where $N$ is the population size. This candidate solution vector follows the binary encoding scheme. Each candidate solution vector is subjected to all of the following three heuristics in each iteration (or generation) of the algorithm (see Fig. 1(a)).

At each generation $t$, the candidate solution vector ($X_{i}^{t}$) within $n$ dimensional search space, is subjected to the mutation operation yielding the mutant vector ($M_{i}^{t}$). The mutation operation is applied as presented in Eq. (1).

$$M_{i}^{t}=X_{i1}^{t}+F*(X_{i2}^{t}-X_{i3}^{t})$$ (1)

where $X_{i1}^{t}$,$X_{i2}^{t}$ and $X_{i3}^{t}$ are three randomly chosen distinct vectors from the current generation $t$. $F$, the mutation factor, is a user-defined parameter, and lies in the range $[0,1]$. After this, the mutant vector may not be binary anymore. Hence, sigmoid based discretization process (see Eq. (2)) is applied to every $m_{ij}^{t}$ ($j^{th}$ member of the $M_{i}^{t}$) thereby converting continuous vector into a binary vector.

$$m_{ij}^{t}=\begin{cases}1,&\text{if }\text{rand}(0,1)<\text{sigmoid}(m_{j}^{t})\\ 0,&\text{else}\end{cases}$$ (2)

Then, the discretized mutant vector is subjected to crossover operation where it is subjected to the mating with the corresponding candidate solution vector to generate the trial vector. The crossover operation is applied to trial vector $U_{i}^{t}$, as presented in Eq. (3).

$$u_{ij}^{t}=\begin{cases}m_{ij}^{t},&\text{if }\text{rand}(0,1)<CR\text{ and }j\neq\text{$randi$}\\ x_{ij}^{t},&\text{if }\text{rand}(0,1)\geq CR\text{ and }j\neq\text{$randi$}\end{cases}$$ (3)

where $j=1,2,…n,u_{ij}^{t}$ is the $j^{th}$ bit of $U_{i}^{t}$, $rand(0,1)$ is the random number generated in the interval $[0,1]$ from a uniform distribution. $randi$ is a randomly chosen index to make sure that the generated trial vector is different from the mutant vector. CR represents the crossover rate, is a user-defined parameter, and lies in the range $[0,1]$.

Finally, the fitness score is computed for the trial vectors. Then, the selection operation is applied by comparing the corresponding target vectors and trial vector to produce an offspring. Better solutions survive and forms the parent population for the subsequent iteration. The selection operation follows the rule as presented in Eq. (4):

$$X_{i}^{(t+1)}=\begin{cases}X_{i}^{(t)},&\text{if }f(X_{i})>f(U_{i})\\ U_{i}^{(t)},&\text{otherwise}\end{cases}$$ (4)

As mentioned earlier, this is continued till the completion of maximum iterations or other convergence criteria, if any, are met.

2.2 Quantum Computing

To exploit the notions of the quantum theory within classical computers, quantum-inspired algorithms are proposed by [19]. It employs quantum mechanics concepts such as quantum measurement, superposition of states, inference, and entanglement. A quantum bit is the basic unit of information in quantum computation defined by the linear combination of the two states as given in Eq. (5).

$$Q=\alpha\ket{0}+\beta\ket{1}$$ (5)

The coefficients of $\alpha$ and $\beta$ are two complex numbers that must satisfy the norm relation as given in Eq. (6).

$$|\alpha|^{2}+|\beta|^{2}=1$$ (6)

where the probability of the observing state $\ket{0}$ is $|\alpha|^{2}$ and the probability of the state $\ket{1}$ is $|\beta|^{2}$. A quantum register is composed of $n$ qubits containing $2^{n}$ possible values simultaneously owing to the superposition of states.

2.3 Overview of Chaos Theory

The theory of chaos originated in the 1800s and was further developed by [7] to tackle challenges in complex non-linear systems [30]. Chaotic systems are dynamic and deterministic, evolving from initial conditions, with trajectories describing the system states in the state space. Chaos heavily relies on the initial conditions and exhibits two properties: ergodicity and intrinsically stochastic nature. Chaotic maps generate the sequence of numbers that exhibit these chaotic properties, aiding EAs in escaping local minima [31, 25, 26]. It is noticed that these chaotic maps are introduced to dynamically adjust the hyperparameters, and enhance adaptability to handle evolutionary dynamics optimally. It also facilitates searching in the regions which are left out by the random sequence. They enable exploration in regions left out by random sequences and have been extensively studied in large-scale continuous optimization problems. Chaotic maps produce a series of numbers from a probability distribution that differs from a uniform (0,1) distribution, exhibiting deterministic randomness. This deterministic nature allows predicting the sequence of numbers generated, as they are governed by differential equations and subject to the initial conditions. Now, a well-known chaotic map logistic map which is used in the current study is discussed below:

Logistic map [28]: The Logistic map exhibits chaotic behaviour in a discrete-time demographic model. This is a polynomial mapping of degree 2. The mathematical representation is defined in Eq. (7).

$$d_{t+1}=\lambda*(d_{t}*(1-d_{t}))$$ (7)

Here, the constant $\lambda$ lies in the range of $[0,4]$ and determines the behaviour of this Logistic map. In the current research study, $\lambda$ = 4 is chosen.

2.4 Literature Review

The first quantum-inspired evolutionary algorithm (QEA) was proposed by Han and Kim for the knapsack problem in the pioneering paper, [19]. Their approach employed qubit notation and a rotation gate together with the migration operation, strategically guiding the solution to reach a global optimum. [22] proposed an adaptive quantum differential evolution (AQDE), which dynamically adjusts the mutation and cross-over rates based on their success streaks. It outperformed QEA in solving the knapsack problem. [39] proposed elitism-based quantum differential evolution for FSS on small datasets. [13] proposed Multi-strategy quantum differential evolution, which integrates multiple strategies, including a mutation strategy with difference vector, multi-population mutation, and adaptive rotation angle state. Other quantum evolutionary algorithms studied in the literature include (i) hybrid of QDE (QDE acronym is not defined anywhere) and grey wolf optimizer [43] tailored for the knapsack problem, (ii) multi-objective quantum-inspired hybrid DE [24], which is a hybrid of genetic algorithm quantum variants and DE for multi-objective next-release problems, and (iii) vector hop algorithm [18], a hybrid of differential evolution and particle swarm optimization. [12] proposed an improved quantum differential evolution by incorporating the principles of the divide-and-conquer concept of a cooperative coevolutionary algorithm, which improved both exploration and exploitation capabilities. (iv) Another hybrid DE where the first stage invoked QDE and the resultant population is passed on to BDE in the second stage, which continues the evolution process. All these hybrids improved the search process and obtained better convergence capabilities than their individual constituents in the standalone mode.

We now briefly survey the parallel and distributed versions of DE [20, 33, 37, 48, 40] developed across varied environments like Spark, CUDA, MPI, and OpenMP. [48] proposed two master-slave-based parallel approaches, namely, (i) a data-based MapReduce model and (ii) a population-based MapReduce model. [40] introduced two parallel strategies, master-slave and island approaches, evaluated on the AWS Spark cloud and tested the performance on benchmark functions. [6] and [5] developed master-slave-based parallel DE approaches for large-scale clustering and cluster optimization problems, respectively. [2] introduced a cost-sensitive DE classifier (SCDE) based on Euclidean distance for imbalanced classification datasets. [1] introduced a fine-grained parallel DE under the OpenMP framework for optimal networking, aiming to reduce the computational load on mappers and reducers. [11] introduced Parallel DE (PDE) under the Spark environment, showing promising speedup. [21] presented SgtDE, a grouping topology model for large-scale optimization, achieving significant speedup. [10] developed parallel DE under CUDA, while [44] designed a self-adaptive DE framework in CUDA for benchmark functions. Further, several parallel versions of DE are proposed to solve resource allocation problems [3, 15, 14], hydro scheduling [16], large-scale clustering [23], optimized workflow placement [41], and multi-objective flow scheduling problems [34] as well. Table A1 of the Appendix captures a brief overview of the literature.

Figure 2: Generic schematic diagram of the island model based wrapper

3 Proposed Methodology

This section introduces the objective function employed in the current study, followed by the algorithm, and an overview of the proposed parallel mechanism, which applies to all algorithms discussed. All the acronyms employed in the study are presented in Table 1.

3.1 Objective Function

The objective function considered in this study is the Area under the receiver operating characteristic Curve (AUC) Eq.(8). It is defined as the mean of specificity and sensitivity. We specifically considered AUC due to its proven robustness while handling imbalanced datasets.

$$\text{AUC}=\frac{\text{Sensitivity}+\text{Specificity}}{2}$$ (8)

Where sensitivity (refer to Eq. (10)) is the ratio of the positive samples that are correctly predicted to be positive to all the predicted positive samples. This is also called the True Positive Rate (TPR).

$$\text{Sensitivity}=\frac{TP}{TP+FN}$$ (9)

Where $TP$ and $FN$ are the true positive and false negative, respectively. Specificity (refer to Eq. (10)) is the ratio of the negative samples that are correctly predicted to be negative to all the predicted negative samples. This is also called True Negative Rate (TNR).

$$\text{Specificity}=\frac{TN}{TN+FP}$$ (10)

Where $TN$ and $FP$ are the true negative, and false positive, respectively.

This study employs LR and LLR as the classifiers for the proposed FSS wrapper because their training converges quickly. Unlike LR, LLR performs regularization and reduces the number of features too.

3.2 Proposed Chaotic Quantum Inspired Algorithm

In this study, we propose two distinct algorithms based on chaotic and quantum principles: (i) non-gate variant, one focusing on information exchange across different qubits without gates (CLQBDE-I), and (ii) gate variant, the second involving information exchange via gates (CLQBDE-II) (see Table 2 and Fig. 1(b), Fig. 1(c)).

Unlike the chaos-based algorithms proposed in the literature [29, 4, 46], the Lyapunov exponent guides our population initialization method, ensuring that the wrapper based FSS operates within a true chaotic regime. Consequently, chaotic numbers are not utilized from the initial time step of the logistic chaotic map, unlike the existing approaches. Instead, they are introduced after a specific number of time steps, which are determined by the Lyapunov exponent. After thorough experimentation, we discovered that the logistic map enters a truly chaotic regime after $5000$ time steps. This is observed offline and employed within the algorithm. Histogram plots comparing the distributions of numbers generated from the uniform random distribution and logistic chaotic map guided by Lyapunov exponent from 1 to 5000 time steps are presented in Fig. A1 of the Appendix. It is crucial to note that in the current study, Lyapunov exponent-guided logistic map is introduced as the sequence generated by chaotic numbers is observed to be in true chaotic regime after first $5000$ timesteps. Consequently, both variants, CQBDE-I and CQBDE-II, incorporate chaotic numbers in the initialization for each qubit’s alpha and beta states for $i^{th}$ solution, satisfying Eq. (11).

$$Q^{(i)}=\begin{bmatrix}\alpha^{(i)}_{1}&\alpha^{(i)}_{2}&\dots&\alpha^{(i)}_{d}\\ \beta^{(i)}_{1}&\beta^{(i)}_{2}&\dots&\beta^{(i)}_{d}\end{bmatrix}$$ (11)

$$\left|\alpha^{(i)}_{j}\right|^{2}+\left|\beta^{(i)}_{j}\right|^{2}=1$$ (12)

where $n$ is the number of dimensions, $i$ is the index of the solution in the population of size $N$.

The quantum representation of each population member is given by Eq. (11) where the components satisfy the condition in Eq. (12). It consists of $N$ candidate solution vectors ($X_{i}$), where $n$ is the population size and each solution vector is of size $2*n$, where $n$ is the number of features. Often thus, formed quantum state vectors need to be converted into binary encoded solutions for combinatorial optimization problems. In the literature, the following rule is generally adopted.

$$x^{(i)}_{j}=\begin{cases}1,&\text{if }\text{rand}(0,1)<\left|\beta^{(i)}_{j}\right|^{2}\\ 0,&\text{otherwise}\end{cases}$$ (13)

where $i$ and $j$ are the solution and feature indices, respectively.

However, we observed that as the number of dimensions ($n$) increases, more features are being selected. This observation prompted us to introduce a constraint-based quantum state-to-binary solution conversion. The modified rule is defined as follows (refer to Eq. (14)). This Equation is called as the threshold trick. This is employed for both the QBDE-I and QBDE-II variants.

$$x_{i}=\begin{cases}1,&\text{if }\text{rand}(0,1)<\left|\beta^{(i)}_{j}\right|^{2}\text{ and }d<\theta\\ 0,&\text{otherwise}\end{cases}$$ (14)

Here, $\theta$ represents a user-defined parameter within the range of $[0,1]$, $d$ is the number generated from rand(0,1).

We proposed two different variants of chaotic quantum algorithms guided by the Lyapunov exponent and threshold trick. The distinction between them is as follows: (i) CLQBDE-I is a non-gate variant, and (ii) CLQBDE-II is a gate-variant where a rotation gate is employed in place of quantum mutation.

3.2.1 CLQBDE-I

The proposed Chaotic Quantum Binary Differential Evolution-I (CLQBDE-I) begins by generating the quantum matrix as described in Eq. (11). For each candidate solution $i$, the $\alpha^{(i)}_{j}$ variables are generated from the chaotic series, while the corresponding $\beta^{(i)}_{j}$ variables are obtained by using the relation Eq. (12). As mentioned earlier, the first step of the chaotic series is guided by the Lyapunov exponent. Once generated, the binary encoded solution is produced using the threshold trick outlined in Eq. (14). The fitness score is calculated by training the classifier, LR or LLR.

$$\small\begin{bmatrix}\alpha^{(i)}_{j}\\ \beta^{(i)}_{j}\end{bmatrix}=\begin{bmatrix}\alpha^{(x1)}_{j}+F\cdot(\alpha^{(x2)}_{j}-\alpha^{(x3)}_{j})\\ \beta^{(x1)}_{j}+F\cdot(\beta^{(x2)}_{j}-\beta^{(x3)}_{j})\end{bmatrix}$$ (15)

Then, the initial population undergoes mutation operation (refer to Eq. 15) effected by each candidate solution’s $\alpha$ and $\beta$ states by randomly selecting the three solutions. This operation is performed for a pre-specified population size.

$$\begin{bmatrix}\alpha^{(u)}_{j}\\ \beta^{(u)}_{j}\end{bmatrix}=\begin{cases}\begin{bmatrix}\alpha^{(m)}_{j}\\ \beta^{(m)}_{j}\end{bmatrix},&\text{if }d<CR\text{ and }j\neq\text{$randi$}\\ \\ \begin{bmatrix}\alpha^{(u)}_{j}\\ \beta^{(u)}_{j}\end{bmatrix},&\text{otherwise}\end{cases}$$ (16)

where $d$ is the number generated from rand(0,1).

Subsequently, the mutated quantum matrix of each of the solutions undergoes crossover operation as described in Eq. 16. Selection is then performed by combining parent and offspring populations and sorting them based on fitness scores. The resulting population becomes the parent population for the next iteration. This iterative process continues for a pre-specified number of iterations, resulting in the convergence of the algorithm, and the resultant population is evaluated on the test data, which remains unchanged throughout the evolution process.

3.2.2 CLQBDE-II

We proposed another variant, gate-based quantum binary differential evolution, and named it chaotic quantum binary differential evolution-II (CQBDE-II). The distinction between them lies in employing the rotation gate in place of the quantum mutation operation to generate the mutated quantum matrix.

The population is initialized The rotation gate is introduced to mimic the behaviour of the mutation gate. These rotation gates are represented as unitary matrices, which are employed to rotate the state of a qubit based on the rotation angle. We followed the look-up table obtained from [19] to decide the rotation angle. The rotation gate is multiplied with each qubit of the corresponding qubit in the candidate solution and thus the mutated vector is generated. Thissolution vector is then subjected to the crossover operation, selection, and threshold trick. These steps are identical to that of CLQBDE-I as discussed in the sub-section 3.2.1. The rotation gate employed in this study is presented in Eq. 17.

$$U(\Delta\theta)=\begin{bmatrix}\cos(\Delta\theta)&-\sin(\Delta\theta)\\ \sin(\Delta\theta)&\cos(\Delta\theta)\end{bmatrix}$$ (17)

Table 3: Encoding Scheme of the population

Key	Value:
$Key_{1}$	$\langle$ $BinaryVector_{1}$, $QuantumMatrix_{1}$, $trainedModelCoef_{1}$, $AUC_{1}$ $\rangle$
$Key_{2}$	$\langle$ $BinaryVector_{2}$, $QuantumMatrix_{2}$, $trainedModelCoef_{2}$, $AUC_{2}$ $\rangle$
…	…
$Key_{N}$	$\langle$ $BinaryVector_{N}$, $QuantumMatrix_{N}$, $trainedModelCoef_{N}$, $AUC_{N}$ $\rangle$

3.3 Proposed Parallel Approach

The proposed parallel algorithm adopts an island approach, which divides the data into several partitions known as data islands.

Table 3 depicts the population schema maintained for all approaches in this work. The population consists of a solution of size N. Each solution has two different fields: (i) Key field having the unique id information to identify the solution uniquely, (ii) The Value field has the following subfields: (a) Binary vector: which is of length number of features, nfeat is a binary vector where the presence of ’0’ says that a particular feature is not selected, and ’1’ says that the particular feature is selected. (b) Quantum matrix: Here, the quantum state information of the corresponding is stored.(c) Trained model coefficients: as we know that in the wrapper methods, a classifier is chosen to evaluate the solution’s performance. Thus trained model coefficients are stored in this sub-field. The main reason for storing them is to use them in the test phase. (d) AUC: The trained model results in an AUC for each solution. That information is stored in the sub-field. This kind of schema makes sure to preserve the solution-related information in a single space.

The proposed parallel approach comprises three main phases: Initialization, Training phase, and Test phase. As illustrated in Fig. 2, there are $k$ sub-populations and $k$ data islands. The parallel algorithm operates at two levels: driver and worker. The corresponding algorithms are presented in Algorithms 1 and 2 (see Appendix), respectively. All these phases, along with migration rule invocation, happen at the driver node, while the evolution of the parallel EAs happens during the training phase, and evaluation during the test phase is executed at the worker node. However, filtering top solutions post-migration and the aggregation of the final results during the test phase is performed at the driver node.

Phase-I: Initialization In this phase, the quantum matrix for the population size is chaotically initialized using the biased sampling method [42] depending on the underlying metaheuristic (i.e., DE / EA). This population follows the data structure as presented in Table 3. Subsequently, binary-encoded solutions are obtained using a threshold trick. Thus initialized population is treated as a global population and broadcasted along with hyperparameters from the driver to the worker nodes.

Phase-II: Training Phase A miniature EA is evolved in parallel in each data island at each worker node. The block diagram of miniature wrapper is depicted in Fig. 1(d). A sub-population/local population of size $localN$ ($<N$) is initially extracted by following random sampling with replacement. The $train-and-update$ phase begins, where fitness scores are evaluated using the classifier, LR or LLR, followed by the updation of the quantum matrix, binary encoded solutions, classifier coefficients, and AUC in the respective fields (see Table 3). This $train-and-update$ process is common for all the algorithms. Heuristics specific to each algorithm are applied to the quantum matrix of the population, resulting in the generation of an offspring quantum matrix. $Train-and-update$ operations follow, and the top $localN$ solutions are retained after selection. This process is repeated for a pre-specified maximum number of generations ($mGen$), resulting in $k*lps$ solutions from each of the $k$ data islands. Upon completing these steps, the control returns to the driver, where the migration policy is invoked.

Migration Policy Solutions are sorted based on the fitness scores, and the top $N$ solutions are selected. After invoking the migration policy, the worker algorithm is again executed using the population obtained post-migration which is distributed to worker nodes. Subsequently, the sub-population corresponding to a single island is selected by following random sampling with replacement and the evolution process continues. This process is repeated for a pre-specified number of migrations ($mMig$).

Phase-III: Test Phase In this phase, the converged population obtained in the training phase is evaluated on the test dataset. That means, AUC is computed on the test dataset using the coefficients corresponding to each solution. Subsequently, test fitness scores are computed, and the results are reported accordingly. Notably, while the proposed parallel metaheuristics differ in their respective heuristics, they adhere uniformly to the same train-and-update step, migration policy invocation, and test phase.

4 Results & Discussions

All experiments are conducted on a 5-node cluster configuration, comprising a driver node and 4 slave nodes, each equipped with Intel i7 $9^{th}$ generation processors and 32GB of RAM. Notably, the driver node serves the dual role as the driver node as well as the worker node. The benchmark datasets analyzed in this study are briefly described in Table 4, while Table 5 presents the optimal hyperparameters identified through meticulous fine-tuning following grid-search. A train-test split ratio of 80%:20% is followed using stratified random sampling to ensure equal proportion of representation of classes in both datasets. All algorithms are executed for 20 runs in order to mitigate the impact of random seed variation on the algorithm, which is the standard practice followed in the evolutionary algorithms (EAs) literature. For all the experiments, the number of migrations is fixed at one. The top solution in each run is identified as the one achieving the highest AUC. Thus, we get 20 top solutions, one each for 20 runs. Mean AUC and mean cardinality of these 20 top solutions are computed and presented in Tables 6-10.

The proposed different quantum algorithms are compared against a non-quantum counterpart, BDE, and the Chaotic Quantum-Inspired Evolution Algorithm (CQIEA)[35]. To assess the effectiveness of the proposed method, we developed the variants listed in Table 2.

The reasons for developing these variants are as follows: The algorithms QBDE-I and QBDE-II are meant to see the effectiveness of introducing quantum operators; CQBDE-I and CQBDE-II are designed to see the effectiveness of chaos; CLQBDE-I and CQBDE-II are designed to see the role of the Lyapunov exponent when coupled with chaotic maps. Further, CTQIEA and CLTQIEA are introduced to see the effectiveness of the threshold trick and threshold trick in combination with the Lyapunov exponent, respectively.

Table 4: Dataset Information

Dataset Name	# Objects	# Features	Size of dataset
Epsilon	500,000	2,000	10.8 GB
IEEE Malware	1,500,000	1,000	3.2 GB
OVA_Omentum	1,584	10,935	108.3 MB
OVA_Uterus	1,584	10,935	108.3 MB

Table 5: Hyper parameters employed in the current study

Dataset Name	$N$ ($localN$)	Total generations $mGen$ per migration	CR	MR	$\theta$
Epsilon	30 (15)	20 (10)	0.90	0.80	0.10
IEEE Malware	30 (15)	10 (5)	0.90	0.80	0.15
OVA_Omentum	300 (200)	20 (10)	0.90	0.8	0.01
OVA_Uterus	300 (200)	20 (10)	0.90	0.8	0.01

Table 6: Results of BDE variants

	Datasets
Algorithm	Epsilon		IEEE Malware		OVA_Omentum		OVA_Uterus
	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$
BDE+LR	1321.3	0.835	646.50	0.802	876.26	0.836	63.35	0.788
BDE+LLR	1323.7	0.846	505.7	0.817	112.65	0.833	108.65	0.802

Table 7: Results of QBDE variants

	Datasets
Algorithm	Epsilon		IEEE Malware		OVA_Omentum		OVA_Uterus
	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$
QBDE-I+LR	204.7	0.742	121.25	0.753	159.25	0.800	69.90	0.807
QBDE-II+LR	214.75	0.745	127.30	0.758	70.2	0.846	70.65	0.808
QBDE-I+LLR	248.95	0.753	168.5	0.826	168.5	0.826	107.75	0.796
QBDE-II+LLR	140.85	0.783	140.6	0.793	107.0	0.944	108.4	0.799

Table 8: Results of CQBDE variants

	Datasets
Algorithm	Epsilon		IEEE Malware		OVA_Omentum		OVA_Uterus
	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$
CQBDE-I+LR	269.6	0.522	142.05	0.80	794.7	0.5	830.05	0.5
CQBDE-II+LR	143.85	0.780	144.5	0.791	107.5	0.905	106.45	0.842
CQBDE-I+LLR	397.95	0.797	162.55	0.818	159.25	0.808	108.75	0.796
CQBDE-II+LLR	138.65	0.784	145.85	0.799	110.35	0.855	101.7	0.808

Table 9: Results of CLQBDE variants

	Datasets
Algorithm	Epsilon		IEEE Malware		OVA_Omentum		OVA_Uterus
	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$
CLQBDE-I+LR	273.85	0.523	143.8	0.806	798.7	0.5	798.7	0.5
CLQBDE-II+LR	141.5	0.777	140.6	0.793	107	0.944	107.15	0.853
CLQBDE-I+LLR	145.95	0.755	145.95	0.796	156.95	0.923	151.0	0.865
CLQBDE-II+LLR	141.70	0.784	140.9	0.801	110.25	0.954	112.2	0.812

Table 10: Results of CQIEA variants

	Datasets
Algorithm	Epsilon		IEEE Malware		OVA_Omentum		OVA_Uterus
	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$	$f_{1}$	$f_{2}$
CQIEA+LR	1968.10	0.841	885.5	0.802	5033.0	0.788	6942.9	0.792
CTQIEA+LR	421.0	0.790	217.65	0.779	240.6	0.947	241.1	0.884
CLTQIEA+LR	412.7	0.790	208.3	0.812	290.0	0.819	280.95	0.794
CQIEA+LLR	888.0	0.853	677.9	0.810	8356.25	0.837	6942.90	0.792
CTQIEA+LLR	256.20	0.815	190.8	0.798	261.50	0.842	318.8	0.806
CLTQIEA+LLR	258.50	0.815	192.05	0.801	258.70	0.893	332.3	0.803

Table 11: Results of the paired t-test

Dataset	Model (Top1 vs Top2 w.r.t AUC)	t-statistic	p-value
Epsilon	CTQIEA + LLR* vs	0.421	0.675
	CLTQIEA+LLR
IEEE Malware	QBDE-I+LLR* vs	1.204	0.235
	CQBDE-I + LLR
OVA_Omentum	CLQBDE-II + LLR* vs	1.14	0.260
	CTQIEA + LR
OVA_Uterus	CTQIEA + LR vs	3.84	0.0004
	CLQBDE-I + LLR**
*Better method based on cost-benefit analysis.
**Better method based on statistical significance.

4.1 Comparative Analysis Ablation Study

We conducted our ablation study in a boosted manner. That means in the first instance, parallel BDE is developed as a baseline and then quantum-inspired versions of BDE are developed (boosted wrapper-I), and finally, chaotic versions of these quantum-inspired algorithms are developed (boosted wrapper-II).

Initially, non-quantum variants, i.e., BDE variants, are developed as a baseline, and the corresponding results are presented in Table 6. With an aim to remove irrelevant features, we employed LLR in place of LR as a classifier. LLR demonstrated its ability to weed out unimportant features in the following datasets: (i) In the IEEE Malware dataset, BDE with LLR as the classifier, obtained mean AUC, and there was a reduction of 150 features. (ii) Similarly, in the OVA_Omentum dataset, LLR led to a great reduction (around 600 features) in terms of mean cardinality, with slightly improved mean AUC. However, this observation is not consistently noticed across the remaining two datasets. In the other datasets (Epsilon and OVA_Uterus), the mean AUC is slightly improved. However, it is accompanied by a modest increase in the mean cardinality.

The strength of LLR over LR is that the latter retains unimportant features during the evolution process, leading to increased cardinality of the solutions.

In the boosted wrapper-I, to obtain better optimal results, we developed their quantum-inspired counterparts, namely, QBDE-I and QBDE-II, where the QBDE-I variant has gates, and the QBDE-II variant has no gates. Consequently, we achieved lower mean cardinality accompanied by higher mean AUC across all datasets except the Epsilon dataset (see Table 7).

Then, in the boosted wrapper-II, to further improve the exploration capability, we introduced chaotic initialization in two different ways: i.e., with and without Lyapunov-guided chaotic series (see Table 8, (see Table 9)) resulting in CLQBDE and CQBDE respectively. In the latter case (i.e., without the Lyapunov exponent), we noticed a reduction in mean AUC in all the datasets than their corresponding non-chaotic counterparts (i.e. QBDE variants). This is likely due to the fact that we do not operate in a true chaotic regime.

To balance both the exploration and exploitation capability of an algorithm, we introduced Lyapunov-guided chaotic series (see Table  9)in the initialization phase. It turned out that these variants achieved efficient solutions meaning higher mean AUC with decreased mean cardinality. The chaotic numbers introduced in the first step (of CQBDE variants) are not generated in a truly chaotic regime. However, after introducing the Lyapunov exponent-guided chaotic numbers (i.e. CLQBDE variants), we discarded the first 5000 time steps, ensuring the chaotic series is in the true chaotic regime. This transition helped the algorithm improve exploration and exploitation capability as evidenced by the results in Table 9). This is spectacularly observed in the case of both high-dimensional datasets, i.e., OVA_Omentum and OVA_Uterus. However, the true chaotic regime alone did not suffice in the other two datasets.

Further, to create a level-playing field to the algorithm we would like to compare with (i.e. variants of CQIEA [35] and CLQIEA), we demonstrated the effectiveness of the threshold trick and the influence of Lyapunov exponent-guided chaotic numbers by invoking them in the variants of CQIEA [35] and CLQIEA in Table 10. We noticed that the original CQIEA [35], i.e., without the threshold trick and the Lyapunov exponent-guided chaotic numbers, yielded unacceptably high mean cardinality, impacting its performance adversely. However, after introducing the threshold trick combined with Lyapunov exponent-guided chaotic numbers, we could effectively control the cardinality while improving the mean AUC. This behaviour is particularly evident when LASSO LR (LLR) is employed as the classifier.

4.2 Statistical Testing

A two-tailed t-test at 5% level of significance and 38 (=20+20-2) degrees of freedom is conducted on the mean AUC obtained from 20 runs across the top two algorithms (w.r.t mean AUC) in each dataset (refer to Table 11). The t-test demonstrates that in three out of four datasets (except OVA_Uterus), the top 2 best algorithms solely based on AUC, turned out to be statistically similar. The top 2 best algorithms corresponding to each dataset are as follows: (i) in the Epsilon dataset (CTQIEA+LLR and CLTQIEA+LLR), (ii) in IEEE Malware dataset (QBDE-I+LLR and CQBDE-I+LLR), (iii) in OVA_Omentum dataset (CLQBDE-II+LLR and CTQIEA+LR) and (iv) in OVA_Uterus dataset (CTQIEA+LR and CLQBDE-I+LLR).

In the case of statistical similarity, to break the tie, preference is accorded to the algorithm that selected less mean cardinality. In other words, higher preference is accorded to the algorithm that yielded a great reduction in mean cardinality and an insignificant reduction in mean AUC. After performing this type of cost-benefit analysis, the better-performing algorithm is presented in bold face in Table 11).

It is noticed that in three datasets, the top-performing algorithm yielded a higher mean AUC with almost similar mean cardinality. In these cases, preference is accorded to the one with a higher mean AUC. For example, in the IEEE Malware dataset, QBDE-I+LLR obtained a mean AUC of 0.826 with a mean cardinality of 168.5. However, the next-best algorithm (CQBDE-I+LLR) obtained a mean AUC of 0.818 (which is $<$ 0.008 than that obtained by QBDE-I+LLR) with a mean cardinality of 162.55. Here, the reduction in cardinality is very minimal ($<$ 6 features). This makes the QBDE-I+LLR win over the CQBDE-I algorithm. The same cost-benefit analysis is followed for the other datasets as well, resulting in the Epsilon and OVA_Omentum datasets; the QBDE-I + LLR and CLQBDE-II+LLR are the winners, respectively. However, in the OVA_Uterus dataset, CLQBDE-I + LLR turned out to be statistically significant when compared to the CTQIEA-I + LR. Accordingly, CLQBDE-I + LLR is presented in boldface in Table 11

It is important to note that when algorithm A and algorithm B are statistically similar, the winner is chosen based on the cost-benefit analysis. However, if statistical significance is observed, the preference is accorded to the algorithm, which is the statistically significant algorithm.

In summary, the insights derived from the current study are as follows:

• In the Epsilon dataset, CTQIEA + LLR emerged as the winner due to the significant impact of the threshold trick and chaos.

• In the IEEE Malware dataset, QBDE-I+LLR was the best algorithm, primarily due to the effectiveness of the threshold trick.

• In the OVA_Omentum dataset, CLQBDE-II + LLR outperformed others due to the introduction of Lyapunov-based chaos and the threshold trick.

• In the OVA_Uterus dataset, CLQBDE-I + LLR was the winner largely due to the guidance by Lyapunov.

• The CLQBDE variants with LLR as a classifier turned out to be the best wrappers no matter a quantum gate is adopted or not in both high-dimensional datasets.

5 Conclusions

This study proposes CLQBDE-I, where a Lyapunov exponent-guided chaotic map-based initialization method is incorporated into the quantum-inspired BDE algorithm for FSS. Our results show its superiority on high-dimensional datasets and its competitive performance on Epsilon and IEEE Malware datasets compared to the alternative algorithms. CLQBDE-I outperformed not only our baseline, namely, QBDE but also other baselines CLTQIEA-I and CTQIEA-I in all but the IEEE Malware dataset in the latter algorithm. Overall, working in conjunction with QBDE variants, integrating Lyapunov exponent-guided chaotic dynamics into them yielded better solutions than simple, chaotic-based initialization methods.

Future research directions include extending this approach to multi-objective environments, designing hybrid chaotic mapping techniques, and exploring chaotic quantum hybrid EAs. This methodology also holds promise for large-scale clustering and feature selection applications across various domains like finance and economics.

Appendix A Details of the proposed approach

Figure 3: Histograms

Table 12: Summary of parallel evolutionary algorithms in different environments

Authors	Algorithm	Environment	Solved Problem
Zhou [17]	DE	Spark	Discussed Pros and cons of various parallel approaches
Teijeiro et al. [18]	DE	Spark + AWS	Benchmark functions
Cho et al. [19]	DE	Spark	Clustering
Chen et al. [20]	Modified DE	SPMD	Cluster Optimization
Al-Sawwa and Ludwig [21]	DE	Spark	DE based classifier
Adhianto et al. [22]	DE	OpenMP	Optical Network problem
Deng et al. [23]	DE	Spark	Benchmark functions
He et al. [24]	Five variants of DE	Spark + Cloud	Ring topology model applied to benchmark functions
Veronse & Krohling [25]	DE	CUDA	Large scale optimization
Wong et al. [26]	Self-Adaptive DE	CUDA	Benchmark functions
Cao et al. [27]	DPCCMOEA	MPI	Large scale optimization
Ge et al. [28]	DDE-AMS	MPI	Large scale optimization
Falco et al. [29]	DE	MPI	Resource allocation
Glotic et al. [30]	PSADE	MATLAB	Hydro Scheduling algorithm
Daoudi et al. [31]	DE	Hadoop	Clustering
Thomert et al. [32]	NSDE-II	OpenMP	Cloud work placement
Kromer et al. [33]	DE	Unified Parallel C	Large scale optimization
Qian et al. [34]	MPFPSP	Multithreading	Flow scheduling problem
Vivek et al. [12]	PB-TADE, PB-DETA,PB-DE	Spark	FSS
Current study	CLTQBDE-I, CLTQIEA	Spark	FSS

Algorithm 1 Driver Algorithm

1:Input: $ps$, $lps$, $X$, $Xtrain$, $Xtest$, $mMig$, $mGen$

2:Output: $P$: population evolved after $mMig$ migrations

3:$i\leftarrow 0$

4:$Qt\leftarrow$ Chaotically Initialize the quantum matrix

5:$X\leftarrow$ getBinarySolution($Qt$)

6:$P\leftarrow[\phi,\phi,\ldots,\phi]_{ps\times 4}$

7:for $i=0,1,\ldots,ps$ do

8:     $sol\leftarrow[i,X[i],\phi,0.0]$ $\triangleright$ Encoding Scheme (see Table 3)

9:     $P\leftarrow P\cup sol$

10:end for

11:$(Xtrain_{1},Xtrain_{2},\ldots,Xtrain_{k})\leftarrow Xtrain$ $\triangleright$ Divide data into $k$ islands

12:while $i=0,1,\ldots,mMig$ do

13:     Call IslandMapper($P$) $\triangleright$ Migration Rule

14:     $R\leftarrow$ Collect $\{p_{r_{1}},p_{r_{2}},\ldots,p_{r_{k}}\}$ from $k$ islands

15:     $R\leftarrow$ SortBasedOnFitness($R$)

16:     $P\leftarrow\{R:R\text{ until }|R|<ps\}$

17:end while

18:for each sol in $P_{r}$ do

19:     $model\leftarrow$ sol[2] $\triangleright$ Collect coefficients for each solution

20:     $score\leftarrow$ testModel($model,Xtest$) $\triangleright$ Evaluate AUC on test dataset

21:     $cardinality\leftarrow$ sum(sol[1])

22:     $testFitness\leftarrow(score\times(1-(cardinality/m)))$

23:end for

24:Return $P$, $testFitness$

Algorithm 2 Worker algorithm

1:Input: $lps$, $mGen$, $F$, $CR$

2:Output: $localP$: population evolved after $mGen$ maximum iterations

3:$k\leftarrow 0$

4:localP $\leftarrow$ Randomly pick lps number of solutions from P $\triangleright$ local population

5:for $i=0,1,\ldots,\text{lps}$ do

6:     bVC $\leftarrow$ localP[i].BinaryEncodedVector

7:     auc $\leftarrow$ LLR (bVC, $islandData$)

8:     localP[i][3] $\leftarrow$ updateAUC

9:end for

10:for $k=0,1,\ldots,\text{mGen}$ do

11:     mV $\leftarrow\{\phi\}$

12:     for $i=0,1,\ldots,\text{lps}$ do

13:         loPVec $\leftarrow$ localP[i].QuantumMatrix

14:         mVec $\leftarrow$ Mutation(loPVec)

15:         mV $\cup\{mVec\}$

16:     end for

17:     locOfP $\leftarrow\{\phi\}$

18:     for $i=0,1,\ldots,\text{lps}$ do

19:         mV $\leftarrow$ mutatedP[i].QuantumMatrix

20:         TV $\leftarrow$ Invoke Crossover on mV

21:         locOfP $\cup\{\text{TV}\}$

22:     end for

23:     Invoke train-and-update phase on locOfP

24:     newP $\leftarrow\{\phi\}$

25:     for $i=0,1,\ldots,\text{lps}$ do

26:         nS $\leftarrow$ Invoke Selection (locOfP[i], P[i])

27:         newP $\cup\{\text{nS}\}$

28:     end for

29:     localP $\leftarrow$ newP

30:     $k\leftarrow k+1$

31:end for

32:return localP