Feature selection for high-dimensional data based on scaled cross operator threshold filtering specific memory algorithm

Wavelet threshold filtering algorithm is a de-noising method that is simple to implement, small in arithmetic and good in filtering effect. The basic principle is: the noisy image is decomposed into wavelet coefficients of different amplitudes after wavelet transform, the energy of the image signal is concentrated in a few wavelet coefficients, which makes the wavelet coefficients of large amplitude, while the noise energy corresponds to the wavelet coefficients of smaller amplitude, according to the different sizes of the wavelet coefficients corresponding to them, a threshold is set to distinguish these coefficients, the wavelet coefficients of amplitude lower than the threshold is rejected, and wavelet coefficients of amplitude greater than the threshold is retained or shrinkage is done accordingly. The wavelet coefficients whose magnitude is lower than the threshold are eliminated, the wavelet coefficients whose magnitude is greater than the threshold are retained or shrinkage processed accordingly, and finally the remaining wavelet coefficients are utilized for image reconstruction, which realizes the denoising of the image and retains the detailed information of the image well.

Using the threshold method of filtering, there are two key points: the first is the selection of the threshold parameters; the second is the selection of the threshold function.

The wavelet threshold filtering algorithm relies on the threshold selection, and its size has a direct impact on image denoising. If the threshold is too large, it will make the useful high-frequency information (such as edges, details of texture information) is lost, resulting in image distortion, blur; if the threshold is too small, it will make too much noise residual and lead to poor filtering effect.

Currently commonly used image denoising thresholds are VisuShrink thresholds.VisuShrink thresholds are also unified thresholds or universal thresholds, and their threshold λ is expressed as follows: (5)λ=σ2lnN$$\lambda = \sigma \sqrt {2\ln N}$$

In equation (5), σ is the standardized variance of the noise and N is the size of the noise-containing image or the length of the signal. There is another condition for the use of this threshold, which is that the noise variance must be evaluated first. In practice, the standard variance σ^$$\hat \sigma$$ of the noise is unknown, and since the noise is mainly concentrated in the smallest scale to get HH₁ sub-band, we can utilize the median estimation method, i.e., wavelet coefficients of the HH₁ sub-band to evaluate the noise variance. (6)σ^=Median(|W[xij]|)0.6745|W[xii]|∈HH1$$\hat \sigma = \frac{{Median\left( {\left| {W\left[ {{x_{ij}}} \right]} \right|} \right)}}{{0.6745}}\left| {W\left[ {{x_{ii}}} \right]} \right| \in H{H_1}$$

In Eq. (6), σ^$$\hat \sigma$$ is the estimate of the standardized variance of the noise, and Median(|W[xij]|)$$Median\left( {\left| {W\left[ {{x_{ij}}} \right]} \right|} \right)$$ denotes the median of the absolute values of the wavelet coefficients of the first layer of high-frequency subbands HH₁ after wavelet decomposition.

Threshold function is also called the threshold contraction strategy, threshold function by thresholding the noise-containing high-frequency wavelet coefficients, it applies different processing strategies to wavelet coefficient modes greater than and less than the threshold, with the aim of maximizing the suppression of wavelet coefficients corresponding to the noise and retaining the coefficients corresponding to the image signal.

Let w_{j, k} be the unprocessed coefficients of the noisy image or signal after wavelet decomposition, w^j,k$${\hat w_{j,k}}$$ denotes the wavelet estimated coefficients output after thresholding, and λ is the image denoising threshold. The commonly used threshold function has a hard inter-value function. This function compares the absolute value of the wavelet transform coefficients w_{j, k} with the selected threshold value λ, and only retains the larger wavelet coefficients and sets the smaller wavelet coefficients to zero, with the following expression: (7)w^j,k={ 0 |wj,k|<λ wi,j |wj,k|≥λ$${\hat w_{j,k}} = \left\{ {\begin{array}{c} 0&{\left| {{w_{j,k}}} \right| < \lambda } \\ {{w_{i,j}}}&{\left| {{w_{j,k}}} \right| \geq \lambda } \end{array}} \right.$$

The interval function in the ( −λ, λ) interval is not continuous, there is a break point exists, however, in practical applications often need to carry out the derivation of the threshold function, but due to the ( −λ, λ) interval of the threshold function is not continuous all in the break point can not carry out the derivation of the operation, so there are some limitations; in addition, the threshold function is greater than the interval wavelet coefficients remain unchanged and the wavelet coefficients less than the threshold for processing, but the actual situation is greater than the threshold wavelet coefficients may also exist in the noise interference, so this threshold processing method has defects. The wavelet coefficients larger than the threshold may also have noise interference, so there are defects in this thresholding method.

Select the output measurements y_k of the data segment, the input measurements u_k and the filter parameter values to be optimized ρ After the time update and state update of the improved adaptive threshold filter, the a posteriori estimation of the system state is obtained x^(k)$$\hat x\left( k \right)$$; compute the objective function needed for optimization, and use the genetic algorithm to perform optimal retention, replication, crossover, and mutation operations in the m solution space to obtain the new corrected population, which is reintroduced into the AKF for parameter estimation, and repeat the above steps for the measurement information until the objective function meets the design requirements or the population update of the genetic algorithm ends.

Genetic algorithm is a kind of global search and optimization algorithm independent of the domain of the problem under study, with good parallel computing performance and global search capability. The basic process of genetic algorithm is shown in Fig. 2, with reference to the genetic phenomena of biological natural selection, gene crossover and mutation, the solution object is coded and expressed, and then a reasonable fitness function is set to select the good population individuals to be copied, and the group with greater fitness is inherited to the next generation, and with reference to the mutation phenomenon, the mutation operation is performed randomly on individual genes in order to expand the search space. In this way, continuous reproduction is carried out until the individual with optimal fitness is obtained or the set loop termination condition is satisfied.

Figure 2.

Basic flow of genetic algorithm

The computational steps generally include: 1)

Initialization: the initialization operation generally consists of solving for the binary, floating-point, or symbolic encoding of the parameters, setting the population size, and generating the initial population and evolutionary generations.

In genetic algorithms, the crossover operator is the main feature that distinguishes them from other evolutionary algorithms. By simulating the biological hybrid reproduction process to perform gene crossover operations on individuals, it increases the diversity of individuals, enlarges the range of optimality search, and plays the role of reaching the global optimal solution in the solution space.

The specific operation is as follows: in the selected population, two individuals are selected according to the design pattern, and part of the genes of both are exchanged in a designed way according to the crossover probability P_c, so as to reproduce to get 2 new individuals. When the probability of crossover P_c is large, the crossover operator has a strong ability to reproduce new individuals, and thus has a better ability to explore the global optimal solution in the solution space; comparatively, when the probability of crossover P_c is small, the probability of the good individuals of the population to be destroyed is smaller, and the ability to find the optimal solution in the localization of the current solution is better.

In this paper, floating point coding is used to express the gene values in each individual as a floating point number in the corresponding range, and the coding length of the individual is equal to the number of variables in the coding length. The single-point crossover, two-point crossover and uniform crossover applicable to binary can also be applied to floating-point coding, while arithmetic crossover operator, scaling crossover operator, etc. can also be used based on floating-point coding. Considering the large length of individuals, this paper uses the scaling crossover operator. After the initial population selection parameter values are randomly dispersed in the solution space, the general population size is large, after mixing the weightless selection operator, the existing size can be obtained less than or equal to the initial population size, followed by carrying out the selection operator operation in the existing population.

The scaling crossover probability is set up in such a way that after sorting the fitness size of the existing population individuals, two individuals are selected to generate new individuals by non-uniform linear combination with the proportion of their fitness as the benchmark. Assuming an arithmetic crossover between Individual A and Individual B, the probability of their selection and the two new individuals produced are: (8)Pc=fAfA+fB$${P_c} = \frac{{{f_A}}}{{{f_A} + {f_B}}}$$ (9)A′ = PcB+(1−Pc)A B′ = PcA+(1−Pc)B$$\begin{array}{rcl} A^\prime &=& {P_c}B + \left( {1 - {P_c}} \right)A \\ B^\prime &=& {P_c}A + \left( {1 - {P_c}} \right)B \\ \end{array}$$

The following describes the selection strategy for the two crossover individuals, the first choice is to find the largest and smallest fitness values based on the fitness of the existing population after sorting, which account for fp_max and fp_min of the sum, respectively, and the probability of the two being selected in the crossover is 1 and 0. The probability of the remaining individuals to be selected in the crossover operator is then shown in Equation (10) (10)Pi=fpmin+(fpmax−fpmin)i−1n′−1$${P_i} = f{p_{\min }} + \left( {f{p_{\max }} - f{p_{\min }}} \right)\frac{{i - 1}}{{n^\prime - 1}}$$

In the form of a roulette wheel as shown in Fig. 3, after selecting two individuals, the scaling crossover operator is operated according to Eq. (9), and at the end of the crossover step, the new individual is added to the population, so that the number of existing populations is n″.

Figure 3.

Roulette selects cross individuals

Since Pk/k−1$${P_{k/k - 1}}$$, K_k and P_k contain only the pre-test data, making the filter valuation process, the new measure value correction role declines, the role of the old measure value is relatively rising, in the process of calculating the introduction of error, resulting in the valuation of the error than expected, causing the filter dispersion; therefore, the fundamental measure to prevent dispersion is to increase the weight of the new measure value in the current filtering. To this end, an adaptive threshold filtering algorithm with exponentially weighted decay memory is proposed, and its main idea is that when there is a large error in the model, the memory length of the threshold filter can be limited by the exponentially weighted decay memory factor at each moment, which plays a role in fully utilizing the new measure value to achieve the effect of suppressing the filter dispersion.

Taking k = N then, the optimal gain matrix formula for threshold filtering is: (11)KN=PN/N−1HN′T(HN′PN/N−1HN′T+RN′)−1$${K_N} = {P_{N/N - 1}}H_N^{\prime T}{\left( {H_N^\prime {P_{N/N - 1}}H_N^{\prime T} + R_N^\prime } \right)^{ - 1}}$$

In order to overcome the dispersion, it is necessary to reduce the value of K_k before the N st moment, and the relative prominence of K_N at the N rd moment, so that the corrective effect of the new measurements rises, and the effect of the obsolete measurements declines, and because of the inverse relationship between K_k and Rk′$${{R_k'}}$$, it is necessary to reduce the value of K_k before the N th moment, and it is possible to make the Rk′$${{R_k'}}$$ th moment further and further away from the N th moment; it is gradually getting bigger and bigger.

Using an exponentially weighted approach, P₀; R1′,R2′,⋯,RN′$$R_1^\prime,R_2^\prime, \cdots ,R_N^\prime$$; Q1′,Q2′,⋯,QN′$$Q_1^\prime,Q_2^\prime, \cdots ,Q_N^\prime$$ can be changed to (12)e∑i=0N=1P0;e∑i=0N=1R1′,e∑i=0N=1R2′,⋯,e∑i=0N=1λiRN′,RN′;e∑i=0N=1λiQ1′,e∑i=0N=1λiQ2′,⋯,e∑i=0N=1λiQN′,QN′$${e^{\sum\limits_{i = 0}^{N = 1} {{P_0}} }};{e^{\sum\limits_{i = 0}^{N = 1} {{{R_1'}}} }},{e^{\sum\limits_{i = 0}^{N = 1} {{{R_2'}}} }}, \cdots ,{e^{\sum\limits_{i = 0}^{N = 1} {{\lambda _i}} }}{R_N'},{R_N'};{e^{\sum\limits_{i = 0}^{N = 1} {{\lambda _i}} }}{Q'_1},{e^{\sum\limits_{i = 0}^{N = 1} {{\lambda _i}} }}{Q'_2}, \cdots ,{e^{\sum\limits_{i = 0}^{N = 1} {{\lambda _i}} }}{Q_N'},{Q_N'}$$

Replacing the data prior to moment N with Eq. (12), the derivation of the corresponding threshold filtering becomes (13)X^*k/k−1=ΦkX^k−1*$${\widehat X^*}_{k/k - 1} = {\Phi _k}\widehat X_{k - 1}^*$$ (14)X^*k=X^*k/k−1+Kk(Zk′−Hk′X^k/k−1*)$${\widehat X^*}_k = {\widehat X^*}_{k/k - 1} + {K_k}\left( {{{Z_k'}} - {{H_k'}}\widehat X_{k/k - 1}^*} \right)$$ (15)Kk*=Pk/k−1*Hk′T(Hk′Pk/k−1*Hk′T+Rk′*)−1$$K_k^* = P_{k/k - 1}^*H_k^{\prime T}{\left( {H_k^\prime P_{k/k - 1}^*H_k^{\prime T} + R_k^{\prime *}} \right)^{ - 1}}$$ (16)Pk/k−1*=eλkΦk/k−1Pk−1*Φk/k−1T+Qk−1*$$P_{k/k - 1}^* = {e^{{\lambda _k}}}{\Phi _{k/k - 1}}P_{k - 1}^*\Phi _{k/k - 1}^T + Q_{k - 1}^*$$ (17)Pk*=[I−Kk*Hk′]Pk/k−1*$$P_k^* = \left[ {I - K_k^*{{H_k'}}} \right]P_{k/k - 1}^*$$

The initial value is X0*=E[X0];P0*=var[X0]=P0$$X_0^* = E\left[ {{X_0}} \right];P_0^* = \operatorname{var} \left[ {{X_0}} \right] = {P_0}$$

In the formula, the variable with “*” is the value of the variable at the N moment of optimal filtering, compared with the standard threshold filtering, the first term in the exponentially weighted one-step prediction variance array Pk/k−1*$$P_{k/k - 1}^*$$ expression is multiplied by a factor eλk$${e^{{\lambda _k}}}$$ greater than 1, which makes Pk/k−1*$$P_{k/k - 1}^*$$ larger, and thus Kk*$$K_k^*$$ larger, as can be seen from Eq. (14), and the increase in Kk*$$K_k^*$$ increases the weight of the new measure at the moment of k of the new measurement value at the moment of 8, which effectively prevents the dispersion of the observation value.

The threshold filtering dispersion criterion is (18)e^kTe^k≤γ⋅trE(e^ke^kT)$$\widehat e_k^T{\widehat e_k} \leq \gamma \cdot trE\left( {{{\widehat e}_k}\widehat e_k^T} \right)$$

where e^k=Zk′−Hk′X^k/k−1*$${\hat e_k} = Z_k^\prime - H_k^\prime \hat X_{k/k - 1}^*$$ is the sequence of residuals, γ ≥ 1 is the reserve coefficient, and γ = 1 is taken as the strictest convergence condition, we have (19)eλk=e^kTe^k−tr[Hk′Qk−1*Hk′T+Rk′*]tr[Hk′Φk/k−1Pk−1*Φk/k−1THk′T]$${e^{{\lambda _k}}} = \frac{{\hat e_k^T{{\hat e}_k} - tr\left[ {H_k^\prime Q_{k - 1}^*H_k^{\prime T} + R_k^{\prime *}} \right]}}{{tr\left[ {H_k^\prime {{{\Phi }}_{k/k - 1}}P_{k - 1}^*{{\Phi }}_{k/k - 1}^TH_k^{\prime T}} \right]}}$$

In order to suppress dispersion and increase the weight of the new measurement value, eλk$${e^{{\lambda _k}}}$$ is generally taken to be a value greater than 1, and is calculated as 1 when eλk$${e^{{\lambda _k}}}$$ is less than 1.

In practice, since different data channels require different attenuation factors to achieve the effect of suppressing dispersion in the calculation process, multiple attenuation memory factors are added to realize the fading attenuation operation at different rates for different data channels, respectively, which can result in multiple attenuation memory factor filters with higher performance.

Equation (16) can be improved as (20)Pk/k−1=SkΦk/k−1Pk−1*Φk/k−1T+Q*k−1$${P_{k/k - 1}} = {S_k}{\Phi _{k/k - 1}}P_{k - 1}^*\Phi _{k/k - 1}^T + {Q^*}_{k - 1}$$

where S_k is an exponentially weighted matrix of multiple decaying memory factors, Sk=diag(eλ1k,eλ2k,⋯,eλ(1−1)k,eλnk)$${S_k} = diag\left( {{e^{{\lambda _{1k}}}},{e^{{\lambda _{2k}}}}, \cdots ,{e^{{\lambda _{(1 - 1)k}}}},{e^{{\lambda _{nk}}}}} \right)$$.

All the experimental results were conducted on Windows 11 with Intel Core i7-8750H with 4.10 GHz CPU and 8-GIGABYTE (GB) memory. In the experiments, two statistical conditional independence (CI) tests were used, e.g., the G² test for the discrete dataset and the Fisher z test for the continuous dataset, with significance levels of 0.02 and 0.06, respectively.In Table 1, five benchmark BNs with a sample size of 500 are depicted.These BNs have between 30 and 802 features (variables/nodes), and the benchmark BNs obtain the number of features (variables/nodes) for each of the network’s feature (variable/node) of the MB.

Table 2 shows the experimental results of the optimization method in this paper with IAMB, MMMB, STMB, and HITON-MB on five synthetic datasets with a sample size of 500. Figure 4 shows the runtime comparison of the optimization method in this paper with IAMB, MMMB, STMB, and HITON-MB on five synthetic datasets with a sample size of 500. The experiments use the F1, Precision, and Recall metrics, as well as classifiers from MATLAB’s built-in functions or tools, to evaluate the accuracy of the generated MBs with respect to the target MBs and to determine their similarity to real MBs.

Figure 4.

Running time (seconds) of the algorithm on five base BNS

Analysis of Table 2 and Fig. 4 reveals that in Child network, this paper’s optimization method is more accurate than IAMB, MMMB, STMB, and HITON-MB in terms of F1, Precision, and Recall rates, whereas the running time is faster than MMMB, STMB, and HITON-MB, but slower than IAMB. On Child10 network, the optimization method in this paper achieves higher accuracy than IAMB, MMMB, STMB in terms of F1 and Precision. In Recall, the optimization method in this paper outperforms IAMB, MMMB, and STMB.HITON-MB is more accurate than the optimization method in this paper in F1 and Recall, however, the accuracy of the optimization method in this paper is comparable to that of HITON-MB considering the accuracy of 0.06 significant level. The running time of the optimization method in this paper is higher than MMMB, STMB, HITON-MB and lower than IAMB.

The optimization method in this paper is more accurate and runs faster than MMMB, STMB, and HITON-MB on large networks like Pig and Gene networks with sample size of 500. In addition, the optimization method in this paper is more accurate than IAMB but runs slower. On dense networks like Hailfinder, the optimization method in this paper is faster than the other three algorithms but slower than IAMB. The optimization method in this paper outperforms all tested methods in terms of F1, Precision and Recall.

The running time of the optimization method in this paper is generally faster than MMMB, STMB, and HITON-MB, but slower than IAMB because the optimization method in this paper reduces the number of independence tests for the entire candidate MB using the target T, which are the conditional independence tests that are currently selected at each estimation. However, IAMB requires more data samples for each test because the number of data samples required is exponentially related to the size of the condition set. Thus, the IAMB learning technique is time-efficient but not data-intensive. When the sample size of the data set is insufficient, the IAMB algorithm is unable to reliably find MBs. due to the limited number of samples, the quality of the MBs learned by the IAMB algorithm is greatly affected in practice, and is not as effective as the optimization method in this paper overall.

Comprehensive experimental results show that compared with other existing algorithms, the optimization method in this paper generally has better performance in terms of F1, Precision, Recall, and Run-Time, and is able to find the MB of the target feature T more reliably.

In terms of parameter settings, this paper uses the default parameters of weka for the optimization methods in this paper. Table 3 shows the parameter settings of the other four algorithms. All algorithms use the performance of the classifier as the evaluation criterion for the feature subset. To be fair, KNN (k=1) is used as the classifier in this paper, and the maximum number of iterations and population size of all algorithms are set to 60 and 350, respectively. It should be noted that all the algorithms have both single-objective and multi-objective optimization algorithms, but the classification error rate is a common concern for all feature selection algorithms, so this paper selects the individual with the lowest training error rate from the final solution set of each algorithm at the last iteration for comparison.

After setting up the experimental parameters, this paper compares the classification error rates of the five algorithms in the datasets. In a total of 20 comparisons in 5 datasets, the optimization method in this paper won 16 times, drew 1 time, and lost 3 times. In terms of classification error rate, the optimization method of this paper ranks first in 4 out of 5 datasets. Therefore, the experimental results verify that the optimization method presented in this paper is effective in dealing with high-dimensional datasets. The advantage of this paper’s optimization method over other more advanced feature selection algorithms is mainly the combination of two search operators with causal relationships using the scaling crossover operator.

To further validate the effectiveness of the optimization method in this paper, Table 4 gives the results of comparing the nonparametric effect sizes of the optimization method in this paper with those of other algorithms using Cliff’s delta, which denotes the degree of overlap between the two sets of samples, ranging from -1 to 1. In general, when |δ| is less than 0.146, it means that the two groups of samples are close to each other. When |δ| is between 0.146 and 0.32, it indicates that the difference between the two sets of samples is small. When |δ| is between 0.32 and 0.473, it indicates that there is a moderate difference between these two groups of samples. When |δ| is greater than 0.473, it indicates that there is a large difference between these two groups of samples. From the results listed in Table 4, it can be seen that the optimization method in this paper has moderate or large differences with the more advanced algorithms on most of the data sets.

In order to further validate the effectiveness of the scaling crossover operator proposed in the optimization method of this paper, three scaling crossover operators of the same class are used in this part for performance comparison, namely arithmetic crossover operator, selection operator and mutation operator. Figure 5 shows the performance of the scaling crossover operator, arithmetic crossover operator, selection operator and mutation operator on “Child10”, “Pig”, “Gene” and “Hailfinder” dataset, respectively.

Figure 5.

Box diagram of feature dimensions obtained by four operators

As shown in Fig. 5, the scaled crossover operator in the optimization method of this paper shows a clear advantage in selecting the number of features compared with the other three operators. Specifically, the scaled crossover operator can obtain a smaller subset of features in the dataset. And the experimental results show that different selection probabilities of different operators lead to different performance in different datasets. The higher the selection probability, the smaller the size of the selected features, and the higher the error rate, and vice versa. This shows that the selection probability of the scaled crossover operator has a key impact on the performance of the optimization method in this paper on different datasets. A lower selection probability of the scaling crossover operator for selecting features and obtaining a smaller subset of features will also be able to improve the correct MB selection rate of this paper’s optimization method in the target feature T.