{'fraction_non_alphanumeric': 0.04369889018691589, 'fraction_numerical': 0.012120327102803738, 'mean_word_length': 4.087852897473997, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Gaussian Mirror for Controlled Feature Selection in Neural Networks

Deep neural networks (DNNs) have become increasingly popular and achieved outstanding performance in predictive tasks. However, the DNN framework itself cannot inform the user which features are more or less relevant for making the prediction, which limits its applicability in many scientific fields. We introduce neural Gaussian mirrors (NGMs), in which mirrored features are created, via a structured perturbation based on a kernel-based conditional dependence measure, to help evaluate feature importance. We design two modifications of the DNN architecture for incorporating mirrored features and providing mirror statistics to measure feature importance. As shown in simulated and real data examples, the proposed method controls the feature selection error rate at a predefined level and maintains a high selection power even with the presence of highly correlated features.


Introduction
Recent advances in deep neural networks (DNNs) significantly improve the state-of-the-art prediction methods in many fields such as image processing, healthcare, genomics, and finance. However, * The first two authors contributed equally to this paper. † The first two authors contributed equally to this paper.
complex structures of DNNs often lead to their lacks of interpretability. In many fields of science, model interpretation is essential for understand the underlying scientific mechanisms and for rational decision-making. At the very least, it is important to know which features are more relevant to the prediction outcome and which are not, i.e., conducting feature selection. Without properly assessing and controlling the reliability of such feature selection efforts, the reproducibility and generalizability of the discoveries are likely unstable. For example, in genomics, current developments lead to good predictions of complex traits using millions of genetic and phenotypic predictors. If one trait is predicted as "abnormal," it is important for the doctor to know which predictors are active in such a prediction, whether they make medical and biological sense, and whether there exist targeted treatment strategies.
In the classical statistical inference, some pvalue based methods such as the Benjamini-Hochberg [5] and Benjamini-Yekutieli [6] procedures have been widely employed for controlling the false discovery rate (FDR). However, in nonlinear and complex models, it is often difficult to obtain valid p-values since the analytical form of the distribution of a relevant test statistics is largely unknown. Approximating p-values by bootstrapping and other resampling methods is a valid strategy, but is also often infeasible due to model non-identifiability, (i.e., the existence of multiple equivalent models, such as a DNN with permuted input or internal nodes), complex data and model structures, and high computational costs. In an effect to code with these difficulties, [2] proposed the knockoff filter, which constructs a "fake" copy of the feature matrix that is independent of the response but retains the covariance structure of the original feature matrix (also called the "design matrix"). Important features can be selected after comparing inferential results for the true features with those for the knockoffs. To extend the idea to higher dimensional and more complex models, [9] further proposed the Model-X knockoff, which considers a random design with known joint distribution (Gaussian) of the input features. Recently, [34] develops the PCS inference to investigate the stability of the data by introducing perturbations to the data, the model, and the algorithm. Dropout [29] and Gaussian Dropout can be viewed as examples of model perturbations, which are efficient for reducing the instability and elevating out-of-sample performances.
For enhancing interpretability, feature selection in neural networks has attracted some recent attention. To this end, gradients with respect to each input are usually taken as an effective measurement of feature importance in convolutional and multi-layer neural networks as shown in [18]. [27] selects features in DNNs based on changes of the cross-validation (CV) classification error rate caused by the removal of an individual feature. However, how to control the error of such a feature selection process in DNNs is unclear. A recent work by [22] utilizes the Model-X Knockoff method to construct a pairwisely-connected input layer to replace the original one, and the latent "competition" within each pair separates important features from unimportant ones by the measurement of knockoff statistics. However, the Model-X method requires that the distribution of X be either known exactly, or known to be in a nice distribution family that possesses simple sufficient statistics [20]. Also, the construction procedure for the knockoffs in [22] is based on the Gaussian assumption on the input features, which is not satisfied in many applications such as genome wide association studies with discrete input data. In addition, the DeepPINK does not fit exactly the purpose of network interpretation since the input layer is not fully connected to hidden layers, and thus some of the connections in the "original" structure may be lost in this specific Pairwise-Input structure.
In this paper, we introduce the neural Gaussian mirror (NGM) strategy for feature selection with controlled false selection error rates. NGM does not require any knowledge about or assumption on the joint distribution of the input features. Each NGM is created by perturbing an input feature explicitly. For example, for the jth feature, we create two mirrored features as standard Gaussian random variables. The scalar c j is chosen so as to minimize the dependence between X + j and X − j conditional on the remaining variables, X −j . To cope with nonlinear dependence, we introduce a kernel-based conditional independence measure and obtain c j by solving an optimization problem. Then, we construct two modifications of the network architectures based on the mirrored design. The proposed mirror statistics tend to take large positive values for important features and are symmetric around zero for null features, which enables us to control the selection error rate.

Controlled feature selection
In many modern applications, we are not only interested in fitting a model (with potentially many predictors) that can achieve a high prediction accuracy, but also keen in selecting relevant features with controlled false selection error rate. Mathematically, we consider a general model with an unknown link function F illustrating the connection between the response variable Y and p predictor variables X ∈ R p : Y = F (X)+ , where denotes the random noise. We define S 0 to be the set of "null" features: ∀ k ∈ S 0 , we have Y ⊥ ⊥ X j | X −j , where X −j = {X 1 , . . . , X p }\{X j } and define S 1 = {1, . . . , p}\S 0 as the set for relevant features. Feature selection is equivalent to recovering S 1 based on observations. When the estimated set S 1 is produced, the false discoveries can be denoted as S 1 ∩ S c 1 , where S 1 is the true set of important features. The false discovery proportion (FDP) is defined as and its expectation is called the false discovery rate (FDR), i.e., FDR = E[FDP].
Technically, controlled feature selection methods are different from conventional feature selection methods such as Lasso in that the former require additional efforts in assessing the selection uncertainty by means of estimating the FDP, which becomes more challenging in DNNs.

Gaussian mirror design for linear models
To discern relevant features from irrelevant ones, we intentionally perturb the input features via the following Gaussian mirror design: for each feature X j , we construct a mirrored pair (X + j , X − j ) = (X j + c j Z j , X j − c j Z j ) with Z j ∼ N (0, 1) being independent of all the X. We call (X + j , X − j , X −j ) the j-th mirror design. Regressing Y on the j-th mirror design, we obtain coefficients via the ordinary least squares (OLS) method. Let β + j and β − j be the estimated coefficients for X + j and X − j , respectively. The mirror statistics for the jth feature is defined as: For important features j ∈ S 1 , β + j tends to be similar to β − j , which helps cancel out the perturbation part leading to a large positive value of Figure 1: Histogram of mirror statistics: red represents important variables and blue stands for null variables M j . From another point of view, the large value of M j indicates that the importance of the jth features is more stable to the perturbation, i.e. β + j ≈ β − j . For null features j ∈ S 0 , the mirror statistics can be made symmetric about zero by choosing a proper c j as detailed in Theorem 2.1.
i.e., the partial correlation of X + j and X − j given X −j . For any j ∈ S 0 , if we set

3)
then we have P (M j < −t) = P (M j > t) for any t > 0.
A detailed proof of the theorem can be found in [31]. It shows that the symmetric property of the mirror statistics M j for null features is satisfied if we choose c j that minimizes the magnitude of the partial correlation I L j (c), which is equivalent to solving I L j (c) = 0 in linear models. In other words, the perturbation c j Z j makes X + j and X − j partially uncorrelated given X −j . We note that in linear models with OLS fitting, the above construct is equivalent to generating an independent noise feature c j Z and regress Y on (X, c j Z). The mirror statistic is simply the magnitude difference between the estimated coefficient for X j and that for c j Z. However, this simpler construct is no longer equivalent to the original design for nonlinear models or for highdimensional linear models where the OLS is not used for fitting. Also, the conditional independence motivation for choosing c j is lost in this simpler formulation. Figure 1 shows the distribution of the M j 's. It is seen that relevant features can be separated from the null ones fairly well and a consistent estimate of the false discoveries proportion based on the symmetric property can be obtained.
For a threshold t, we select features as Thus we obtain an estimate of the FDP as In practice, we choose a data adaptive threshold to control the FDR at a predefined level q. As shown in Theorem 4 of [31], under weak dependence assumption of M j s, we have as p goes to infinity, i.e., the FDR is asymptotically controlled under the predefined threshold q.

Neural Gaussian mirror
We here describe a model-free mirror design for controlling variable selection errors in DNNs. In linear models, a key step of our mirror design is to choose a proper perturbation level c j so as to annihilate the partial correlation between the mirror variables X + j and X − j . Since the partial correlation being zero for two random variables implies that they are conditional independence under the joint Gaussian assumption, for more complex models such as DNNs, we consider a general kernel-based measure of conditional dependence between X + j and X − j and choose c j by minimizing this measure. Then, we can construct the mirror statistics and establish a data-adaptive threshold to control the selection error rate in a similar way as in linear models. We call this whole procedure the neural Gaussian mirror (NGM).

Model-free mirror design
Without any specific model assumption, we construct the mirrored pair as , and c j is chosen as Here [I K j (c)] 2 is a kernel-based conditional independence measure between X + j and X − j given X −j . We give a detailed expression of [I K j (c)] 2 in Sections 3.2 and 3.3. Compared with (2.3), the kernel-based mirror design can incorporate nonlinear dependence effectively.

Decomposition of log-density function
For notational simplicity, we write (U, V, W ) = (X j + cZ, X j − cZ, X −j ) and denote their logtransferred joint density function as η, which belongs to a tensor product reproducing kernel where H U , H V and H W are marginal RKHSs and "⊗" denotes the tensor product of two vector spaces. We decompose η as The uniqueness of the decomposition is guaranteed by the probabilistic decomposition of H.
For simplicity, we use the Euclidean space as an example to illustrate the basic idea of tensor sum decomposition, which is also known as ANOVA decomposition in linear models. For example, for the d-dimensional Euclidean space, we let f be a vector and let f (x) be its x-th entry, is a probability vector (i.e., δ i ≥ 0, and δ i = 1). The tensor sum decomposition of the Euclidean space R d is where the first space is called the grand mean and the second space is called the main effect. Then, we construct the kernel for However, in RKHS of infinite dimension, the grand mean is not a single vector. Here, we set the average operator A l as is the kernel function in H l and the first equality is due to the reproducing property. E x K l x plays the same role as δ in Euclidean space. Then we have the tensor sum decomposition of marginal RKHS H l defined as 3) Following the same fashion, we call H l 0 as the grand mean space and H l 1 as the main effect space. Note that E x K l x is also known as the kernel mean embedding, which is well established in the statistics literature [7]. We show the kernel functions for H l 0 and H l 1 in Lemma A.3 (see SM for details).
Following [17], we apply the distributive law and have the decomposition of H as We show the kernel functions for each subspace in Lemma A.4 (see SM). Each component in (3.2) is the projection of η on the corresponding subspace in (3.4). Thus, the decomposition of the log-density function in (3.2) is unique.
We introduce the following lemma to establish a sufficient and necessary condition for the conditional independence of U and V given W based on the decomposition in (3.2).
Lemma 3.1 Assume that the log-joint-density function η of (U, V, W ) belongs to a tensor product RKHS. U and V are conditional independent if and only if η U,V + η U,V,W = 0.
Here η C is the function only of variables in the set C. The proof of this lemma is given in SM. Lemma 3.1 implies that the following two hypothesis testing problems are equivalent: and where a key advantage of (3.6) is that we are able to specify the function space of the log-transferred density under both null and alternative.

Kernel-based Conditional Dependence Measure
Let where P H 0 is a projection operator from H to H 0 . Using the reproducing property, we rewrite (3.7) as where K 1 is the kernel for H 110 ⊕H 111 . The score function 1 n K 1 t i is the first order approximation of the likelihood ratio functional. We define our kernel-based conditional measure as the squared norm of the score function of the likelihood ratio functional as which, by the reproducing property, can be expanded as The construction of [I K j (c)] 2 is related to the kernel conditional independence test [13], which generalizes the conditional covariance matrix to a conditional covariance operator in RKHS. However, the calculation of the norm of the conditional covariance operator involves the inverse of an n × n matrix, which is expensive to obtain when the sample size n is large. But we will show that our proposed score test statistics only involve matrix multiplications.
We introduce a matrix form of the squared norm of score function to facilitate the computation process. In (3.11) where H = I n − 1 n 11 T , I n is a n × n identity matrix and 1 n is a n × 1 vector of ones, and [A] ++ = n i=1 n j=1 A ij .The most popular kernel choices are Gaussian and polynomial kernels. With parallel computing, the computational complexity for using either is approximately linear in n.

Mirror Statistics in DNNs
In order to construct mirror statistics in DNNs, we introduce an importance measure for each input feature as a generalization of the connection weights method proposed by [10]. We consider a fully connected multi-layer perceptron (MLP) with k hidden layers denoted as (H where g t+1 is the activation function. We let t = 0 and t = k + 1 denote the input and output layers, respectively.
For a path l = (X j , H i k , y) through the network, We define its accumulated weight as Let Ω jn 1 ) be the weight vector connecting X j with the first layer. We define the feature importance of X j as where P j is the set consisting of all paths connecting X j to y with one node in each layer and ij ∈ R nt×n t+1 for i = 1, . . . , n t and j = 1, . . . , n t+1 . The gradient w.r.t. each input feature ∂y/∂X j in [18] is a modified version of L(X j ) (see the SM for details).
Based on (3.13), we define mirror statistics as and use (2.4) to estimate FDP. By setting a predefined error rate q, a data adaptive threshold is given by τ q = min t {t > 0 : F DP (t) ≤ q} and S 1 = {j : M j ≥ τ q } is the set of selected features.

Two forms of neural Gaussian mirror
The NGM can be realized in two ways: individual or simultaneous. As shown in Figure 2, the individual neural Gaussian mirror (INGM) is constructed by treating one feature a time. More specifically, for feature X j , we create X + j and X − j with c j minimizing [I K j (c)] 2 . Then, we set the input layer as (X 1 , . . . , X + j , X − j , . . . , X p ) and fit the MLP as shown in Figure 2. For any predefined error rate q, we select feature by Algorithm 1.

4:
Create mirrored pair(X + j , X − j ) = (X j + c j Z j , X j − c j Z j )
To increase the efficiency, we construct the simultaneous neural Gaussian mirror (SNGM), by mirroring all the p features at the same time as shown in Figure 3. The input layer include all mirrored pairs, i.e., (X + 1 , X − 1 , . . . , X + p , X − p ). The detailed algorithm is given in Algorithm 2.

Algorithm 2 SNGM: Feature selection with simultaneous-mirrored Pairs
Input: Fixed FDR level q, (x i , y i ), i = 1, ..., n with x i ∈ R p , y i ∈ R 1 Output: S 1 1: for i = 1 to p do 2: Generate Z j ∼ N (0, I n )

3:
Compute c j = arg min c [I K j (c)] 2

Screening
To further reduce the computational cost, we propose a screening step based on the rank of feature importance measure defined in (3.13). The screening procedure is inspired by the RANK method proposed by [11,22], which uses part of the data for estimating the precision matrix and subset selection and leaves the remaining data for controlled feature selection.
In the screening procedure, we first randomly select [n/3] samples (X (1) , y (1) ) to train a neural network. We calculate feature importance measure {L(X

Numerical simulations
We choose an MLP structure with two hidden layers, with N 1 = 20log(p) and N 2 = 10log(p) hidden nodes respectively.We compare the selection power and FDR of NGMs with the DeepPINK, and conduct experiments in both simulated and real-data settings with 20 repetitions. In the simulation studies, we consider data from both linear models and nonlinear models. In each setting we consider two covariance structures for X: the Toeplitz partial correlation structure and the constant partial correlation structure. The Toeplitz partial correlation structure has its precision matrix (i.e., inverse) as Ω = (ρ |i−j| ). The constant partial correlation structure has its precision matrix as Ω = (1 − ρ)I p + ρJ p where J p = 1 p 1 p 1 T p , where ρ = 0.5.

Linear models
First, we examine the performance of all methods in linear models: ∼ N (0, 1), for i = 1, . . . , n. We randomly set k = 30 elements in β to be nonzero and generated from N (0, (20 log(p)/n) 2 ) to mimic various signal strengths in real applications.
As shown in Figure 4, NGMs control the FDR at q = 0.1 and have a higher power than the Deep-PINK. In the constant partial correlation setting with partial correlation ρ = 0.5, DeepPINK shows a power loss when p is larger than 1500 since the minimum eigenvalue of the correlation matrix is 1/(1 + (p − 1)ρ), which approaches 0 as p → ∞ and makes important features highly correlated with their knockoff counterparts. By introducing a random perturbation, which reduces the correlation between mirror statistics, NGMs control FDR at the level of q = 0.1 and maintain power around 0.8 in both correlation structures.
In high-dimensional cases when p ≥ 2000, the NGMs, especially the SNGM and S-SNGM, are computationally more efficient than DeepPINK. In addition, we note that the difference in time consumption between DeepPINK and SNGM

Single-index models
We further test the performance of the proposed methods and DeepPINK in nonlinear cases such as single-index models. In our experiment, we choose three nonlinear link functions: f 1 (t) = t + sin(t), f 2 (t) = 0.5t 3 and f 3 (t) = 0.1t 5 . We try both high-dimensional cases with p = 2000 and low-dimensional cases with p = 500. As shown in Table 1, NGMs and DeepPINK all have desirable performances similar to that in linear  is dominated by its linear term. However, in the latter two cases where f is a polynomial, Table 1 shows that DeepPINK suffers a loss of power. In all the three cases, NGMs are capable of maintaining the power at a higher level than DeepPINK.
To study the influence of sparsity, we further consider an experiment with constant partial correlation ρ = 0.5 and varying the number of important features from 10 to 50 skipping by 10. For varying sparsity levels, as shown in Table  1, NGMs maintain high power with controlled error rate in both the low dimensional case with p = 500 and the high dimensional setting with p = 2000, whereas DeepPINK tends to loss power when the sparsity level is high.

Real data design: tomato dataset
We consider a panel of 292 tomato accessions in [4]. The panel includes breeding materials that are built and characterized with over 11,000 SNPs. Each SNP is coded as 0, 1 and 2 to denote the homozygous (major), heterozygous, and the other homozygous (minor) genotypes, respectively. To evaluate the performance of NGMs and Deep-PINK under non-Gaussian designs, we randomly select p ranging from 100 to 1000 SNPs to form our design matrix, and generate the response y i in the same way as with the simulated examples. This simulation is replicated 20 times indepen- dently, to which each method in consideration is applied.
We set the pre-defined FDR rate at q = 0.2 and compare the empirical FDR and power of NGMs with those of DeepPINK. As shown in table 2, the S-INGM method shows the lowest FDR and highest power when p = 100 and 500. When p =1000 and 2000, S-INGM is more powerful than DeepPINK but the its FDR inflates slightly. S-SNGM is computationally more efficient, but with slightly inferior performances, than S-INGM.
To provide a better visualization of the tradeoff between the error rate and power, we plot the ROC curve for the three methods. Specifically the ROC curve is plotted with true positive rate

Discussion
In this paper, we propose NGMs for feature importance assessment and controlled selection in neural network models, which is an important  aspect for model interpretation. Even in situations with no need of explicit feature selections, having a good rank of the predictors according to their influence on the output can be very helpful for practitioners to prioritize their follow-up work. We emphasize that our method does not require any distributional assumption on X, thus is widely applicable to a broad class of neural network models including those with discrete or categorical features. In addition, the mirror design can be generated to convolutional neural networks (CNN) to measure the stability of the filters. For example, for each original CNN slice in the input tensor, we can create its corresponding mirrored slice by adding and subtracting a slice of random perturbations, and contrast the influences of the original and mirrored slices. detailed study in this direction is deferred to a future work. The data used in the real example is available in (ftp.solgenomics. net/manuscripts/Bauchet_2016/).

Supplementary materials A Test statistics for conditional independence
To simplify the notations, let Proof Suppose f is the joint density function for (U, V, W ) = (X + j , X − j , X −j ), then it can be decomposed into 2 3 = 8 factors: where η represents log-transformed density and η A (α) means the part as the function only of α. Then To assess the dependence between (U, V ) conditioned on W , we derive the conditional density as where C(w) denotes the denominator as the marginal density of W . Therefore U ⊥ ⊥ V |W if and only if any interaction of U and V is constants, i.e. η U,V + η U,V,W = 0.

B Importance measurement in DNN
In order to construct mirror statistics in DNNs, we first adopt an importance measure for each input feature. We consider a fully connected multi-layer perceptron (MLP) with k hidden layers denoted as (H (t) 1 , . . . , H (t) nt ) for t = 1, . . . , k. The output of the (t + 1)-th layer can be described as where g t+1 is the activation function. Note that t = 0 denotes the input layer and t = k + 1 denotes the output layer.
We define δ(l(X j , H as the accumulated weight of the path l. Let jn 1 ) as the weight vector connecting X j with the first layer. We define the feature importance of jth feature as where P j is the set consisting of all paths connecting X j to y with one node in each layer and ij ∈ R nt×n t+1 for i = 1, . . . , n t and j = 1, . . . , n t+1 .
For t = 1, . . . , k, define diagonal matrices Then, the gradient w.r.t. X j can be written as which is a weighted version of L(X j ).

C Perturbation with linear kernel
Suppose n observations (x i , y i ), i = 1, 2, . . . , n are drawn from the model where φ(·) is a feature map and the kernel function k(·, ·) is define as In the classical linear model y = x β + , φ is chosen to be the identity with a constant interception and thus k(·, ·) = ·, · + const.
To specify the test statistics when adopting linear kernels We take gradient of G(c) with respect to c and due to the community of operators: Therefore, ∇ s G(s) = 0 leads to the solution as follow: with A(X j , Z j ) = (Z j X j ) • (Z j X j ) + (X j Z j ) • (X j Z j ). And c * j is the minimum point by computing the second-order derivatives w.r.t. c.
For the space span{W }, let P W be the projection operator, then c j in the Gaussian Mirror can be written as We should note that, without the term W W , Therefore, comparing the form of C.4 and C.5, c * j and c GM j are of the same scale which is validated in the simulation and they are equivalent with orthogonal designs. ∼ N (0, 1), for i = 1, . . . , n. We randomly set k = 30 elements in β to be nonzero and generated from N (0, (20 log(p)/n) 2 ) to mimic various signal strengths in real applications. We set n = 1000 and study both the low-dimensional case with p = 500 and the high-dimensional case with p = 2000 respectively.
As shown in the Table 3, S-SNGM and S-INGM have exact FDR control under q = 0.1 and meanwhile have higher power than the DeepPINK. Besides, SNGM controls FDR under q = 0.1 except for the case with p = 2000, ρ = 0.8 that is diffcult for mirror construction. In the highdimensional setting with high partial correlation, DeepPINK undergoes an obvious power loss since the minimum eigenvalue of correlation matrix is 1/(1 + (p − 1)ρ) ≈ 0, which imposes great obstacles for the Knockoff construction, but NGM with screening procedure is more stable in power maintenance.