{'fraction_non_alphanumeric': 0.049127220062504764, 'fraction_numerical': 0.010328531138044059, 'mean_word_length': 4.224810832337714, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variable Selection for Nonparametric Learning with Power Series Kernels

In this paper, we propose a variable selection method for general nonparametric kernel-based estimation. The proposed method consists of two-stage estimation: (1) construct a consistent estimator of the target function, (2) approximate the estimator using a few variables by l1-type penalized estimation. We see that the proposed method can be applied to various kernel nonparametric estimation such as kernel ridge regression, kernel-based density and density-ratio estimation. We prove that the proposed method has the property of the variable selection consistency when the power series kernel is used. This result is regarded as an extension of the variable selection consistency for the non-negative garrote to the kernel-based estimators. Several experiments including simulation studies and real data applications show the effectiveness of the proposed method.


Introduction
Variable selection is quite important in various machine learning tasks, to improve the performance, select more cost-effective subset of the features and guarantee the interpretability [9]. Many variable selection methods has been developed for parametric learning including linear models. On the other hand, it is desirable to be able to do variable selection in nonparametric learning although such study has been done only in very limited problems. As a motivating example, in the clinical prognostic problem at multiple hospitals [11], learning of prediction models that will be used at target hospital is performed using integrated data from multiple hospitals. However, since the distribution of covariates among hospitals can differ, adaptation by importance weighting using density-ratio [23] is required to appropriately learn a model. Although the density-ratio estimation can be done accurately by the kernel method [10], all covariates must be observed at all hospitals. Since this is very expensive to use, it is natural that we want to select the variables which contribute to the density-ratio in advance. In this paper, we propose a general variable selection method for kernel-based nonparametric learning.
In this paper, we consider kernel methods not only for regression problems but also density ratio estimation, density estimation, density-ridge estimation, etc. In the kernel methods, we employ the power series kernel [34]. We show that the power series kernel has a desirable property, i.e., the invariant property under the variable scaling on the corresponding reproducing kernel Hilbert space. We prove the variable selection consistency of the kernel methods using the power series kernel and NNG under milder assumptions than [6] in which a kernel variant of the irrepresentable condition is assumed.
Our result is regarded as an extension of the variable selection consistency for the original NNG [33,Corollary 2] or adaptive lasso to the kernel-based estimators.
The rest of the paper is organized as follows. In Section 2, we explain the problem setup and give some definitions. Several kernel nonparametric learning problems are formulated in an unified way.
Section 3 provides the explanation of the proposed method. Section 4 is devoted to the main results.
The variable selection consistency is shown in this section. In Section 5, we show the experimental results for both synthetic and real datasets. Finally in Section 6, we conclude the paper with the discussion on future works.

Related Work
The conventional approach to the variable selection was the information criterion such as AIC or the sequential test called forward and backward step-wise selection. For high-dimensional models, however, the information criterion does not work because of the computational issue, i.e., the combinatorial complexity appears in the choice of variables. Also, the statistical test needs to repeat the computation of sample statistics many times. For high-dimensional linear models, [24] proposed the so-called Lasso estimator, in which the 1 -norm of the coefficients was incorporated into the squared loss.
Also, the non-negative garrote (NNG) was proposed by [3] as a modification of the standard least square (LS) estimator for linear models. In the NNG, each coefficient of the LS estimator is shrunk towards zero, and its intensity is controlled by the 1 -regularization. The variable selection consistency of Lasso was proved under the so-called irrepresentable condition [31,28]. On the other hand, [33] and [30] proved that the NNG has the variable selection consistency without the irrepresentable condition.
For some kernel-based estimators, [1] and [8] employed the NNG as the adaptive scaling for the variable selection. Also, [6] proved the variable selection consistency of the scaled kernel-ridge estimator under a variant of the irrepresentable condition.
Rosasco et al. [15] considered the variable selection problem in the kernel ridge regressions. Instead of the scaling parameters, the authors proposed the regularization based on the derivatives. They proved only the selected variables include the target variables with high probability. However, it is not clear whether the extra variables can be removed with high probability. As the author of [15] pointed out, the variable selection consistency is not completely proved, and that was postponed to the future work.
Several literature deal with similar but different problem setting [7,16,18,29]. In these papers, the sparsity is incorporated into the coefficients in the linear sum of the kernel functions for the kernel learning. Hence, the variable selection concerning the covariates was out of the scope of these papers.

Problem Setup
We briefly introduce kernel methods. See [2,21] for details. Let Z be the domain of the d-dimensional data and k : Z × Z → R be a kernel function. The reproducing kernel Hilbert space (RKHS) associated with k is denoted as H or H k . The RKHS is a linear space consisting of real-valued functions on Z which is used as the statistical model. The inner product and the norm on H are represented as f, g In many learning algorithms, the target function is estimated by minimizing an empirical loss function L(f ) with a regularization term λR(f ) as follows: where L(f ) depends on training data and λ is a positive parameter controlling the capacity of the statistical model. Throughout the paper, we assume that R(f ) = 1 2 f 2 and the empirical loss is expressed as the quadratic form, where a linear operator C : H → H and an element g ∈ H depend on training data. Several examples are shown below. Suppose that L(f ) converges to L(f ) in probability for each f ∈ H as the sample size goes to infinity. We assume that L(f ) has the form of where C : H → H and g ∈ H will depend on the probability distribution of the training data. Suppose that the target function f * ∈ H is the minimum solution of (3).
In this paper we focus on the variable selection problem for kernel methods. The kernel-based estimator f (x) usually depends on all variables in x = (x 1 , . . . , x d ). However, the target function f * may depend on only a few variables. The goal of the variable selection is to detect these variables.

Kernel-Ridge estimator
In regression problems, i.i.d. samples (x 1 , y 1 ), . . . , (x n , y n ) ∼ p(x, y) = p(y|x)p(x) are assumed to be generated from the model y i = f * (x i ) + ε i , where f * is the true regression function and ε i is i.i.d.
random observation noise. The variable selection is important in order to model the relationship between x and y.
The RKHS H is commonly used as the statistical model. The expected squared loss and empirical squared loss for f ∈ H are respectively defined as 1 2 (y − f (x)) 2 p(x, y)dxdy and 1 2n n i=1 (y i − f (x i )) 2 . Let us define C = k(x, ·) ⊗ k(x, ·)p(x)dx and g = yk(x, ·)p(x, y)dxdy ∈ H. Here, ⊗ denotes the tensor product defined as g ⊗ h, f = h, f g for f, g, h ∈ H. Then, the loss function L(f ) = 1 2 f, Cf − g, f equals the expected squared loss up to constant terms. Likewise, L(f ) is expressed The kernelridge estimator f is given as the minimizer of L(f ) + λ 2 f 2 subject to f ∈ H. It is well-known that f − f * = o P (1) holds when the output y is bounded and the regularization parameter λ is appropriately chosen; the detail is found in the proof of Theorem 4 in [4].

Kernel-based density-ratio estimator
Density ratio is defined as the ratio of two probability densities [23]. The density ratio is an important versatile tool in statistics and machine learning, since it appears in many learning problems including regression problems under the covariate-shift, two-sample test, outlier detection, etc. Suppose that we have training samples, x 1 , . . . , x n1 ∼ i.i.d. p and y 1 , . . . , y n2 ∼ i.i.d. q, where p and q are probability densities on the domain Z. Our goal is to estimate f * (z) = q(z)/p(z).
Let us consider the variable selection of the density ratio. For the vector z = (z a , z b ), p(z) and q(z) are decomposed into conditional probabilities and marginal ones such as p(z a |z b )p(z b ) and q(z a |z b )q(z b ).
The variable selection of the density ratio is closely related to the identification of the conditional probability that p and q have in common.
A kernel-based density-ratio estimator called KuLSIF was proposed in [10]. The empirical loss of f is n2 n2 j=1 f (y j ). As n = min{n 1 , n 2 } tends to infinity, the empirical loss converges to L(f ) = Z { 1 2 f (z) 2 p(z)−f (z)q(z)}dz due to the law of large numbers. The minimizer of L(f ) is nothing but f * = q/p. The empirical quadratic loss is expressed by C = 1 n1 n1 i=1 k(x i , ·)⊗k(x i , ·) and g = 1 n2 n2 j=1 k(y j , ·), and the expected loss L(f ) by C = k(x, ·) ⊗ k(x, ·)p(x)dx and g = k(y, ·)q(y)dy up to constant terms. [10] proved the statistical consistency of the KuLSIF estimator in the L 2 norm. In Appendix A, we prove the statistical consistency in the RKHS norm.

Kernel-based density estimator and density-ridge estimator
The kernel-based density estimator using infinite dimensional exponential models has been studied in [20]. The problem is to estimate the function f * of the probability density p(z) ∝ exp(f * (z)) using the i.i.d. samples z 1 , . . . , z n ∈ Z from p(z). The variable selection of the probability density on the bounded domain Z is nothing but the identification of variables whose marginal probability is the uniform distribution on the domain.
The estimator f of f * is obtained by minimizing the Hyvärinen score with the regularization. Let ∂ a f be the derivative of the function f (z) w.r.t. z a , i.e., ∂f ∂za (z). Likewise, ∂ aa f denotes the second-order derivative of k w.r.t. the a-th argument. The empirical loss L(f ) for the Hyvärinen score is defined . The reproducing property for the derivative, ∂ a f (z) = f, ∂ a k(z, ·) , is useful to conduct the calculation of the operator C. Note that ∂ a k(z, ·) ∈ H holds for the kernel function k [32]. Likewise, the expected quadratic loss L(f ) is expressed by C = d a=1 ∂ a k(z, ·) ⊗ ∂ a k(z, ·)p(z)dz and g = − d a=1 ∂ aa k(z, ·) p(z)dz ∈ H. When a proper boundary condition is assumed, the expected quadratic loss derived from C and g equals 1 2 d a=1 |∂ a f (z) − ∂ a f * (z)| 2 p(z)dz up to constant terms. Thus, the Hyvärinen score is regarded as the mean square error for derivatives. The estimator using the Hyvärinen score and an appropriate regularization parameter λ has the statistical consistency f − f * = o P (1) under a mild condition [20].
The density-ridge estimator is related to the above density estimator. The target is to estimate where ∂ I is the differential operator ∂ k ∂zi 1 ,...,∂zi k with the set (or multiset) of non-negative integers I = {i 1 , . . . , i k }. The above density estimator corresponds to the case that the set I is a singleton. The ingredients of the quadratic loss, C, g, C and g, are defined from the derivative of the kernel function, ∂ I k(z, ·). The estimated function f is used to extract the "ridge structure" of the probability density p(z) that usually has a complex low-dimensional structure.
The variable selection for the density-ridge estimator is important to boost the estimation accuracy of the ridge and to reduce the computational cost.

Kernel Methods with Adaptive Scaling
For variable selection, we incorporate adaptive scaling parameters to the variables in the RKHS model.
As the adaptive scaling parameter, we employ Breiman's Non-Negative Garrote (NNG) [3]. The original NNG is used to estimate the sparse vector β of the linear regression model y = x T β + , where is the observation error. The least mean square estimator β 0 is mapped to ξ • β 0 , where the non-negative parameter ξ is called the garrote parameter and the operator • denotes the element-wise product of two vectors, i.e., the Hadamard product. The optimal garrote parameter ξ is found by minimizing the empirical squared loss with the non-negative constraint ξ ≥ 0 and the 1 -regularization ξ 1 ≤ c.
Eventually, the sparse vector ξ • β 0 is obtained as the estimator of the coefficient vector.
We incorporate the NNG into kernel methods. In order to induce the sparse estimator, the NNG seems to be more adequate than lasso-type estimator, since lasso estimator is basically available to the feature selection of linear regression models [25]. The linear model with the garrote parameter is expressed as (ξ • x) T β 0 . Likewise, given the kernel-based estimator f (z) for f ∈ H, the garrote parameter ξ is introduced as the form of f ξ (z) := f (ξ • z). Both f ∈ H and ξ can be found by minimizing the empirical loss in which f is replaced with f ξ . Here, we propose a simplified two-stage kernel-based estimator with NNG. The detail is presented in Algorithm 1. In the algorithm, η is a positive regularization parameter that controls the sparsity of the variable selection. Using the representer theorem, the estimator is typically expressed as f ξ (z) = i α i k(z i , ξ • z), where z i is a data point and α i is the estimated parameter in Step 1.
Step 1: Find the kernel-based estimator f by solving (1).
we introduce the upper constraint ξ ≤ 1 that does not appear in the original NNG. This is because we need the contraction condition in Section 3.2 to ensure that the domain of f ξ is properly defined from that of f . The estimated function f ξ depends only on the variables having positive garrote parameter.
In step 2 of Algorithm 1, one can use the standard optimization methods such as the limited-memory BFGS method with the box-constraint. Usually, the objective function is not convex. Practical methods including the multi-start technique should be implemented.

Power Series Kernels and its Invariant Property
The statistical model of the learning algorithm is expressed as In general, the class of power series kernels (PS-kernels) [34] is available in our method. The The polynomial and exponential kernels are included in the class of the PS-kernels. The native in the summation runs on multi-indices with w α > 0. As shown in [34], H k is the RKHS with the inner product The invariant property holds for the RKHS endowed with the power series kernel. Indeed, for Moreover we have f ξ ≤ f . Some other properties of the power series kernels are presented in Appendix B.

Theoretical Results : Variable Selection Consistency
Let us consider the variable selection consistency of Algorithm 1. All proofs are found in Appendix C.
Suppose that the target function f * (z) for z = (z 1 , . . . , z d ) essentially depends on the variables The variable selection consistency means that for the estimated garrote parameter ξ = ( ξ 1 , . . . , ξ d ), the probability of the event {j : ξ j > 0} = {1, . . . , s} tends to one as the sample size n goes to infinity.
Let us define some notations and introduce some assumptions. In the following we assume that the kernel function and its derivatives are bounded by the constant κ > 0, i.e., sup z∈Z k(z, z) ≤ κ and sup z∈Z ∂ a ∂ d+a k(z, z) ≤ κ hold for a = 1, . . . , d. We use the similar inequality up to required order of derivatives. Then, the convergence of the RKHS norm Indeed, the generalized reproducing property [32] leads Under the above assumptions, we prove the variable selection consistency of ξ. The proof is in Appendix C.1 follows the standard argument of the statistical consistency of Mestimators shown in [27,Theorem 5.7]. The order of η n should be greater than δ n and δ n in order to draw ξ 0 to 0. All kernel-based estimators in Section 2 satisfies Assumptions 1. In Appendix C.2, we prove that for each estimator in Section 2, the condition (a) in Assumptions 2 holds under a mild assumption.
We show sufficient conditions of Assumptions 2 (b) and 3, when the quadratic loss functions in (2) and (3)  Loss functions for the estimators in Section 2 are expressed using the above C, C, g and g. The stochastic convergence property of C − C and g − g is guaranteed from Theorem 7 of [14]. Hence, the sufficient condition in Proposition 1 is satisfied.
Let us show another sufficient condition of Assumption 3. We deal with more general operators C and C, while we need an additional smoothness assumption on the target function f * . The operators C, C and the elements g and g in Proposition 2 satisfy the inequalities of the assumption in Lemma 3.

Experimental Results
In this section we show the empirical performance of the kernel-based density-ratio estimator. In section 5.1, we conduct the synthetic data analysis to mainly evaluate the variable selection consistency.
In section 5.2, we analyze three real datasets, PCAS dataset [11], diabetes dataset [22] and Wisconsin breast cancer dataset. The latter two datasets are published at the UCI machine learning repository.
We see that our method can select a practically interpretable subset of features.

Synthetic Data Analysis
We report numerical experiments using synthetic data. The purpose of this simulation is to confirm the statistical accuracy of the KuLSIF with NNG in terms of the variable selection in density ratio.
Suppose that x 1 , . . . , x n were generated from the d-dimensional normal distribution with mean µ and the variance-covariance matrix I d , i.e., N d (µ, I d ). This distribution corresponds to p in the denominator of the density ratio. Likewise, suppose that y 1 , . . . , y m were generated from the probability Our methods were compared with the sparse logistic regression (SLR) estimator, that is given as the logistic regression with the 1 -regularization. Originally, the SLR estimator is the learning algorithm for sparse classification problems. Suppose that the samples from p(x) (resp. q(x)) have the label +1 (resp. −1). Then, the ratio Pr(−1|x)/ Pr(+1|x) of the estimated logistic model Pr(+1|x) = 1/(1 + exp(−(β 0 + β T x))) is approximately proportional to the density ratio q/p if the statistical model is specified. The 1 -regularization β 1 to the weight vector induces the sparse solution as well as the lasso estimator. The regularization parameter for the 1 regularization in SLR was determined by the cross validation. In our method, the hyper-parameter λ of KuLSIF in step 1 was set to 1/(n ∧ m) 0.9 , which guarantees the statistical consistency of KuLSIF under a proper assumption [10]. In step 2 of KuLSIF with NNG, the regularization parameters η varies from 0.01 to 1. For the computation, we used the R language. The glmnet package of the R language was used for the SRL. For the optimization in step 2 of KuLSIF with NNG, we used the limited-memory BFGS method with the box-constraint implemented in optim function in R.
We report the result of the numerical experiments. The test loss of the estimated density ratio, w, was measured by the shifted squared loss, L( w) = 1 2 p(x) w(x) 2 dx − q(y) w(y)dy. The variable selection accuracy was evaluated by the false positive rate (FPR) and false negative rate (FNR). In each problem setup, the simulation was repeated 30 times. The Table 5.1 shows the test loss, the FPR, and FNR averaged over the repetitions.
The result indicates that the KuLSIF with NNG using the exponential kernel performs better than the other method in terms of the variable selection if the regularization parameter η is properly chosen.
When we use the Gaussian kernel with our method, we cannot expect the ability of the variable selection.
The variables selected by the SLR had relatively high FPR. In the experiment, the regularization parameter of SLR was chosen by the cross validation under the empirical 0/1 loss. In this case, a bigger subset of variables tends to be chosen. The test loss of the SLR is larger than the other methods, since the classification-oriented learning algorithms are considered not to be necessarily suit to the density ratio estimation.

Real Application : Covariate-shift Adaptation for Binary Classification
In this section, we show the results of applying the proposed method to learning under the covariate shift problem [19,23]. The covariate shift is a phenomenon in which the distribution of covariates changes between learning phase and test phase, and it has been proved that the covariate shift can "adapt" by weighting the sample with density-ratio. We consider the binary classification problem under the covariate-shift, and compare the test accuracy of the five scenaros:(i) no adaptation, (ii) adapting with the density-ratio using all variables, (iii) adapting with the density-ratio using selected variables by kernel NNG (proposed method), (iv) adapting with the density-ratio using selected variables by Lasso and (v) adapting with the density-ratio using selected variables by sequential forward selection (SFS). Here we employ the logistic regression to learn a classification model. In this study, we analyzed four medical datasets: post-cardiac arrest syndrome (pcas) [11], chronic kidney disease (ckd), cervical cancer (ccancer) and cortex nuclear (cnuclear) [5]. The characteristics of each dataset are summarized in Table 2. In the column of the sample size (#samples) and the number of features (#fearures), those after preprocessing, such as missing value correction, are entered. In each dataset, we stratified the data into two groups by a certain feature, and learning was performed on one group, and the classification accuracy was evaluated on the other. The stratification factor is summarized in the last column of  Table 2. The learning phase consists of kernel density-ratio estimation (KuLSIF) step and weighted empirical risk minimization (weighted logistic regression) step. All experiments are implemented in Python 3.6.1. We use scikit-learn [12] for Lasso and mlxtend [13] for SFS. Similar to the previous section, we used the limited-memory BFGS method with the box-constraint in minimize function in Scipy. The results of the experiments are summarized in Figure 1, 2 and Table 3. Figure 1 shows the differences of the density ratio estimated by selected variables (red, proposed method) and by full variables (blue) for the pcas dataset. The differences of each histogram are measured by root mean The bold symbol represents the best performance. As shown in the results, the proposed method appropriately selected a small number of variables without deteriorating the classification accuracy. Figure 2 shows the comparison of the estimated density-ratio between selected variables (by kernel NNG, Lasso and SFS) and full variables. The differences of each histogram are also measured by RMSE.
Although the Lasso and the SFS achieved the small RMSE in ckd and ccancer datasets, both methods could not select the features appropriately. In pcas dataset, SFS shows the best performance in both

Concluding Remarks
This paper provided a unified variable selection method for nonparametric learning by the power series kernel. The proposed method can be applied to kernel methods using the quadratic loss such as kernel-ridge regression and kernel density-ratio estimation. Theoretically, we proved the variable selection consistency under mild assumption thanks to the property of the power series kernel.
Experimental results showed that our method efficiently worked for the variable selection for kernel methods on both synthetic and real-world data. In the second step of our method, we need to solve nonconvex optimization problem for variable selection. In numerical experiments, non-linear optimization algorithms such as the limited-memory BFGS method showed a good performance. Future work includes the development of more efficient optimization methods to deal with learning problems with massive data sets.
A Statistical consistency of kernel-based density-ratio estimator The proof follows the convergence analysis developed by [4,20]. The estimator f is the minimum solution of the function J(f ) over H, where Thus, we have f = ( C + λ n I) −1 g. When the sample size n = min{n 1 , n 2 } tends to infinity, J(f ) converges to J(f ) that is defined as The minimizer of J(f ) is expressed as f λn = (C + λ n I) −1 g. Note that the true density ratio is given We define , Proposition A.4 in [20] is used to bound S 1 , S 2 and S 3 as follows, Using the above bounds, we obtain The asymptotic order of A(λ n ) was revealed by Proposition A.3(i) in [20]. Indeed, (I) the proposition shows that A(λ n ) → 0 as λ n → 0 if f 0 ∈ Range(C); (II) if f 0 ∈ Range(C β ) for β > 0, it follows that
In terms of the derivative ∂f ξ ∂ξi , we prove some formulae. Proof. Suppose As shown in the proof of (4), the derivative of f ξ is given by the sum of the derivative of each term c α ξ α x α . Let B be a positive number such that 0 ≤ α 1 ξ α1−1 is confirmed by In general, we have Lemma 5. Let us consider the RKHS H k endowed with the PS kernel k. We define e i ∈ R d as the unit vector with one in the i-th position and zeros otherwise. For f ∈ H k and ξ ∈ (0, 1) d , we have for i = 1, . . . , d.

C.2.2 Kernel-based density-ratio estimator
The condition (a) in Assumption 2 is confirmed in the same way as the kernel-ridge estimator. Hence, we omit the details.

C.2.3 Kernel-based density estimator
For f * ∈ H k , the probability density p(z) ∝ exp(f * (z)) is strictly positive on the compact domain.