Latent regularization for feature selection using kernel methods in tumor classification

By a Multiple Kernel Learning process built to optimize the KTA score a subset of features is selected and a kernel is obtained. This subsection defines how to select features with MKL by a supervised approach.
Given $\mathcal{X}\subseteq\mathbb{R}^{n}$ a n-dimensional space , a function $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ which is symmetric and semi-positive definite meets the condition

where $\phi$ is a mapping function from $\mathcal{X}$ to a high dimensional feature Hilbert space $\mathcal{H}$ such that

and $\mathcal{H}$ is defined as the Reproducing Kernel Hilbert Space (RKHS) of $k$ [29]. The function $k$ computes the inner product of a pair of vectors in $\mathcal{H}$ the Hilbert Space. Considering a set of $m$ labeled samples such that

where $x_{i}\in\mathcal{X}$ and $y_{i}\in\left\{-1,+1\right\}$ the Kernel Matrix or Gram Matrix $G$ is defined as a $M\times M$ matrix with entries $K_{ij}$. Every entry of the Kernel or Gram Matrix is defined as

The kernel function can be thought as an information bottleneck that concentrates the information required to perform a learning task [30]. Kernels can be thought as similarity functions where an output close to $0$ means orthogonal samples or high dissimilar and an output close to $1$ means similar samples in the Hilbert space. In this work we will use the Radial Basis Function (RBF) Kernel which is defined as

where $\gamma$ is a scale parameter.
Given two valid kernels $K_{1}$ and $K_{2}$ over a set of samples $M$, the alignment $A$ between both kernels is defined as

and measures the similarity between the two kernels using the same sample set $M$ [20][31]. Here the operation $\left\langle.\right\rangle_{F}$ is the Frobenious Inner product. The alignment can be thought as how similar both kernels maps the samples from $M$. If tumor labels are used, then $K_{2}$ represents an ideal Kernel matrix or target $K_{yy}$ where $K_{yy}=+1$ if $y_{i}=y_{j}$ and $K_{yy}=0$ if $y_{i}\neq y_{j}$. The alignment of a kernel $K$ built on $x$ and the target $K_{yy}$ can be expressed as

known as Kernel Target Alignment (KTA) score. The higher the KTA between a given kernel matrix $K$ and its target $K_{yy}$ the higher the inter-cluster separation between the two classes.
Multiple Kernel Learning (MKL) allows the combination of simple kernels to build a more complex one with higher KTA. MKL is defined as the linear combination of multiple kernels to build a final one [7] and can be expressed as

where the vector parameter $\mu$ corresponds to the weight $\mu_{i}>0$ of each kernel $k_{i}$ and it is directly related to the importance of each kernel in the final solution. Despite there are different objective functions to calculate the weights of the MKL model, like the classification accuracy of a support vector machine, in this work the MKL model is built with the objective to maximize the KTA. This means that the resulting kernel from the combination of the initial ones will present a higher inter-class separability than each individual kernel on its own. The resulting KTA for the kernel $\boldsymbol{K_{\mu}}$ is $\mathit{A}\left(\boldsymbol{K_{\mu}},K_{yy}\right)$ and computed as detailed in equation 6.
This work uses a greedy strategy [16] to compute the weights $\mu$ of the resulting $\boldsymbol{K_{\mu}}$ kernel by combining just two kernels at each iteration while maximizing the KTA of the resulting kernel. This strategy allows a simpler computation of the $\mu$ vector by obtaining the weights values where the derivative of the alignment becomes $=0$ as optimal condition for each partial derivative. Solving the MKL process with two kernels at each iteration $t$ means that the weight vector is $\boldsymbol{\mu}=\left[\mu_{\alpha},\mu_{\beta}\right]$. Since the solution with two kernels is convex [32] the values of $\mu_{\alpha}$ and $\mu_{\beta}$ can be obtained from

subject to $\mu_{\alpha},\mu_{\beta}\geq 0$. From a set of $N$ kernels, the greedy approach starts at iteration $t=0$ by choosing the kernel $K_{i}$ from $N$ with highest $KTA_{i}$. Then assign $K_{i}$ to a empty list $P$ containing the kernels which are part of the solution and makes it the current solution $K_{\mu}^{(t)}=K_{i}$. Then for the following iterations $t=1...p$ the method combines iterativelly the current $K_{\mu}^{(t)}$ and a new kernel $K_{j}$ from $N$ with $i\neq j$ to obtain a new solution $K_{\mu}^{(t+1)}=\mu_{\alpha}K_{\mu}^{(t)}+\mu_{\beta}K_{j}$ that maximizes the overall KTA as $A(K_{yy},K_{\mu}^{(t+1)})$. This process prioritizes in adding first to the solution $P$ the kernels that increase the most the overall KTA, thus the kernel selected first will have the highest weight $\mu_{i}$. The following kernels to be selected will have a lower weight $\mu_{i}$ since their contribution to maximize the overall KTA as $\Delta KTA$ is going to be lower than the previous added kernel. The KTA increment at each iteration will be $\Delta KTA^{(t)}>\Delta KTA^{(t+1)}...>\Delta KTA^{(p)}$. This means that when a number $p$ of kernels is requested to be in the solution weights tend to decrease with $p$ as $\mu_{0}>\mu_{1}>....>\mu_{p}$.
Given a dataset of $m$ samples characterized by $n$ features, this work proposes to build $n$ feature-wise kernels $K_{i}$. Feature-wise kernel $\gamma_{i}$ parameters are selected based on the criterion of maximizing the alignment of each kernel. Then using the proposed MKL method a subset of feature-wise kernels $P$ is iterativelly selected and combined to increase the overall KTA of the resulting kernel $\boldsymbol{K}_{\mu}$. Only the feature-wise kernels that increases the KTA are included in the final kernel. This approach leads to a sparse solution where the number of selected features $p$ associated to the selected feature-wise kernels is $p<<n$. The desired output of the greedy MKL approach is a reduced set of $p$ features and a kernel $\boldsymbol{K}_{\mu}$ that improve the inter-cluster distance between samples of different tumor classes and thus improve the support vector classification. The positive values of the weighting vector $\mu_{i}>0$ indicates the feature importance on the result and the selected features.
Once $\boldsymbol{K}_{\mu}$ is built it is used as a custom kernel function in a support vector classifier [29] for binary classification in tumor stage and patient survival problems.

In the previous section a feature selection method based on multiple kernel learning and kernel target alignment is described. This method selects features using a supervised kernel matrix $K_{yy}$ where the separation between classes is ideal as target. But this objective is very specific and since the number of samples is small it can lead to significant different solutions when changing samples in the training set. Thus it may be interesting to be less specific and more conservative about the structure of the data. For this reason a relaxation of this supervised constraint is proposed in this work by mixing the target kernel $K_{yy}$ with a $K_{z}$ kernel built on the latent space $Z$ learned by a dimensionality reduction algorithm based on a non-linear transformation $\phi_{z}(x)$. Note that the $K_{yy}$ is built from the tumor labels and the $K_{z}$ from the extracted latent variables, thus the first one is built by a supervised approach and the latter by an unsupervised one. The mixture of both kernels forms a new target kernel $K_{\delta}$ that has supervised and unsupervised information from the training samples. For the unsupervised and nonlinear dimensionality reduction transformation $\phi_{z}(x)$ any algorithm can be used such as Kernel-PCA, Autoencoders or t-SNE. The key aspect is to learn using an unsupervised criteria, a latent space from the training samples to capture the structure of the data and to mix it with the supervised labels as a mean to relax the supervised constraint.
Most supervised feature selection models $f$ are expressed as

where a subset $p$ of selected features minimizes a loss function between the model output $f(x)$ and the true labels $y$. This work proposes a feature selection model that learns not only from the true labels $y$ but also from the latent structure of the training data as

where $z=\phi(x)_{z}$ is the latent space obtained from a nonlinear and unsupervised mapping function. The mixture between the labels $y$ and $z$ creates a hybrid target composed by a supervised and unsupervised approach. In this work to extract a latent space $z$ the function $\phi(x)_{z}$ is learned by a Kernel-PCA method [8].
This method is based on a nonlinear mapping $\phi_{z}(x)$ from $\mathcal{X}$ to $\mathcal{H}$ where the standard Principal Component Analysis (PCA) algorithm is performed.
The PCA is a linear dimensionality reduction method and is defined as the orthogonal projection of the training samples into a low dimensional space such that the variance of the projected samples is maximized [33]. From a set of samples ${x_{m}}$ where $m=1,...,M$ characterized in $n$ dimensions, the goal of PCA is to map $x_{m}$ into a low dimensional space of dimension $p$ such that $p<n$. The PCA projects the data into the p low dimensional space that maximizes the variance. The PCA transformation is obtained by applying the spectral decomposition to the covariance matrix $C$ of the training data $x_{m}$ such that $U\Lambda=CU$. The $U$ matrix contains stacked vectors $u_{1},...,u_{p}$ where $u_{i}$ is the ith eigenvector corresponding to the ith largest eigenvalue $\lambda_{i}$ in the $\Lambda$ matrix. Then by selecting the first $p$ eigen-vectors a new $n\times p$ matrix ${U}^{\prime}$ is obtained and works as the desired subspace. Then by computing $X{U}^{\prime}=Z$ the original data is projected by ${U}^{\prime}$ into $Z$ of dimension $p$.
To make PCA nonlinear the Kernel Methods introduced in the previous section are used to define the Kernel-PCA in order to perform an implicit PCA in the projected Hilbert space $\phi(x)$ via the kernel trick [34]. In this case the spectral decomposition is applied to a Kernel Matrix $K$ such that

where $U$ and $\Lambda$ are the matrices containing the eigenvectors and eigenvalues respectively and $m$ the number of samples. Then the resulting coordinates $(l_{1},...,l_{p})$ known as kernel principal components are calculated by

and this is equivalent to perform a nonlinear PCA in the original input space.
Once the kernel-PCA is computed with the training samples $\boldsymbol{x}_{tr}$ these are mapped to $\boldsymbol{z}$ and a kernel $k_{z}$ is built on the latent space. $k_{z}$ captures the unsupervised structure of the training samples in the latent space. Then by a linear combination of the supervised $k_{yy}$ and unsupervised $k_{z}$ kernels a new hybrid target kernel $K_{\delta}$ can be created

This new kernel contains a linear combination of both the supervised labels and unsupervised latent variables ruled by a new parameter named as Mixture Coefficient $\delta\in[0,1]$.

Figure 1 shows the KLR-FS pipeline. From equation 10 three scenarios are possible. When $\delta=1$ then $K_{\delta}=K_{yy}$ and corresponds to the supervised kernel described in the previous section and the feature selection process is completely supervised. On the other hand when $\delta=0$ then $K_{\delta}=K_{z}$ and it represents the unsupervised structure provided by the kernel-PCA of the training samples. When $\delta=0$ the target $K_{\delta}$ does not contain any supervised information and the MKL process is unsupervised. Finally, every value of $0<\delta<1$ corresponds to a kernel $K_{\delta}$ that has a mixture of supervised and unsupervised components of the problem. Then by this approach a latent regularization of the supervised feature selection problem is possible. The hypothesis stated in this work is based on the assumption that the features selected by MKL in a mixture of supervised labels and unsupervised latent variables results in a higher generalization ability and thus in higher classification performance.
Finally, once the MKL step is done a $\boldsymbol{K}_{\mu}$ kernel and a subset $P$ of selected features are obtained as result. Then the $\boldsymbol{K}_{\mu}$ kernel is used as the kernel in a support vector classification [35] to classify tumor profiles.

This work is focused on binary classification of tumor samples by using reduced feature subsets selected by KLR-FS. The classification task gives an idea about the predictive power of the selected features. Each binary classification problem has two classes: a Positive and a Negative class. To evaluate a classification model on new independent test samples four outcomes are possible: True Positive (TP), False Negative (FN), True Negative (TN) and False Positive (FP). Given the Positive class the TP counts the number of samples well classified while the FN counts the wrong classified by the model. Given the Negative class the TN counts the number of samples well classified while the FP counts the wrong classified ones.
Another way to measure the performance of a classifier is to estimate the Area Under the Receiving Operation Characteristic Curve (ROC) [36] or Area Under the Curve (AUC). The ROC represents the TP rate as a function of the FP rate. The AUC is a score used to measure the overall performance of a binary classifier. It ranges from $0.5$ to $1$. An AUC $=0.5$ represents the performance of a random classifier while an AUC $=1$ corresponds to a perfect classifier. Unlike Accuracy or other classification metrics based on a specific decision threshold defined a priori by the classifier the AUC is independent of the classification threshold and is a measure of a continuous output. This makes the AUC a robust score. The AUC is used to evaluate the classification performance of a binary classifier trained on the selected features by KLR-FS. For all the experiments a Support Vector Classifier is trained using the selected features of each method.
To evaluate the robustness and quality of the feature selection process the Redundancy Rate (RED) [37] [14] metric is used. The RED score measures the mean value of absolute correlation between features is computed as

where $\rho_{ij}$ is the correlation score between the ith and the jth selected variables. The RED metric takes values between $0$ and $1$. A low value of RED indicates that the majority of the selected features within the subset $P$ have low linear correlation between each other and thus a low redundancy is expected within the selected feature set which can be interpreted as a high quality of selection. On the other side values of RED close to $1$ correponds to a subset of features with high redundancy which is a non desired output.

The proposed feature selection method for classification is compared with four other supervised feature selection methods. These methods are Minimum Redundancy Maximum Relevance (mRMR)[38], the High-Dimensional Feature Selection by Feature-Wise Kernelized Lasso (HSIC-Lasso) [14], the SVM Recursive Feature Elimination (SVM-RFE)[13] and the univariate ANOVA filter [39].
The mRMR method attempts to discard redundant features while keeping the features with highest relevance to the target labels $y$. The objective is to select a feature subset that best characterizes the statistical property of a target label [38]. The selection has the constraint that selected features are mutually as dissimilar to each other as possible, but marginally as similar to the target label.
The HSIC-Lasso method [14] is a feature-wise kernelized Lasso that captures non-linear dependency between input features and target labels.
The SVM-RFE is a greedy feature selection method [40][41] that generates a ranking list of features and selects a subset of the top-ranked features. The ranking is built by a feature weight vector $w$ obtained from the parameters of the hyperplane decision function of a SVM classifier and the top $p$ features are selected.
Finally, ANOVA t-test is an univariate feature selection method that gives a score to each feature based on the p-value obtained from a statistical test between the feature $i$ and the class label $y$.
During support vector classification for benchmark methods a RBF Kernel has been used with a grid search on gamma hyperparameter $\gamma=[0.01,0.1,1,10]$. For the KLR-FS the resulting $k_{\mu}$ kernel is used. A grid search for the C hyperparameter $C=[0.1,1,10,100]$ has been used for all kernels. Hyperparameters have been selected by 5-fold cross validation on train set to select the best model.

The initial sample set has been randomly split between $80\%$ as train and $20\%$ as test. By using just the training samples all the gene expression features have been auto-scaled with 0 mean and unit variance. Test samples have been auto-scaled using the transformation learned from train samples. Then by using training samples the feature selection methods are applied and a subset of features are selected. With the selected features a support vector classifier is trained and tuned by 5-fold cross validation within train set for hyperparameter selection and evaluated on the test set.

This subsection details how the feature selection process of KLR-FS behaves with different settings of the mixture coefficient $\delta$. As explained in the materials and methods section by varying the values of $\delta$ different target kernels $K_{\delta}$ are obtained, each one with a different mix between the unsupervised $K_{z}$ and the supervised $K_{yy}$ kernels. A set of six values of the Mixture Coefficient $\delta=[0,0.2,0.4,0.6,0.8,1.0]$ are evaluated by increasing the dependency to the labels. Each value of $\delta$ determines a feature selection model. Then the resulting kernel $K_{\delta}$ is used for support vector classification and evaluated on the test set.
Figure 2 shows the AUC results on classification by using features selected for different values of $\delta$ on the three datasets. The AUC score peaks between $\delta=0.4$ and $\delta=0.6$ in all the cases. This results evidence the importance of the latent regularization in classification and means that the maximum AUC score using the KLR-FS features is obtained with a mixture between supervised labels and unsupervised latent structure.
Figure 3 shows how the RED score varies for different values of the Mixture Coefficient $\delta$ on the corresponding selected features and how the RED score is reduced by the latent regularization. For the BRCA-US and PACA-CA datasets the lowest values of RED are obtained for $\delta=0.4$ and $\delta=0.6$ and match with the same range of $\delta$ values with maximum AUC. For the SMK-CAN dataset the lowest RED scores are obtained between $\delta=0.6$ and $\delta=1$ and only when the number of selected features $p=20$ and $p=30$ the lowest RED score is at $\delta=1$. Figure 3 shows that when the number of selected features $p$ is small the mixture parameter $\delta$ has a stronger impact.

In this subsection the KLR-FS is compared with the benchmark methods for different number of selected features $p$. For performance estimation the results of the classification and selection tasks are averaged among all the random iterations. The feature selection is applied on train samples and the classifier is trained on the selected features to finaly classify the test samples. The mixture coefficient used for KLR-FS for benchmark is $\delta=0.6$ since it appear to be the best combination according to AUC and RED scores in the previous section.

Figure 4 shows the classification performance on test set for different numbers of $p$ selected features by each method. For the BRCA-US and PACA-CA datasets the KLR-FS shows the highest mean AUC score and the lowest classification variance for every number of selected features. For the SMK-CAN dataset the KLR-FS has the highest classification performance with $p=10$ features and for larger number of features it is outperformed by the RFE method.
In the three datasets the classification results using the KLR-FS features are the most constant as the value of $p$ increases. According to section 3.1 the resulting $K_{\mu}$ assigns a larger weight $\mu_{i}$ to the variables that are selected first since these are the ones that increments more the overall alignment $\Delta A(K_{\delta},K_{\mu})$, while to those variables selected later the corresponding weight is much smaller. The KLR-FS results on Figure 4 suggest that these classification problems can be solved by selecting the first 10 variables since selecting more variables does not introduce a significant improvement in the classification results.

Figure 5 shows the RED score of each subset of selected features for different values of $p$. For the BRCA-US and PACA-CA datasets the lowest score is obtained by the mRMR method with a $RED<0.1$ followed by RFE, KLR-FS and HSIC-Lasso. For the SMK-CAN dataset the lowest RED score is obtained by the RFE method followed by HSIC-Lasso and KLR-FS while the highest corresponds to mRMR and ANOVA. The ANOVA shows the highest RED score in all datasets. KRL-FS and HSIC-Lasso shows a similar behaviour with a $RED<0.2$ in the BRCA-US and PACA-CA datasets and a $RED<0.4$ for the SMK-CAN datasets.
This section compares the KLR-FS method with $\delta=0.6$ and four feature selection methods. The main objective of this work is to select features to improve the tumor classification. It is observed that the classification performance of the KLR-FS is the highest for the majority of the cases.