Identifying causal relationships is a profound and hard problem pervading the sciences. While randomized controlled intervention studies are considered the gold standard, they are in many cases ruled out by financial or ethical concerns [26, 36] . In order to improve understanding of systems and help design relevant interventions, the subset of causes which have a direct effect (direct causes / direct causal parents) need to be identified.
Let us consider the setup described in Figure 1 , corresponding to a linear system only for the response Y ,
where U is an independent noise variable with zero mean, θ, X ∈ R d , Y ∈ R and ·, · denotes the inner product. We investigate how to find the direct causes of Y among a high-dimensional vector of covariates X, where the covariates can have arbitrary non-linear relationships between them. From our formulation, a given entry of θ should be non-zero if and only if the variable corresponding to that particular coefficient is a direct causal parent [28] , e.g., X 1 and X 2 in Figure 1 . We restrict ourselves to the setting of linear direct causes of Y (LDC) and no feature descending from Y (NFD). LDC is justified as an approximation when the effects of each causal feature are weak such the possibly non-linear effects can be linearized; NFD is justified in some applications where we can exclude any influence of Y on a covariate. This is for example the case when X are genetic factors, and Y is a particular trait/phenotype. Figure 1 : Graphical representation of Causal Feature Selection in our setting. Out of variables {X 1 , · · · , X 11 }, only X 1 and X 2 are direct causal parents of Y such that Y = θ 1 X 1 + θ 2 X 2 + U , U being an independent zero-mean noise. We propose an approach to find X 1 and X 2 in this paper under assumptions discussed in the text. An example of this setup in the real-world is finding genes which directly cause a phenotype.
While applicable to full graph discovery rather than the simplified problem of finding causal parents, state of the art methods for causal discovery often rely on strong assumptions or the availability of interventional data, and have prohibitive computational cost as explained in section 1.1 in more detail. In addition to and despite their strong assumptions, causal discovery methods may perform worse than simple regression baselines [16, 18, 42] .
While plain regression techniques have appealing computational cost, they come without guarantees. When using unregularized least-square regression to estimate θ, there can be infinitely many possible choices for θ recovered with equivalent prediction accuracy for regressing Y , in case of over-parametrized models. However, none of these choices provide any information about the features which, when intervened upon, directly cause the output variable Y . On the other hand, when using a regularized method such as Lasso, a critical issue is the bias induced by regularization [20] .
When knowing the distinction between covariates and direct causes, Double ML approaches [7] have shown promising bias compensation results in the context of high dimensional observed confounding of a single variable. In the present paper, we generalize them to the problem of finding direct causes.
Our key contributions are:
• We show that under the assumption that no feature of X is a child of Y , the Double ML [8] principle can be applied in an iterative and parallel way to find the subset of direct causes, even in the challenging case of only having observational data. • Our approach has low computational complexity requirement (O(d) in terms of number of variables), outperforming the usual O(2 d ). • Our method provides asymptotic guarantees that the set can be recovered exactly. Importantly, this result is valid without additional assumptions, and especially neither requires linear interactions among the covariates, nor faithfulness. • Extensive experimental results demonstrate state-of-the-art performance of our method. Our approach significantly outperforms all other methods, especially in the case of non-linear interactions between covariates, despite relying only on linear projection.
The question of finding direct causal parents is also addressed in the literature as mediation analysis [3, 15, 34] . Several principled approaches have been proposed (relying for instance on Instrumental Variables (IVs)) [1, 2, 4] to test for a single direct effect in the context of specific causal graphs.
Extensions of the IV based approach to generalized IVs based approaches [5, 39] are the closest known result to discover direct causal parents. However, no algorithm is provided in [5] to identify the instrumental set. Subsequently, an algorithm is provided in [39] for discovering the instrumental set in the simple setting where all the interactions are linear and the graph is acyclic. In contrast, our method allows non-linear interaction amongst the variables.
Several other works have also tried to address the problem of discovering causal features. The authors review work on causal feature selection in [14] . More recent papers on causal feature selection have appeared since [6, 25, 40] , but none of those claims to recover all the direct causal parents asymptotically or non-asymptotically as we do in our case.
There has been another line of works on inferring causal relationships from observational data most of which require strong assumptions, such as faithfulness [21, 26, 36] . Classical approaches along these lines include the PC-algorithm [36] , which can only reconstruct the network up to a Markov equivalence class. Another approach is to restrict the class of interactions among the covariates and the functional form of the signal-noise mixing (typically considered additive) or the distribution (e.g., non-Gaussianity) to achieve identifiability (see [17, 30] ); this includes linear approaches like LiNGAM [33] and nonlinear generalizations with additive noise [29] . For a recent review of the empirical performance of structure learning algorithms and a detailed description of causal discovery methods, we refer to [16] . Recently, there have been several attempts at solving the problem of causal inference by exploiting the invariance of a prediction under a causal model given different experimental settings [13, 27] . The computational cost to run both of the algorithms is exponential in the number of variables, when aiming to discover the full causal graph.
Our method mainly takes inspiration from Debiased/Double ML method [7] which utilizes the concept of orthogonalization to overcome the bias introduced due to regularization. We will discuss this in detail in the next section. Considering specific example, the Lasso suffers from the fact that the estimated coefficients are shrunk towards zero which is undesirable [37] . To overcome this limitation, an approach of debiasing the Lasso was proposed in several papers [19, 20, 41] . However, unlike our approach, Debiased Lasso methods don't recover all the non-zero coefficients of the parameter vector θ under generic assumptions we have in the paper.
Before describing the proposed method, we quickly discuss Double ML and Neyman orthogonality in the next section which will be helpful in building the theoretical framework for our method.
Given a fixed set of policy variables D and control variables X, an unbiased estimator of the parameter θ 0 for the partial regression model in Equation (2) can be obtained via the orthogonalization approach as in [7] . 
where Y is the outcome variable, U, V are disturbances and g 0 , m 0 : R d → R are non-linear functions. Orthogonalization is obtained via the use of "Neyman Orthogonality Condition" which we describe below. Figure 2 : A generic example of identification of a causal effect θ in the presence of causal and anti-causal interactions between the causal predictor and other putative parents, and possibly arbitrary nonlinear and non-invertible assignments for all nodes except Y (see Proposition 2) We have X −k = Z 1 ∪ Z 2 .
Neyman Orthogonality Condition: the traditional estimator of θ 0 in Equation (2) can be simply obtained by finding the zero of the empirical average of a score function ψ such that ψ(W ; θ, g) = D (Y − Dθ − g(X)). However, the estimation of θ 0 is sensitive to the bias in the estimation of the function g. Neyman [24] proposed an orthogonalization approach to get the estimate for θ which is more robust to the bias in the estimation of nuisance parameter (m 0 , g 0 ). For a moment, if we assume that the true nuisance parameter is η 0 (which represents m 0 and g 0 in Equation (2)) then the orthogonalized "score" function ψ should satisfy the property that the Gateaux derivative operator with respect to η vanishes when evaluated at the true parameter values:
Orthogonalized or Double/Debiased ML estimatorθ 0 solves 1 n n i=1 ψ(W ;θ 0 ,η 0 ) = 0, whereη 0 is the estimator of η 0 and ψ satisfies condition in Equation (3) . For the partially linear model discussed in Equation (2), the score function ψ is, ψ(W ; θ, η) = (Y − Dθ − g(X))(D − m(X)) (4) with η = (m, g).
From Double ML to Causal Discovery: The distinction between policy variables and confounding variables is not always known in advance, which motivates us to consider the more general setting of causal discovery. To this end, we consider a set of variables X = {X 1 , X 2 , · · · X d } which includes direct causal parents of the outcome variable Y as well as other variables. We also reiterate our assumption that the relationship between outcome variable and direct causal parents of the outcome variable is linear. The relationship among other variables can be nonlinear. We now provide a general approach to scanning putative direct causes scaling "linearly" with their number, based on the application of a statistical test and a parameter estimation based on Debiased machine learning method. We describe first the algorithm and then provide theoretical support for its performance.
We provide pseudo-code for our proposed method in Algorithm 1. Intuitively, the idea is to do a one-vs-rest split for each variable in turn, and try to estimate the link between that particular variable and the outcome variable using Debiased approach. To do so, we decompose eq. 1 to single out a variable D = X k as policy variable and take the remaining variables Z = X −k = X\X k as multidimensional control variables, and run Double ML estimation assuming the partial regression model presented in Sec. 2.1, which now takes the form
The step-wise description of our estimation algorithm goes as follows: (a) Select one of the variables X i to estimate its (hypothetical) linear causal effect θ on Y . (b) Set all of the other variables X −i as the set of possible confounders. (c) Use the Double ML approach to estimate the parameter θ i.e. the causal effect of X i on Y . (d) If the variable X i is not a causal parent, the distribution of the conditional covariance χ i (Proposition 3) is a Gaussian centered around zero. We use a simple normality test for χ i to select or discard X i as one of the direct causal parents of Y .
We iteratively repeat the procedure on each of the variables until completion. Pseudo-code for the entire procedure is given below in Algorithm 1.
Note that only for (5) this corresponds to the Double ML approach by estimating θ k . For all other cases (for instance D not corresponding to direct causes) it is a none trivial result to show that one obtains a non-zero parameter. Remarks on Algorithm 1: X
i is a vector which corresponds to the samples chosen in the k th subsampling procedure, X In general the subscript i represents the estimation for the i th variable and super-script k represents the k th subsampling procedure. K represents the set obtained after sample splitting. m [\k] i are (possibly nonlinear) parametric functions fitted using (1 st , . . . , k − 1 th , k + 1 th , . . . , K th ) subsamples.
Now we describe the execution of our algorithm for a simple graph with 3 nodes. Let us consider the following linear structural equation model as an example of our general formulation:
Y := θ 1 X 1 + θ 2 X 2 + ε 3 , X 2 := a 12 X 1 + ε 2 , and X 1 := ε 1 .
Example 1. Let us consider the system whose structural equation model is given in (6) . If ε 1 , ε 2 and ε 3 are independent uncorrelated noise terms with zero mean, then Algorithm 1 will recover the coefficients θ 1 and θ 2 .
A detailed proof is given in the Appendix. While the estimation of the parameter θ 1 is in line with the assumed partial regression model of eq. (5), the estimation of θ 2 does not follow the same. However, it can be seen from the proof that θ 2 can also be estimated from the orthogonal score (Equation (4)).
Algorithm 1 Efficient Causal Structure Search 1: Input: response Y ∈ R N , covariates X ∈ R N ×d , significance level α, number of partitions K. 2: Split N observations into K-fold random partitions, I k for k = 1, 2 . . . , K, each having n = N/K observations. 3: for i = 1, . . . , d do 4: for Subsample k ∈ [K] do 5: D k ← X [k] i and Z k ← X Fit g
ij (Z kj ) 10 :χ
ij (Z kj ) −D kj g [\k] ij (Z kj ) + m [\k] ij (Z kj ) g [\k] ij (Z kj ) +V [k] ij D kj 11:
j m [\k] ij (Z kj ) −D kj g [\k] ij (Z kj ) + m [\k] ij (Z kj )g [\k] ij (Z kj ) +V [k] ij D kj −χ [ We now show that this result holds for a more general graph structure given in figure 2, allowing also for non-linear interactions among features. Proposition 2. Assume the partially linear Gaussian model of Fig. 2 , denote X −k = [Z 1 , Z 2 ] the control variables, γ = (γ 1 , γ 1 , γ 12 ) the parameter vector of the (possibly non-linear) assignments between putative parents of Y , and β = (β 1 , β 2 ), the vector of causal coefficients for encoding linear effects of X −k on outcome Y . Then, independently of the γ parameters and of the functional form of the associated assignments between parents of Y , the score
with r XX −k = E[X k X −k ]E[X −k X −k ] −1 , follows the Neyman orthogonality condition for the estimation of θ with nuisance parameters η = (β, γ) which reads
Comparing the score in Equation (7) with the score in Equation (4), there are two takeaways from Proposition 2: (i) the orthogonality condition remains invariant irrespective of the causal direction between X k and Z, and (ii) the second term in Equation (15) replaces function m by the (unbiased) linear regression estimator for modelling all the relations; given that the relation between Z and Y is linear, even if relationships between Z and X k are non-linear. Combining with the Double ML theoretical results [7] , this suggests that regularized predictors based on Lasso or ridge regression are tools of choice for fitting functions (m, g).
Now that we have illustrated and justified the fitting procedure of our algorithm, we provide a theoretically justified statistical decision criterion for the direct causes after the model has been fitted. Consider (Y, X), Y ∈ R, X ∈ R d , satisfying
E(Y 2 ) < ∞, E(U 2 ) < ∞, E(U ) = 0, E(U | X j ) = 0, ∀j, and E( X 2 2 ) < ∞,
where U is noise variable and X −j represents all the variables except X j . The assumptions made with the above formulation are standard in the orthogonal machine learning literature [31, 35, 8] .
Let us define the quantity
for j ∈ {1, · · · , d}, which is the expected conditional covariance of X j given X −j . Proposition 3. Let P A Y = {j ∈ {1, . . . , p} : θ j = 0} . For each j ∈ {1, . . . , p} let X −j be the vector equals to X but excluding coordinate j and define θ −j similarly. Define for j ∈ {1, . . . , p}
which also has the double robustness property ( [8, 31] ) then under the conditions given in Equa-
There are two main implications of the results provided in Proposition 3. (i) χ j is non-zero only for direct causal parents of the outcome variable and χ j has double robustness property as shown in [31, 35, 8] , hence one can use regularized methods like ridge regression or Lasso to estimate the function m. Afterwards, one can perform statistical tests on top of it to decide between zero or non-zero test. (ii) In line with the above orthogonal score results, we see that this quantity can be estimated using linear (unbiased) regression to fit the function m, although interactions between features may be non-linear. Next, we discuss the variance of our estimator so that later a statistical test can be used to differentiate between zero and non-zero test.
Variance of Empirical Estimates of χ j :
Randomly split the data in two halves, say D n1 and D n2 . Take j ∈ {1, . . . , p}. For k = 1 let k = 2, for k = 2 let k = 1. For k = 1, 2, compute estimates of
can be considered as regularized regression problems. We use Lasso as the estimator for conditional expectation (eq. (17)) in the experiments. Now, we compute the empirical estimates of χ j . Theorem 1 of [35] provides conditions under which (see also [8] ), when the estimators
are Lasso-type regularized linear regressions, it holds that asymptotically χ j ≈ N χ j ,
In this case, the test that rejects χ j = 0 when | χ j | ≥ 1.96
σj √ n will have approximately 5% level. The probability of rejecting the null when it is false is
Expressions for χ j and σ j are provided in the Appendix A.2. In order to account for multiple testing, we use Bonferroni correction. Comments about Estimator: In this paper, we use Lasso for the nuisance parameter estimation as the variance of the conditional covariance is known [35] . One can also use other estimators instead, assuming one obtains a good enough estimate of the nuisance parameter (upto N −1/4neighbourhood [7] ) with the right variance term, which is beyond the scope of this paper.
In this section, we perform extensive empirical evaluation for our method.
For every combination of number of nodes (#nodes), sparsity (p s ), noise level (σ 2 ), number of observations (z), and non-linear probability (p n ) (see Table C For each DAG, z number of samples are generated. We provide more details about the parameters (#nodes, p s , p n and z) and data generation process in Appendix C. For future benchmarking, the generated files with the code will be made available later.
The baselines we compare our method against are: LINGAM [33], order-independent PC [11] , rankPC, MMHC [38], GES [9] , rankGES, ARGES (adaptively restricted GES [23]), rankARGES, FCI+ [10] and Lasso which are suitable for observational data. CompareCausalNetworks Package 1 is used to run the baselines methods. We use 10-fold cross validation to choose the parameters of all approaches. Recall, Fall-out, Critical Success Index, Accuracy, F1 Score, and Matthews correlation coefficient [22] are considered as metric for the evaluation. These metrics are described in Appendix C.
Regression Technique and Hyper-parameters: We use Lasso as the estimator of conditional expectation for our method because the variance bound for χ j with Lasso type estimator of conditional expectation (Equation (17)) is provided in Equation (16). Further, using more splits than 2 splits in the experiment increases the performance of parameter estimation. Plots for parameter estimation are provided in Appendix F. It can be observed that the accuracies of the competing baselines significantly drop with increasing noise level, while our method is more robust to it. More plots are given in Appendix G and Appendix H for several other combinations of varying parameters in the simulation. We also extensively compare all the metrics (Recall, Fall-out, Critical Success Index, Accuracy, F1 Score and Matthews correlation coefficient) for all the methods in Appendix E. According to these metrics, our approach performs better than baselines in most of the cases regardless of the set of parameters used for generating data. Our method shows in particular stability in performance w.r.t. the number of nodes (Table  E. 1), partially non-linear relationships (Table E. 3), sparsity (Table E. 2), number of observations (Table E.4) , and noise level (Table E. 5) . We also show the plot of parameter estimation for direct causal parents vs. non-causal parents in Figure 4 . In the plots and tables, we denote our approach as CORTH Features. Figure 5 shows the runtime of the method in seconds as a function of the graph's size. Notice that the runtime of our algorithm in the log-log plot is roughly linear. Since we used 5000 observations, any additional overhead is coming from cross-validation. In these simulations sparsity, number of observations, nonlinaer prob., and noise level are set to 0.3, 5000, 0, and 1 respectively.
We also apply our algorithm on a real world data sets and evaluate the performance of our algorithm on a recent COVID-19 Dataset [12] where the task is to predict COVID-19 cases (confirmed using RT-PCR) amongst suspected ones. For an existing and extensive analysis of the dataset with predictive methods we refer to [32] . We apply our algorithm to discover the features which directly cause the diagnosed infection. We found that the following were the most common causes across different runs of our approach: Patient age quantile, Arterial Lactic Acid, Promyelocytes, and Base excess venous blood gas analysis. Lacking medical ground truth, we report these not as a corroboration of our approach, but rather as a potential contribution to causal discovery in this challening problem. It is encouraging that some of these variables are consistent with other studies [32] . Details on data preprocessing and more results are available in Appendix D.
A recent empirical evaluation of different causal discovery methods highlighted the desirability of more efficient algorithms [16] . In the present work, we provide identifiability results for the set of direct causal parents including the case of partially nonlinear models, as well as a highly efficient algorithm that scales linearly in the number of variables and exhibits state-of-the-art performance across extensive experiments. Our approach builds on the Double ML method of [7] , however, we stress that it goes beyond this since the distinction between policy variables and confounding variables is not always known in advance. Whilst not amounting to full causal graph discovery, identification of causal parents is of major interest in real-world applications, e.g., when assaying the causal influence of genes on the phenotype. A natural direction worth exploring is to extend this approach for discovering direct causal parents in the case when non-linear relation exists between output variable and its direct causal parents.
[ 
In this section, we use a generic example shown in Figure 2 which we show again in Figure A .1 to illustrate the role of interactions between the covariates on the proposed causal discovery algorithm.
The estimator discussed can simply be derived from the Neyman orthogonality condition. We now provide the below the proof for Proposition 2. For the sake of completeness, we also rewrite the statement of the proposition again. Proposition (Restatement of Proposition 2). Assume the partially linear Gaussian model of Fig. A.1 , denote X −k = [Z 1 , Z 2 ] the control variables, γ = (γ 1 , γ 1 , γ 12 ) the parameter vector of the (possibly non-linear) assignments between putative parents of Y , and β = (β 1 , β 2 ) the vector of causal coefficients for encoding linear effects of X −k on outcome Y . Then, independently from the γ parameters and of the functional form of the associated assignments between parents of Y , the score
with
follows the Neyman orthogonality condition for the estimation of θ with nuisance parameters η = (β, γ) which reads
Proof of Proposition 2. Using the Markov factorization
due to linearity and gaussianity of the assignement of Y , we obtain a negative log likelihood of the form (up to additive constants)
where f stands for the negative log likelihood of the second factor. Following [7] [eq. (2.7)], this leads to the Neyman orthogonal score
Following eq. (2.8) of the same paper, we derive the expression of µ as
Reintroducing µ in the expression of ψ leads to the result.
Proposition (Restatement of Proposition 3). Let P A Y = {j ∈ {1, . . . , p} : θ j = 0} . For each j ∈ {1, . . . , p} let X −j be the vector equals to X but excluding coordinate j and define θ −j similarly. Define for j ∈ {1, . . . , p}
then under the conditions given in Equations (8)-(10)
Proof. Take j ∈ P A Y . Then, from (8)
Thus
For c), we rewrite
Let G the sub-sigma algebra generated by X −j , under our assumptions,
is orthogonal to any elements of L 2 (Ω, G). Noticing that r XX −k X −j − E(X j | X −j ) is an element of L 2 (Ω, G), the second right-hand side term of the above equation vanishes and we get the result.
Suppose we have n i.i.d. observations D n = {(X i , Y i ), i = 1 . . . , n}. Randomly split the data in two halves, say D n1 and D n2 . Take j ∈ {1, . . . , p}. For k = 1 let k = 2, for k = 2 let k = 1. For k = 1, 2, compute estimates of E k (Y | X −j ) and E k (X j | X −j ) using the data in sample k.
Computing E k (Y | X −j ) and E k (X j | X −j ) can be considered as regularized regression problems.
We use Lasso as the estimator for conditional expectation (eq. (17)) in the experiments. Now, we compute the empirical estimates of χ j . Let,
P nk here denotes empirical average and σ k j denotes empirical variance of χ j . Finally, let
Theorem 1 of [? ] provides conditions under which (see also [8] ), when the estimators
are Lasso-type regularized linear regressions, it holds that asymptotically χ j ≈ N χ j ,
In this case, the test that rejects χ j = 0 when | χ j | ≥ 1.96
σj √ n will have approximately 5% level. The probability of rejecting the null when it is false is
The result discussed in Proposition 2 is not directly intuitive. In simple words, there are two takeaways from Proposition 2: (i) the orthogonality condition remains invariant irrespective of the causal direction between X k and Z, and (ii) the second term in equation 15 suggests to use a linear estimator for modeling all the relations, given that the relation between Z and Y is linear.
To generate more intuition, we provide a few examples. Let us go back again to the three variable interaction assuming the following structural equation model:
where f is a nonlinear function and ε 1 , ε 2 and ε 3 are zero mean Gaussian noises.
• Consider the case when f (x) = x 2 . The goal is to estimate the parameter θ 1 which we call θ 1 . We follow the standard double ML procedure assuming policy variable X 1 and control X 2 , although the ground truth causal dependency X 1 → X 2 in contradiction with such setting (see equation 2). The estimate of θ 2 following the double ML procedure, which we callθ 2 =
. Similarly, we want to estimate X 1 = αX 2 + η from which we get,
It is easy to see that E[X 1 X 2 ] = E[X 3 1 ] = 0. Hence, α = 0 and it is easy to seê θ 1 = θ 1 .
• Consider now the more general case where f is any nonlinear function. As in the previously discussed example, the goal is to estimate θ 1 . We haveθ 2 ] . We substitute these estimates into the orthogonality condition (15):
From the above two examples it is clear that even though the internal relations between the variables are nonlinear, all we need is an unbiased linear estimate to estimate the causal parameter in a directed acyclic graph.
For every combination of number of nodes (#nodes), sparsity (p s ), noise level (σ 2 ), number of observation (z), and non-linear probability (p n ) (look at (19)) with variance σ 2 starting from root of the DAG. For future benchmarking, the generated files will be made available with the code later on.
We generate DAGs (Direct Acyclic Graphs) in multiple steps: i) a random permutation of nodes is chosen as a topological order of a DAG. ii) Based on this order, directed edges are added to this DAG from each node to its followers with a certain probability p s (sparsity). iii) For each observation, values are assigned to nodes according to the topological order of the DAG in such a way that each node's value is determined by summing over transformations (linear or nonlinear with a certain nonlinear probability p n ) of values of its direct causes with the addition of Gaussian distributed noise. The non-linear transformation used is α tanh(βx), with α = 0.5 and β = 1.5. If the set of parents for the node X is denoted as P A X as before then value assignment for a node X is as follow:
where ε ∼ N (0, σ 2 ) in which σ 2 represents noise level. ι (X) is an indicator functions which decides between linear or non-linear contribution of X in X . We decide the value of ι (p n ) by generating a binary randon number which is 1 with probablity p n and 0 with probability 1 − p n . The value of θ is set to 2 for the small DAGs (number of nodes equal to 5 or 10) and 0.5 for large DAGs (number of nodes equal to 20 or 50) due to the value exploitation that might happen in large graphs.
We vary and investigate the effect of non-linear relationships, the number of nodes, number of observations, effect of sparsity and noise level while simulating the data. We summarize the factors in the data generation in table C.1.
Let the total number of true positives, false positives, true negatives ,and false negatives denoted by TP, FP, TN, and FN, we evaluate our method using following metrics:
• Recall (true positive rate):
T P R = T P T P + F N • Fall-out (false positive rate):
F P R = F P F P + T N • Critical Success Index (CSI): also known as Threat Score.
T P T P + F N + F P • Accuracy: .
In some rare cases, we encountered zero-divided-by-zero and divided-by-zero cases for some of these metrics. In these situations, scores are reported 1 and 0 respectively while Fall-out is reported 0 and 1.
The preprocessing stage for this dataset is the same as [32] except that, for each target variable upsampling is used to resolve data imbalance. Compared to other methods, our approach is stable wrt. noise level, sparsity and even partially non-linear relationships, significantly outperforming the majority of baselines.