Evolutionary Multitasking-Based Multiobjective Optimization Algorithm for Channel Selection in Hybrid Brain Computer Interfacing Systems

Hybrid-modality brain-computer Interfaces (BCIs), which combine motor imagery (MI) bio-signals and steady-state visual evoked potentials (SSVEPs), has attracted wide attention in the research field of neural engineering. The number of channels should be as small as possible for real-life applications. However, most of recent works about channel selection only focus on either the performance of classification task or the effectiveness of device control. Few works conduct channel selection for MI and SSVEP classification tasks simultaneously. In this paper, a multitasking-based multiobjective evolutionary algorithm (EMMOA) was proposed to select appropriate channels for these two classification tasks at the same time. Moreover, a two-stage framework was introduced to balance the number of selected channels and the classification accuracy in the proposed algorithm. The experimental results verified the feasibility of multiobjective optimization methodology for channel selection of hybrid BCI tasks.


INTRODUCTION
Brain-computer interface (BCI) develops a communication means between brains and external devices (Alcaide-Aguirre and Huggins, 2014). This technology is used for developing assistive devices for helping paralyzed patients or control gaming. Due to its convenience, electroencephalography (EEG) is widely adopted in non-invasive BCIs for multi-channel signal acquisition (Kevric and Subasi, 2017). In general, more channels would provide more information that can improve the performance of BCI classification. However, channel selection is essential to filter redundant information and simplify experimental manipulation. For example, MI tasks may only need a number of channels as few as 3-5 without lowering the classification accuracy according to temporal-spatio analysis (Pfurtscheller and da Silva, 1999). Furthermore, using more channels is inconvenience for clinical operation and signal processing. Therefore, how to make a compromise between channel number and classification accuracy is a worth-probing question, especially for hybrid BCIs in real-world applications.
In previous works, several effective methods have been proposed for channel selection, such as filter-based methods (Yang et al., 2016), wrapper-based methods (Qiu et al., 2016), correlationbased methods (Jin et al., 2019), machine learning methods (Su et al., 2019), Heuristic searching methods (Sun et al., 2021), and so on. As mentioned before, channel selection needs to make a compromise between two objectives: the number of selected channels and classification accuracy. This perfectly coincides with the goal of multiobjective evolutionary algorithms (MOEAs), which is to make a compromise among multiple optimization objectives. In recent years, multiobjective evolutionary algorithms (MOEAs) have demonstrated their effectiveness in channel selection. For example, the classic Non-dominated Sorting Genetic Algorithm-II (NSGA-II), Multiobjective Evolutionary Algorithm Based on Decomposition (MOEA/D), and Multiobjective Particle Swarm Optimization (MOPSO) have been successfully applied for channel selection in the task of single-modality -based BCIs Moubayed et al., 2010;Kee et al., 2015). Current studies show that hybrid BCI, which combines motor imagery (MI) signals and steady-state visual evoked potentials (SSVEPs), can achieve certain goals more efficiently than singlemodality -based BCI systems (Long et al., 2012;Ko et al., 2014). However, relevant works only focus on either the MI classification task (Jin et al., 2020;Sun et al., 2021) or the SSVEP classification task (Zhang et al., 2017;Ravi et al., 2019). Few works conduct channel selection for hybrid BCI tasks.
In this paper, a multitasking-based multiobjective evolutionary algorithm (EMMOA) is proposed to perform channel selection for MI and SSVEP tasks at the same time. In EMMOA, the problems of channel selection for both MI and SSVEP tasks can be regarded as multi-objective optimization problems (MOPs). In order to make a balance between classification accuracy and channel number, the problem of channel selection for MI task can be formulated as a twoobjective optimization problem, which contains two objectives: the classification accuracy for MI task and the number of selected channels. Evolutionary multitasking mechanism, which is inspired by bio-cultural models of multifactorial inheritance (Gupta et al., 2015(Gupta et al., , 2016, provides the investigation of solving multiple tasks via a single population concurrently. In the evolutionary multitasking mechanism, different tasks will experience information transfer during the evolution process since they use the same population. Therefore, if multiple tasks are related then the searching process of solving one task may offer help in solving the other tasks (Gong et al., 2019;Bai et al., 2021). As the evolutionary multitasking mechanism is proved to be an efficient way to optimize multiple tasks, it has been introduced in EMMOA to optimize MI and SSVEP tasks simultaneously in this paper. Furthermore, EMMOA adopts a two-stage framework to improve searching efficiency. The first stage is based on an evolutionary multitasking mechanism and aims to obtain the Pareto-optimal solutions (PS) for MI and SSVEP tasks by one single population. The second stage is local searching. The second stage constructs a three-objective optimization problem, which used classification accuracy for MI task, classification accuracy for SSVEP task, and the number of selected channels as the optimization objectives. The constructed three-objective problem is optimized according to decision variable analysis based on the results obtained by the first stage.
The rest of this paper is organized as follows: Section 2 elaborates the data set, pre-processing and optimization methods used in this paper. Section 3 presents experimental and result analysis, followed by further discussion of this paper.
In this experiment, 7 healthy volunteers, aged from 21 to 30, were selected for data acquisition. Subjects had no prior experience with the operation of hybrid BCI. They all gave informed consent approved by the Ethics Committee.

Feature Selection and Classification
In our work, MI tasks including left-and right-hand imagery are performed for pattern recognition. EEG signals are band pass filtered from 5 to 30 Hz firstly. Then the features are extracted by a well-known common spatial pattern (CSP) algorithm. The CSP method is useful for discriminating two populations of EEG dataset. The detailed methodology can be reviewed in Muller-Gerking et al. (2004). And a radial basis function kernel support vector machine (RBF-SVM) classifier is used for training these feature vectors.
From Lin et al. (2006), we utilize the canonical correlation analysis (CCA) method for SSVEP detection. CCA reflected the correlation relationship between EEG response signals and classical Fourier series at the stimulus frequency and its harmonics. It is widely used for spectral analysis in the biosignal process.

Multi-Objective Channel Selection Problem Formulation
Suppose the number of total channels is K, then a solution for the channel selection problem can be defined as a K-dimensional vector x, as shown in Equation (1). x i (1 ≤ i ≤ K) is called a decision variable. x i = 0 means the corresponding channel is not selected. Otherwise, the corresponding channel is selected.
The purpose of channel selection problems is to select a smaller set of channels with as little sacrifice as possible in classification accuracy. For the MI classification task, the first objective function is the MI classification accuracy rate (MAR). The second objective function is defined as NC = K − C, where C is the number of selected channels (i.e., the number of elements whose value is 1 in solution x). Therefore, the ideal situation for an optimization algorithm is obtaining an optimal solution x * that has the maximum values for both MAR and NC. However, it is impossible to get the abovementioned x * , since the confliction between MAR and NC. Thus, evolutionary optimization algorithms will obtain a Pareto Set (PS), which contains a set of Pareto-optimal solutions (nondominated solutions) (Wang et al., 2018). One solution x is called a Pareto-optimal solution, if and only if there not exist any other solutions that are better than x for all objective functions. The Pareto optimal solutions are incomparable, since solution x may be better than solution y for one objective function but worse than y for the other one. Similarly, the two objectives for SSVEP classification task are the SSVEP classification accuracy rate (SAR) and NC. In the proposed algorithm, the second stage aims to improve the optimization performance for both MI and SSVEP tasks by local searching, therefore a three-objective optimization problem, which uses MAR, SAR, and NC as the objective functions, is adopted in this stage.

Evolutionary Multitasking-Based Multiobjective Optimization Algorithm (EMMOA)
The framework of EMMOA is given in Figure 2. As shown in Figure 2, EMMOA adopts a two-stage framework. In EMMOA, the first stage is designed according to evolutionary multitasking mechanism and uses one single population to optimize two tasks (MI task and SSVEP task) simultaneously. In this case, the two tasks will experience information transfer during the evolution process since they use the same population. The information transfer can share underlying similarities between the two tasks thereby facilitating improved optimization performance for both MI and SSVEP tasks. Based on the output of the first stage, the second stage aims to obtain the final output of the algorithm. With the help of decision variables analysis, the local searching strategy in the second stage makes a better compromise among MI classification accuracy, SSVEP classification accuracy, and the number of selected channels. More specifically, in the second stage, a decision variable analysis operator will be carried out on the Pareto-optimal sets obtained by the first stage and aims to divide the decision variables into different groups. The local searching operator will determine the searching directions of individuals according to the type of each decision variable and help the algorithm search more efficiently. The detailed description of the two stages is given below.
In the first stage, suppose the population size is N and the total number of channels is K, then the individual population can be initialized by generating N individuals (solutions) and each individual contains K elements. In this case, the individual population can be initialized as a N × K matrix. Each element in the individual population is uniformly generated from [0, 1]. If the value of an element is larger than 0.5, then the element will take the value of 1. Otherwise, the element will be set to 0. By associating each individual with a task label by the task assignment operator. In the initialization step, the task assignment operator is carried out by giving each individual the label, which is generated from {1, 2} randomly. In the evolution process, the task labels of individuals are inherited from their parents. Each individual will optimize the task that is specified by its task label. For a task label, 1 and 2 represent the corresponding individual is assigned to optimize MI task and SSVEP task, respectively. The individual population is updated by using the population evolution step, which includes the tournament selection operator , partial-mapped crossover operator (Ismkhan and Zamanifar, 2015), and flip-bit mutation operator (Chicano et al., 2015). The tournament selection operator, partial-mapped crossover operator, and simple mutation operator are classic operators, which are widely adopted in evolutionary algorithms to generate offspring populations (lines 1-21 in Algorithm 1). The task labels of offspring individuals are inherited from their parents in the evolution process (line 7 and line 19 in Algorithm 1). After updating the individual population (lines 22-23 in Algorithm 1), the Pareto-optimal sets for MI and SSVEP tasks will be updated accordingly (line 24 in Algorithm 1). Therefore, the first stage will output two Pareto-optimal sets (namely PS_MI and PS_SSVEP) for MI and SSVEP tasks separately. PS_MI contains a set of non-dominated solutions which aim to optimize MAR and NC, while PS_SSVEP contains a set of non-dominated solutions which aim to optimize SAR and NC. In this paper, the maximum size of both PS_MI and PS_SSVEP is set to 100. Take PS_MI for an example, if the size of PS_MI exceeds 100, then the non-dominated solutions in PS_MI will be decreased according to Crowding Distance (Deb et al., 2002). Crowding Distance, which is proposed to estimate the density of solutions, is widely adopted to maintain population diversity in EAs. The solution with a larger Crowding Distance value can be considered to make great contributions for maintaining population diversity. Therefore, if the number of the solutions in PS_MI exceeds 100, then the members with larger Crowding Distance values will be remain.
The second stage is implemented based on the output of the first stage, i.e., the current Pareto-optimal sets for MI and SSVEP tasks (PS_MI and PS_SSVEP). The purpose of the second stage is to obtain a set of Pareto-optimal solutions that make a compromise among MI classification accuracy, SSVEP classification accuracy, and the number of selected channels. Therefore, in this stage, the optimization problem becomes a maximum three-objective optimization with MAR, SAR, and NC as its objective functions. In this case, the Pareto-optimal set outputted by the whole algorithm consists of a set of nondominated solutions which aim to make a compromise among MAR, SAR, and NC. To improve the searching efficiency, a local searching operator, which is based on a decision variable analysis strategy, is introduced in this stage. The main idea of the proposed decision variable analysis operator is to divide the decision variables into three groups: add_group, delete_group, and invalid_group. The division of the first 11 and the last 4 decision variables are based on PS_SSVEP and PS_MI, respectively. Take the jth (1 ≤ j ≤ 11) decision variable for an example. At first, initialize flag to 0. For a non-dominated solution (x i ) in PS_SSVEP, if the classification accuracy is raised with x i (j) = 0, then one can consider the performance of x i will be improved without selecting the jth channel. In this case, flag = flag − 1. If x i (j) = 1 gets a better classification accuracy, then flag = flag + 1. If the classification accuracy is not affected by x i (j), then the value of flag will remain unchanged. All the members in PS_SSVEP will be executed to get the final value of flag.The jth decision variable will be divided into different groups according to the value of flag. The detailed procedure of the proposed decision variable analysis operator is given in Algorithm 2.
After decision variable analysis, the local searching operator is implemented on PS_MI and PS_SSVEP. The detailed procedure of the proposed local searching operator is shown in Algorithm 3. For the Pareto-optimal solutions in PS_MI, the searching directions are obtained according to the types Frontiers in Neuroscience | www.frontiersin.org

of the last four decision variables (lines 2-11 in Algorithm 3).
For the Pareto-optimal solutions in PS_SSVEP, the searching directions are obtained according to the types of the first eleven decision variables (lines 12-21 in Algorithm 3). The current Pareto-optimal set of EMMOA will be acquired by selecting the non-dominated solutions in the original Pareto-optimal set and newPOP obtained by the local searching operator.

Experimental Setup
This section contained two experiments. The first one aimed to demonstrate whether the proposed channel selection algorithm (EMMOA) could obtain a smaller set of channels with as little sacrifice as possible in MI and SSVEP classification accuracy. The second one aimed to evaluate the effectiveness of EMMOA by comparing with three widely used multiobjective optimization algorithms, including NSGA-II (Deb et al., 2002), MOEA/D (Zhang and Li, 2007), and MOPSO (Coello et al., 2004). The number of total function evaluations was 10000 in a single run for all 4 algorithms. The detailed parameter settings of algorithms were given in Table 1.

Results and Analysis
The purpose of the proposed EMMOA is to find the optimum number of channels that would result in a promising performance for both MI and SSVEP tasks. Take subject 3 as an example, Figure 3 illustrated the classification accuracies for MI and SSVEP tasks with different numbers of selected channels. It can be observed from Figure 3, the classification accuracy for MI task tends to become higher at first and then become lower as the number of selected channels increases. The classification accuracy for SSVEP task tends to be stable when the number of  Figure 4 demonstrated that the larger number of selected channels cannot lead to better performance for MI and SSVEP tasks. This may indicate the necessity of selecting a promising subset of channels. Table 2 demonstrated the classification accuracy for MI and SSVEP tasks by using all channels and the reduced set of EEG channels selected by the proposed EMMOA. As Table 2 shows, EMMOA can find a smaller set of channels without depredating the classification performance for all 7 subjects. For example, the classification accuracies for MI and SSVEP tasks using 5 channels Algorithm 3 Detailed procedure of local searching operator. Input: add_group, delete_group, and invalid_group n1 (number of the solutions in PS_MI), n2 (number of the solutions in PS_SSVEP) Output: newPOP 1: newPOP=∅; 2: for i = 1 TO n1 do 3: x is the i th solution in PS_MI, x1 = x;  selected by EMMOA are better than those using all channels. For subject 6, the classification performance using all channels is same to that of using 5 channels. However, a smaller number of channels is certainly more convenient for real-life applications. Table 3 gives the average classification accuracies for all subjects by using all channels and the Pareto-optimal set obtained by EMMOA in a single run. It can be observed from Table 3, the classification performance by using 6, 7, and 8 channels selected by EMMOA is better than using all channels. Therefore, the  experimental results demonstrated the effectiveness of using the proposed channel selection algorithm. For further comparison of the proposed EMMOA with other classical multi-objective evolutionary algorithms, three algorithms, including NSGA-II, MOEA/D and MOPSO, are adopted in this section. All the comparative algorithms have been running 30 times to get the statistical results, which are shown in Table 4. In Table 4, the numbers in brackets indicate the rank of the corresponding algorithm in terms of average MI accuracy, average SSVEP accuracy, and average number of selected channels, respectively. A smaller rank value of an algorithm means a better performance achieved by the algorithm. In the last column of Table 4, EMMOA obtained the smallest (best) average rank. In other words, the proposed channel selection algorithm achieved the best performance in terms of average rank. Table 5 gives the maximum average classification accuracies for MI and SSVEP tasks in the PS obtained by all the 4 algorithms with different numbers of selected channels. The numbers in boldface give the best accuracy rates in different cases. It can be observed from Table 5, EMMOA obtained the best classification accuracies for all 3 cases.
The results in Tables 4, 5 demonstrated the effectiveness of the proposed algorithm. This may be because the multitasking mechanism adopted in the first stage of EMMOA help the algorithm find more promising PS for each task efficiently. On the other hand, with the help of decision variables analysis, the local searching strategy in the second stage makes a better compromise among MI classification accuracy, SSVEP classification accuracy, and the number of selected channels. Figures 5, 6 illustrated the convergence of four algorithms for MI and SSVEP tasks, respectively. In Figures 5, 6, the x-axis represented the number of function evaluations and the y-axis displayed the maximum classification accuracy averaged over all the 7 subjects. As Figure 5 showed, the proposed EMMOA obtained the second-best and the best convergence speed for MI and SSVEP tasks among four algorithms, respectively. This is because the decision variable analysis strategy adopted in EMMOA helped the algorithm searching more efficiently.

DISCUSSION
It can be observed from Figure 4, the best classification accuracy can be achieved when the number of selected channels is not very large. This observation can be   Table 3. The classification accuracies for both MI and SSVEP tasks are higher than 0.65 when the number of selected channels is larger than 3. As Figure 4 showed, the improvement in classification accuracy is not obvious by adopting 4 channels when compared with using more than 4 channels. Especially for MI task, the classification accuracy by using 4 channels is even better than using all 15 channels. The phenomenon is in agreement with the conclusion in Moubayed et al. (2010), which showed that some channels may bring artifacts and then lead to a degeneration in classification accuracy.
In the proposed algorithm, the decision variable analysis operator played a vital role in the second stage of EMMOA, since it provided searching directions for the next local searching operator. Table 6 showed the results outputted by the decision variable analysis operator in the middle and last stages of the optimization process by using EMMOA. add_group, delete_group and invalid_group contains a set of channels that are considered to be profitable, disadvantageous, and invalid for the corresponding classification task. As Table 6 shows, most results obtained by the decision variable analysis operator in the middle stage overlap with those obtained in the last stage. However, the results in different stages are not exactly the same. This phenomenon may indicate the uncertainty of the decision variable analysis operator to some extent. As described in section 2.4, the performance of the second stage of EMMOA depends on the results of the proposed decision variable analysis operator. Therefore, the second stage may not obtain a promising output if the decision variable analysis operator cannot get an appropriate division of the decision variables, especially in the early stage of the evolution process. How to improve the stability of the proposed decision variable analysis operator will be the next work of this paper.
As shown in Figure 7, the distribution diagrams were illustrated in conditions of different numbers of selected channels. Evidently, the electrodes were first eliminated from motor cortex when the number was sufficient for MI and SSVEP recognition. And the quantity of channels were close between central area and occipital area along with the decline of total number. Meanwhile, the priority of selecting was unclear in these two regions because of the alternate quantitative superiority. However, the feature information needed to be reserved those extracted from left-and right-motor cortex for spatial filtering. Thus, the limitation of channel number in the occipital area might be firstly considered as the kernel factor of EMMOA.

CONCLUSION
In our study, a multitasking-based multiobjective evolutionary algorithm (EMMOA) was proposed to select appropriate channels for these two classification tasks at the same time. Moreover, a two-stage framework was introduced to balance the number of selected channels and the classification accuracy in the proposed algorithm. The experimental results verified the feasibility of multiobjective optimization methodology for channel selection of hybrid BCI tasks.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of Tongji University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
TL designed the methodology of data process and performed the data analysis. TL and ZX organized the data and wrote the manuscript. LC and GT reviewed and edited the manuscript. All authors read and approved the submitted manuscript.

FUNDING
This work was supported by the National Natural Science Foundation of China (grant nos. 61806122 and 62102242).