Year: 2022

Abstract:
The ability of discovering feature interactions in a black-box model is vital to explainable deep learning. We propose a principled, global interaction detection method by casting our target as a multi-arm bandits problem and solving it swiftly with the UCB algorithm. This adaptive method is free of ad-hoc assumptions and among the cutting-edge methods with outstanding detection accuracy and stability. Based on the detection outcome, a lightweight and interpretable deep learning model (called ParaACE) is further built using the alternating conditional expectation (ACE) method. Our proposed ParaACE improves the prediction performance by 26% and reduces the model size by 100+ times as compared to its Teacher model over various datasets. Furthermore, we show the great potential of our method for scientific discovery through interpreting various real datasets in the economics and smart medicine sectors. The code is available at https://github.com/zhangtj1996/ParaACE.

heading: 1 INTRODUCTION
text: Explainable machine learning is an active research field that focuses on providing interpretable models, transparent explanations, and confident decisions to practical AI systems. Investigating feature interaction is vital to model interpretability. Interaction detection should be able to reveal which subset of features influence the output jointly, and what the corresponding nonlinear transformation is. We aim to find the underlying interactions from data, such that we can interpret the model properly. To go one step further, we hope that new models can be built with the aid of the detected interaction knowledge economically to avoid heavy parameterization.
In this paper, we first restrict ourselves to interaction detection. A novel method is derived directly from the most acknowledged definition of feature interaction (Friedman et al., 2008). We further apply the obtained interaction knowledge to design a transparent and refined neural network (NN). This approach fits well into the data science life cycle (Yu & Kumbier, 2019), which was applied, for instance in Tsang et al. (2020a), successfully for interpretable recommender system design.
The main contributions of this paper include: 1) a fast and principled interaction detection method, 2) a lightweight and interpretable neural network model that can surpass its Teacher, 3) thorough theoretical analysis and performance evaluations with real datasets. More details are given below.
1. We propose a generic interaction detection method based on a global statistical metric, namely the expected Hessian,Hij := Ex [ ∂2F (x) ∂xi∂xj ] . Notably, our method is model-agnostic
and applicable to any pre-trained learning model, F (x), being for instance, a deep neural network model or a tree model. It is also flexible to use for multi-way interaction detection.
2. To speed up the detection process, we evaluate the expected Hessian via adaptive sampling using the Upper Confidence Bound (UCB) algorithm (Lai & Robbins, 1985), which can significantly reduce the computational complexity. Besides, thorough analysis of the proposed interaction detection method are conducted.
3. Using the detected interaction pairs, we further design a compressed but interpretable Student model which can surpass its Teacher by 26% in terms of data fitting performance
averaged over various datasets. The compressed model reduces its size over 100 times compared to the baseline fully-connected, over-parameterized neural network (OverparaFC).
4. We demonstrate the linkages between our compressed model and the classic alternating conditional expectation (ACE) model (Breiman & Friedman, 1985).
5. We conduct large-scale performance evaluations and further explain the obtained model interpretability with some real datasets.
The remainder of this paper is organized as follows. Section 2 introduces all related works. Our proposed interaction detection method is introduced in Section 3. By exploiting the detection outcome, a new variant of lightweight and interpretable deep learning model is introduced in Section 4. Experimental results are given in Section 5. Finally, we conclude the paper in Section 6.

heading: 2 RELATED WORKS
text: Interaction Detection: Early works adopt pure statistics for detecting feature interactions, and representatives include ANOVA and GUIDE (Wonnacott & Wonnacott, 1990; Fisher, 1992; Loh, 2002). These works have motivated a plethora of new methods with the aim to enhance the detection accuracy and/or efficiency. The first class of methods centered around the GA2M and tree models, see for instance Lou et al. (2013); Sorokina et al. (2008); Friedman et al. (2008); Lundberg et al. (2020). The second class of methods were built on the so-called factorization machines (Rendle, 2010) as well as its new variants (Xiao et al., 2017; Song et al., 2019) with attention mechanism. The third class exploits the most recent advances in deep learning, including the Neural Interaction Detection (NID) and some new variants (Tsang et al., 2018a;b; Cui et al., 2019): Persistence Interaction Detection (PID) (Liu et al., 2020), Integrated Hessians (IH) (Janizek et al., 2020), Shapley interaction (Zhang et al., 2020; Sundararajan et al., 2020), etc. Although we have witnessed well improved interaction detection performance for many datasets over the decades, the above methods still lead to inconsistent detection results for some other datasets. The reasons are twofold. Firstly, the interaction strength is empirically defined, for instance, NID method computes the interaction strength via summarizing the neural network weights. Secondly, a specific deep learning model is required, for instance, an `1-regularized ReLU network is required by the latest NID and PID methods to maintain high accuracy. In contrast, our proposed method is derived directly from the definition of feature interaction (Friedman et al., 2008) and moreover is not confined to any specific learning model.
Model Interpretability: There are two categories of approaches to address model interpretability, namely the transparency-based and post-hoc approaches (Došilović et al., 2018). Transparency-based approaches require the model itself to be simple and interpretable, like linear models, decision trees, etc. One can directly read off the interpretations from their coefficients or decision rules. But often, they are less accurate due to limited representation power. In contrast, post-hoc approaches extract useful information from a pre-trained model, which is often complex and hard to interpret. Well-known methods such as LIME (Ribeiro et al., 2016) and SHAP (Lundberg & Lee, 2017) fall in this category, but they did not take feature interaction into account. Recently proposed symbolic metamodel (Alaa & van der Schaar, 2019) captures nonlinear interactions by approximating the black-box model with explicit symbolic expressions. There are also some other post-hoc approaches based on sensitivity analysis (Cortez & Embrechts, 2013), which return a quantification of feature importance νi = E[(
∂F (x) ∂xi )2] (Kucherenko et al., 2009) and interactions νij = E[|∂ 2F (x) ∂xi∂xj
|2] (Roustant et al., 2014) by input perturbation. Our work aims to combine the strengths of the two categories. Concretely, we first extract the interaction knowledge by a post-hoc method, and then build a transparent and interpretable learning model as illustrated in Figure 3.
Model Compression and Knowledge Distillation (KD): Model compression (Bucilua et al., 2006) aims to learn a small compressed model (Student) from a large complex model (Teacher) with augmented training data produced by the Teacher. Compared to the Teacher, the Student can make similar or even better predictions. Knowledge distillation (Hinton et al., 2015) mainly deals with multi-class classification problems and extracts “valuable information that defines a rich similarity structure over the data”. Our work introduces a novel viewpoint of knowledge (the interacted relationships) and targets a lightweight but more accurate Student model.

heading: 3 PROPOSED INTERACTION DETECTION METHOD
text: Before diving into in-depth discussions of interaction detection, we need to formally define what feature interaction is. The textbook definition according to Friedman et al. (2008) is given below. Definition 3.1 (Friedman & Popescu 2008). A function F : Rp → R is said to exhibit an interaction between two of its variables xi and xj if the difference in the value of F (x) as a result of changing the value of xi depends on the value of xj .
Equivalently, if Ex [ |∂ 2F (x) ∂xi∂xj |2 ] > 0, namely the partial derivative w.r.t xj turns out to be dependent on xi, then we say xi and xj are interacted. Otherwise, xi and xj have no interaction, if F (x) can be expressed as the sum of two functions f\i and f\j (Sorokina et al., 2008), namely,
F (x) = f\i (x1, . . . , xi−1, xi+1, . . . , xp) + f\j (x1, . . . , xj−1, xj+1, . . . , xp) ,
where f\i(j) is irrespective to xi(j). Similarly, higher-order (multi-way) interaction can be defined. In this paper, we mainly focus on pairwise interaction, while multi-way interaction is only briefly discussed due to space limitations.

heading: 3.1 INTERACTION STRENGTH MEASURE
text: Motivated by the above definition of interaction, it is natural to take advantage of the Hessian matrix H := ∇2xxF (x). The magnitude of its entry | ∂2F (x) ∂xi∂xj
| contains rich information about the local curvature at a data point x. Note that F (x) is a regression function here, not a loss function. This idea was originally considered in the economics community (Ai & Norton, 2003) and rediscovered for sensitivity analysis in Roustant et al. (2014). The goodness of F (x) as an approximator can essentially influence the interaction detection performance, see our Theorem M.2 in the supplement.
We define g(x, i, j) := |∂ 2F (x) ∂xi∂xj |2 to measure the local interaction strength of the i-th and j-th features at point x. We then use f(i, j) := Ex [g(x, i, j)] = Ex [ |∂ 2F (x) ∂xi∂xj |2 ]
as a measure of the global interaction strength. If f(i, j) ≈ 0, we say feature xi and xj have weak interaction; otherwise, if f(i, j) is significantly larger than zero, then xi and xj have strong interaction. In this paper, we focus on the global interaction.
3.2 INTERACTION STRENGTH EVALUATION
The above defined global interaction strength f(i, j) = Ex [ |∂ 2F (x) ∂xi∂xj |2 ]
is mostly unavailable due to the unknown input distribution x ∼ P (x). The Monte Carlo method can be used to approximate it by computing the sample mean over the training data(Roustant et al., 2014).
Analytical Evaluation: The Hessian matrix∇2xxF (x) for neural networks at a certain data point can be calculated analytically and efficiently, by using the automatic differentiation (Paszke et al., 2017). However, using this analytical solution is problematic for some learning models, such as the ReLU network, Random Forest (RF), etc., see supplement A. For example, the landscape of a ReLU network is piece-wise linear as shown in Figure 1, thus the exact Hessian is a zero matrix at almost every point. So, we turn to the following numerical evaluation for broader horizons.
Numerical Evaluation: Finite difference method is a common way to approximate the Hessian on a given data point (Campolongo & Braddock, 1999), i.e.,
∂2F (x)
∂xi∂xj ≈ 1 4hihj [F (x + eihi + ejhj)− F (x + eihi − ejhj)
− F (x + ejhj − eihi) + F (x− eihi − ejhj)], (1)
where ei is a one-hot vector with the i-th element being equal to one and the rest of elements being zeros. We also note that if the computation of gradient ∂F∂x is cheap (e.g., F is a neural network), then Hij can be approximated as 12h [ ∂F (x+ejh) ∂xi − ∂F (x−ejh)∂xi ] to reduce computation, which is similar to the feature interaction score defined in Greenside et al. (2018). The choice of perturbation size hi or hj (abbr. hi(j)) is critical. Generally, we do not want hi(j) to be too small (incurring round-off error) or too large (incurring truncation error) to get a good overall approximation of the derivative (Jerrell, 1997; Baydin et al., 2018). For our problem, we particularly do not want hi(j) to be too small so that the four evaluated points (shown in the numerator of Equation 1) lie on the same hyperplane (e.g., region A in Figure 1), which makes the quantity in Equation 1 always zero. The following theorem reveals that choosing a sufficiently small hi(j) is not necessary.
Theorem 3.1. For any x and y, function F shows no interaction between them, i.e., it can be decomposed as F (x, y) = a(x) + b(y) if and only if, for any h, k > 0, F (x + h, y + k) − F (x + h, y − k)− F (x− h, y + k) + F (x− h, y − k) = 0.
The magnitude of the numerator in Equation 1 tells us whether F is locally separable for variables xi and xj at point x, see our proof in supplement B.
The main issue of the finite difference method lies in the computational complexity for approximating the global interaction strength f(i, j), especially when the function evaluation itself is expensive. Using all the training samples for finding just a few strongest interaction pairs can be a total waste of computation resources. In general, for a dataset with N samples in p-dimensional feature spaces, a total number of 4Np(p− 1)/2 arithmetic evaluations of the surrogate regression function F are needed. If the Hessian happens to be sparse, that is, there are only a few interactions existing in the ground truth function, massive evaluations on those feature pairs with zero interaction strength should be avoided.
3.3 IDENTIFICATION OF THE k-STRONGEST PAIRWISE INTERACTIONS
In practice, we do not have to obtain the interaction strengths for all interaction pairs, because many of them are simply too weak to have an impact on the output, and the top k strongest pairwise interactions are well sufficient for data modeling and prediction. Finding the top k strongest pairwise interactions fits perfectly into the best k-arms identification problem in the context of multi-armed bandits.
Inspired by the recent work (Bagaria et al., 2018b), we can treat each entry of the Hessian Hij as an arm Ar (see Figure 2). Since Hessian is a symmetric matrix, we only need to consider n = p(p− 1)/2 arms in the set {Ar : r ∈ [n]}. For example, if we choose to pull the arm for the i-th row and j-th column of the Hessian, we will evaluate the local interaction strength g(xξ, i, j) on a randomly selected training data point xξ. The strength g(xξ, i, j) can be regarded as a random reward, where xξ is uniformly drawn from the training data {x1, . . . ,xN}. For ease of notation, we let f(r; ξ) := g(xξ, i, j) as the random reward of the arm r, and define µr := Eξ[f(r; ξ)] as the true mean value. Our goal is to find the best k arms with the highest mean reward.
This reformulation is reasonable since the following three essential assumptions for multi-armed bandits (Slivkins et al., 2019) are satisfied.
a. Only the reward will be observed after each pull; b. The reward for each arm is drawn independently and identically from its reward distribution; c. The reward for each round is bounded. Since we only consider functions F (x) defined on a
compact set, if F is a continuous function and twice differentiable, its first- and second-order derivatives are bounded accordingly.
3.4 DETECTING THE k-STRONGEST PAIRWISE INTERACTIONS WITH UCB ALGORITHM
In this section, we briefly introduce the UCB algorithm (Lai & Robbins, 1985) for finding the k-strongest interactions. Here, we design a sequence of estimators {f̂`(r)} for the mean reward of the arm Ar and construct confidence interval C(`). We let f̂`(r) denote the estimator after ` pulls of the arm Ar. Several assumptions need to be made before we give out the primary outcome.
Assumption 1. Finite m evaluations of one arm are sufficient to obtain an accurate reward. Assumption 2. Estimators f̂`(r) are σr-subgaussian.
Assumption 2 is satisfied because the reward is bounded (Hoeffding’s lemma, 1963). Here, σrs are pre-defined parameters. In the following analysis, we will choose σ = maxr{σr}, such that all f̂`(r) are σ-subgaussian.
Let f̂`(r) := 1` ∑` i=1 f(r; ξi), and we construct the 1− δ confidence intervals (see supplement C) of f̂`(r) as,
C(`) :=
{ √ 2σ2 log 2δ
` if ` ≤ m 0 if ` > m , (2)
µr ∈ [ f̂`(r)− C(`), f̂`(r) + C(`) ] , w.p. 1− δ, (3)
where µr denotes the true expectation. Increasing number of pulls of one arm, the uncertainty of the reward for that arm will decrease. When the number of pulls reaches m, we can remove the uncertainty of the arm and set C(`) = 0. Let `r(t) count the number of pulls of arm Ar till iteration t. To simplify the notations, we let µ̂r(t) = f̂`r(t)(r) and Cr(t) = C(`r(t)) be the sample means and confidence intervals at iteration t. Now, we formally introduce the complete procedure in Algorithm 1.
The Algorithm 1 is initialized by evaluating all arms O(log(n)) times, then at each subsequent iteration, it picks the arm with the highest UCB, i.e., maxj{µ̂j(t)+Cj(t)}, to evaluate (see Figure 2). If one arm has been evaluated for m times, we will its true mean value and set the uncertainty to zero. At the end of each iteration, if there exists one arm whose Lower Confidence Bound (LCB) is higher than all the other’s UCB, we will put it into the set of k best arms.
In the following, we give a theoretical analysis on the complexity of this algorithm. We have to further define ∆(k)i := max(0, µi∗k − µi), where i ∗ k is the index for the k-th best arm and ∆ (k) i indicates the gap of true values between the top-k arms and the rest. Theorem 3.2. With probability 1−Θ( 1n2 ), given maximal pulling times m for each arm, Algorithm 1
returns the k-strongest interaction pairs in O (∑n
i=1 log(n)
( σ2 log(nm)
(∆ (k) i )
2 ∧m
)) time, where (· ∧ ·)
is short for min(·, ·).
Remark 3.1. When k n, the complexity of Algorithm 1 for finding the top k interactions is O(n log(mn) log(n) + km log(n)) under some natural assumptions on the distribution of ∆i, being superior to the naive algorithm (every arm is pulled m times) with complexity O(nm). Remark 3.2. For k = 1, namely the strongest interaction pair identification, it requires at most∑n i=1(( 8σ2
(∆ (1) i )
2 log(n3m)) ∧m) pulls. We leave the complete proof in supplement C.

heading: 3.5 EXTENSION TO MULTI-WAY INTERACTION
text: Similarly, the third-order derivatives ∂ 3F
∂xi∂xj∂xk ,∀(i, j, k), and even the higher-order ones can be
approximated by the finite difference method as well, but the computational complexity will increase geometrically. For a remedy, one can detect higher-order interactions from the obtained pairwise ones. For instance, based on the k-strongest pairwise interaction strengths, we can further construct an undirected graph G, in which each node represents a feature. We develop Algorithm 2 to detect the multi-way interactions through checking all the cliques in the graph G, and more details are given in supplement G. Despite the computation complexity, our proposed algorithm is more advantageous to use because the associated learning model is agnostic and the form of interaction is adaptive as compared with some existing methods, such as Min et al. (2014), which relies on a logistic regression model and multiplicative interactions of a particular form.

heading: 4 MODEL COMPRESSION WITH DETECTED INTERACTIONS
text: The detected interactions can be used for designing structured and more transparent models. A similar idea was mentioned in Tsang et al. (2020b). But we propose a new architecture called Parametric ACE (ParaACE), inspired by several seminal works on Generalized Additive Model (GAM) (Hastie, 2017) and (nonparametric) Alternating Conditional Expectation (ACE) method (Breiman & Friedman, 1985). The original ACE method is based on the following model,
θ(Y ) = p∑ k=1 φk(Xk) + , (4)
and solves for the optimal transformation functions {θ, φ1, · · · , φp} nonparametrically in the Hilbert space. However, it only considers transformation of each single feature and is time consuming for large feature dimension p. We build the ParaACE model with the detected interactions as
h(x; {θsi}, {θri},β) := β0 + p∑ i=1 βisi (xi;θsi) + R∑ i=1 βp+iri (xIi ;θri) , (5)
where {Ii}Ri=1 is a set of indices for the detected interaction pairs, si and ri are the transformation functions for a single feature and interacting features with the corresponding parameters {θsi} and {θri}, β is a weight vector for the output of all the transformation functions.
The new architecture is shown in Figure 3 with two meta layers: the optimal feature transformation layer (Equation 5, boxed in red and green) and an optional fix-up layer (boxed in yellow). The optimal transformation layer gives the transformation functions of both the main effects and interacting features explicitly. The fix-up layer can be understood as the inverse function for the optimal transformation of the output, θ(Y ), in the original ACE (Equation 4). Another functionality of the fix-up layer is to alleviate the negative impacts of the wrongly detected interactions on the output. More explanations are given in supplement H.
We highlight the advantages of this model as follows: 1) interpretable, unlike many other models, ParaACE model
gives explicit transformations of features, making it interpretable. 2) lightweight, this model can significantly reduce the number of parameters from the original one, and unlike the pruning technique (Han et al., 2015), it does not need extra memory to store the indices of the preserved weights.
3) flexible, the structure of NN is modularized by blocks, which is easy to plug in/off for other applications. The fix-up layer can be redesigned as block-wise subnetworks to capture high-level interactions, end up in the shape of a hierarchical deep network.

heading: 5 EXPERIMENTS

heading: 5.1 EXPERIMENTAL SETUP
text: Datasets: We generated 10 synthetic datasets (p = 10) as was used in Tsang et al. (2018a) to test the accuracy of different interaction detection methods. We added a scaled Gaussian noise η ∼ 0.1×N(0, 1) to get the noisy data. We also selected 5 real datasets, namely the Elevators for controlling an F16 aircraft (Itorgo, 2019), Parkinsons for predicting the total UPDRS scores (Tsanas et al., 2009), Skillcraft for game player behavior analysis (Thompson et al., 2013), Cal housing for house price prediction (Pace & Barry, 1997), and Bike sharing for predicting hourly bike rental count (Fanaee-T & Gama, 2014). The datasets are preprocessed, and details are shown in supplement D.
Neural Network Configuration: For the baseline OverparaFC (Teacher), we used the architecture of p-5000-900-400-100-30-1 (indicating a deep neural network with p inputs, one output, and 5 hidden layers with 5000, 900, 400, 100, 30 neurons, respectively), which is a standard ReLU network with the last hidden layer being linear. For the ParaACE model (Student), the size of the network mainly depends on the number of candidate interactions we choose; more interactions imply more blocks. However, the configuration for each subnetwork is fixed, the structure of which is 1-50-8-1 for main effects and 2-50-8-1 for pairwise interactions. These subnetworks are ReLU networks with linear output. The fix-up layer is chosen as a single layer ResNet with the number of neurons equal to 15. All the above networks were initialized with Kaiming’s strategy (He et al., 2015).
Optimizer Setup: We chose adam (Kingma & Ba, 2014) as the optimizer, with the default setup in MXNet (Chen et al., 2015), and the batch size is set to be 500 for all the datasets.

heading: 5.2 PAIRWISE INTERACTION DETECTION
text: We first test on 10 different synthetic datasets (p = 10) (Tsang et al., 2018a) to show the superiority of our proposed method. Datasets with higher input dimensions (p = 50, 100) were compared as well, see supplement F. To compare the detection performance of different methods, we adopt the ROC-AUC metric (see supplement E). In Table 1, we compared our proposed method with various existing methods, including ANOVA (Wonnacott & Wonnacott, 1990), AG (Sorokina et al., 2008), IH (Janizek et al., 2020), NID (Tsang et al., 2018a), and PID (Liu et al., 2020). HierLasso (Bien et al., 2013) and RuleFit (Friedman et al., 2008) have been omitted due to their weak performances. For our method, we pick the same number of data samples n = 10000 for the training of OverparaFC and run 5 independent Monte Carlo tests. In Algorithm 1, we pull each arm 3 times for initialization. Each arm will be pulled at maximally m = 100 times, and we terminate when k = 20 strongest interactions stand out. We take µ̂r(t) upon termination as the final scores for our method. In this experiment, our proposed method needs around 1500 pulls of arms to pick out the top 20 interactions, however, naively pulling each arm 100 times needs 100C210 = 4500 pulls in total.
From the results, we can conclude the following: (1) PID outperforms all others; (2) our proposed method is as competitive as the PID method; (3) while the average performance of our method is only slightly worse than the PID method due to the poor result on dataset 7; this is due to the term 1 1+(x4x5x6x7x8)2
in F7(x) has little impact on the output, and such interactions are hard to detect. Apart from the above facts, both NID and PID work by analyzing the connectivity of NN, which strongly relies on specific sparse ReLU networks and empirical interaction strength computations, and thus are risky to generalize and may fail for some other datasets. For instance, for the function F11(x) = x 2 1x2x 2 3x4 + x 2 5x6x 2 7x8, the ROC-AUC of NID degrades to around 0.69, while our method maintains 1.0 robustly. Moreover, our method is model-agnostic and can be directly applied to any pre-trained learning models, such as tree models, kernel regression models, etc., which is impossible for other methods.

heading: 5.3 MODEL COMPRESSION WITH INTERACTION KNOWLEDGE
text: In this section, we construct a Student model with the learned interactions. With the interaction knowledge extracted from the Teacher model, we first show that the Student model trained under the same setup can achieve well-improved prediction accuracy (on average). We use normalized RMSE (NRMSE), which is defined by RMSEymax−ymin to measure the performances of all competing models such that the comparison across different datasets is fair and convenient.
For synthetic data, we use 800 training samples and 200 test samples. The noise η is injected as described earlier. We choose the top 20 pairwise interactions with our new detection method to build the ParaACE model. We compared our proposed ParaACE with OverparaFC, KD for regression (Chen et al., 2017), and two recent pruning methods known as Lottery Ticket Hypothesis (LTH) (Frankle & Carbin, 2019) and SynFlow (Tanaka et al., 2020). The results are shown in Table 2. As we expected, the KD method shows similar performance as the OverparaFC, because the Student model tries hard to approximate its Teacher. However, our ParaACE model can improve the performance of the OverparaFC by 26.0% on average, and the model compression ratio (CR := # parameters in the uncompressed model # parameters in the compressed model ) is around 300. More details can be found in supplement L. We also found that the ParaACE model is sample efficient. To be concrete, the ParaACE trained on a small training set can surpass the over-parameterized NN trained with substantially more data, see Figure 13 in supplement K.
For real-world datasets, we built the ParaACE model with the top 50 pairwise interactions for Elevators, Parkinsons, Skillcraft, Bike sharing, and the top 20 pairwise interactions for Cal housing (since the number of features is smaller). For most datasets, the detected interactions can help improve the predictive performance and significantly reduce the model size, see Table 3 for detailed results. The comparison with KD is omitted as it is supposed to be close to the Teacher.
Due to space limitations, we give more results for classification tasks in supplement I. Furthermore, we show the performance gain induced by interaction detection via ablation study in supplement N.

heading: 5.4 MODEL INTERPRETATION
text: For synthetic data, we take F10(x) = sinh (x1 + x2) + arccos (tanh(x3 + x5 + x7)) + cos(x4 + x5) + sec(x7x9) as an example. Three dominant pairwise interactions {x1, x2}, {x4, x5}, {x7, x9} have been detected successfully and visualized in Figure 4 (left). The ground truth transformations are shown in the upper part of the figure. Besides, the transformations for the single features can also be learned accurately. Taking F3(x) = exp |x1− x2|+ |x2x3| − (x23)|x4|+ log(x24 + x25 + x27 + x28) + x9 +
1 1+x210 for example, wherein feature x6 is not involved, x9 is linearly transformed, and x10 is transformed by an even function. All of them are correctly learned as shown in Figure 4 (right).
For real applications, we first take the Cal housing dataset for further analysis as less expert knowledge is required. From Figure 5, we can easily recognize that {longitude, latitude} and {totalRooms,totalBedrooms} show strong interactions. This makes sense, because {longitude, latitude} indicates the location, and {totalRooms,totalBedrooms} indicates the fraction of bedrooms. Also in the Parkinsons dataset, we find {Age, Sex} are interacted, which is consistent with the result that “age at onset (of Parkinson’s disease) was 2.1 years later in women than in men” in Haaxma et al. (2007). Interaction detection may also have a significant impact on drug combination discovery (Julkunen et al., 2020). We applied our interaction detection method on an RF model trained with the Drug combination dataset. Among the top 34 detected drug-drug interactions, 15 pairwise interactions have been verified in the DrugBank database (Wishart et al., 2018) (see supplement J). These results confirmed the effectiveness of our interaction detection method and strongly call for careful investigation of the unexplored drug pairs, potentially of great value.

heading: 6 CONCLUSION
text: This paper proposed a fast, generic, and model-agnostic interaction detection method. A thorough analysis of the computational complexity of our adaptive method proved its superior efficacy. The detected interaction knowledge can be further exploited to build a novel Student model, called ParaACE, which is compact, more interpretable, and mostly more accurate than its Teacher and other competing methods. We also showed experimentally that our ParaACE model is sample efficient owing to its simple architecture. Its model-agnostic property enables tremendous applications in different domains of finance, medicine, biology, wireless communication, where a significant number of well-performed black-box models have been built but not interpreted yet.

heading: ACKNOWLEDGMENTS
text: This work was supported by Guangdong Provincial Key Laboratory of Big Data Computing and by NSFC under Grant 61731018. The corresponding author is Feng Yin.

heading: A PROBLEMS WITH ANALYTICAL EVALUATION
text: In this section, we first illustrate the conceptual difference between local and global interaction and then state the problems of analytical evaluation.
Illustration of Local vs. Global Interaction. A quick example showing the difference of local/global interaction is the MATLAB symbol1 (Figure 6(a)). There is no local interaction inside the flat region, but two input variables globally interact.
Analytical evaluation works for neural networks with differentiable activation functions, e.g., sigmoid, tanh, etc. However, it is problematic when we use neural networks with piece-wise linear activation functions (PLNN), such as ReLU, Leaky ReLU, and so on. The Hessian will be a zero matrix at every second-order differentiable point, and it will no longer provide information about local interaction. Imagine the function landscape is joint with many facets, and the interaction information is hidden in the boundary of those facets (the set of non-differentiable points).
A simple example is checkerboard-like function, e.g., XOR function F (x1, x2) = 1{x1x2 > 0}, where x1, x2 are uniformly drawn from [−1, 1]. In this case, x1 and x2 are clearly interacted, however, the second-order derivatives ∂F∂x1∂x2 will be zero for almost all sampled points. The interaction information is hidden at the line x1 = 0 and x2 = 0.
B PROOF OF THEOREM 3.1 AND DISCUSSION ON THE CHOICE OF h
Theorem 3.1. For any x and y, function F shows no interaction between them, i.e., it can be decomposed as F (x, y) = a(x) + b(y) if and only if, ∀ h, k > 0, F (x+ h, y + k)− F (x+ h, y − k)− F (x− h, y + k) + F (x− h, y − k) = 0.
Proof. (⇒) is trivial. (⇐) We will directly have, ∀h, k > 0 and ∀x, y,
F (x+ h, y + k)− F (x+ h, y − k) k = F (x− h, y + k)− F (x− h, y − k) k ,
F (x+ h, y + k)− F (x− h, y + k) h = F (x+ h, y − k)− F (x− h, y − k) h .
1With permission of MathWorks.
Let k → 0, we will have ∂F∂y (x+ h, y) = ∂F ∂y (x− h, y) for any positive h; let h→ 0, we will have ∂F ∂x (x, y+ k) = ∂F ∂x (x, y− k) for any positive k. Thus, ∂F ∂y is irrespective of x and ∂F ∂x is irrespective of y.
Remark B.1. Theorem 3.1 can be extended for higher-order interaction with the same machinery.
A guide for the choice of hi(j). We empirically test the sensitivity of hi(j), see Figure 7. For standardized features, we suggest using hi(j) ∈ [0.6, 0.9] as the rule-of-thumb. Further discussions on the choice of hi(j). People may argue that for the interaction detection task, a small number of evaluations for each feature pair are already enough to get pretty good results, and the benefits of using the UCB algorithm is trivial. This is not always true. We address this from two perspectives in the following. First, in many scenarios, the function can only be evaluated sequentially, and the evaluation may be very expensive. The benefits of the UCB algorithm are revealed then. Second, we empirically study the performance of evaluating each arm evenly under different configurations (perturbation size hi(j), and the number of evaluations for each arm), see Figure 8. The ROC-AUC scores were obtained from the OverparaFC model trained as described in Section 5.2. We expect that given the fixed perturbation size h, the ROC-AUC scores will rise as the number of evaluations increases. This is consistent with the observation that the bricks are brighter from top to bottom. Also, we observed that with a properly chosen h ∈ [0.6, 0.9], the interactions can be correctly detected with a relatively small number of evaluations (the bricks are brighter in the middle, see dataset 3, 4, and 5). The corresponding reasons can be that (a) if h is too small, four evaluated points lie in the same flat region with high probability, which is harmful to interaction detection; (b) if h is too large, some evaluated points may be out of the distribution of the training data, which makes the detected interaction strength less representative. Even though the suggestion of using hi(j) ∈ [0.6, 0.9] is purely empirical, the plots here help us better understand the essence of statistical interactions and their detection process.

heading: C PROOF OF THEOREM 3.2
text: Recall that we construct the 1− δ confidence intervals of f̂`(r) as,
C(`) =
{ √ 2σ2 log 2δ
` if ` ≤ m 0 if ` > m ,
µr ∈ [ f̂`(r)− C(`), f̂`(r) + C(`) ] , w.p. 1− δ.
This is derived from concentration inequality (Inequality 2.9 in Wainwright (2019)). We set δ = 2n3m , which is a small number and easy to analyze. We have, Lemma C.1 (high probability event/clean event). With probability (w.p.) 1− 2n2 , all true values µr lie in the their confidence intervals during the run of the algorithm.
Proof. It’s equivalent to proof:
P ({ ∀ t,∀ r ∈ [n], |µr − µ̂r(t)| ≤ Cr(t)}) ≥ 1− 2
n2
The opposite of this event is
∃ t and r, s.t., |µr − µ̂r(t)| > Cr(t). Since we have,
P ({|µr − µ̂r(t)| > Cr(t)}) ≤ δ = 2
n3m ,
and we evaluate at most nm times (n arms, each arm is pulled m times),
P ({∃ t and r, s.t., |µr − µ̂r(t)| > Cr(t)}) ≤ δnm = 2
n2 (union bound)
P ({ ∀ t,∀ r ∈ [n], |µr − µ̂r(t)| ≤ Cr(t)}) ≥ 1− 2
n2 (clean event)
Using Lemma C.1, we restate Theorem 3.2 as follows.
Theorem C.2 (Restating Theorem 3.2). With probability 1 − Θ( 1n2 ), Algorithm 1 returns the kstrongest interactions with at most
M ≤ n∑ i=1 (( 24σ2 ∆2i log(nm) ) ∧m )
local interaction strength evaluations. This takes O (∑n i=1 log(n) (( σ2 log(nm) ∆2i ) ∧m )) time.
Proof. We first analyze the running of algorithm till the first strongest interaction comes out. If we choose to update arm i 6= i∗1 at time t, then we have
µ̂i(t) + Ci(t) ≥ µ̂i∗1 (t) + Ci∗1 (t). (6)
For equation 6 to occur, at least one of the following three events must occur: E1 = { µ̂i∗1 (t) ≤ µi∗1 − Ci∗1 (t) } ,
E2 = {µ̂i(t) ≥ µi + Ci(t)} , E3 = { ∆ (1) i = µi∗1 − µi ≤ 2Ci(t) } .
To see this, note that if none of E1, E2, E3 occurs, we have
µ̂i(t) + Ci(t) (¬E2) < µi + 2Ci(t) (¬E3) < µi∗1 (¬E1) < µ̂i∗1 + Ci∗1 (t).
From Lemma C.1, E1 and E2 do not occur during the run of the algorithm with probability 1− 2n2 , because
w.p. (
1− 2 n2
) : |µi − µ̂i(t)| ≤ Ci(t),∀ i ∈ [n], ∀ t. (7)
It also implies w.p. 1 − Θ( 1n2 ), the algorithm does not stop pulling arm i until event E3 stops occurring.
Let ζ(w)i index the iteration in which Algorithm 1 evaluates arm Ai for the last time before declaring it to be the w-th strongest interaction. Let’s first consider how many pulls are needed for each arm to
pick out the first strongest interaction. If Ci(ζ (1) i ) ≤
∆ (1) i
2 for arm Ai, then we can stop evaluating arm Ai, that is,
∆ (1) i 2 ≥ √ 2σ2 log n3m
Ti(ζ (1) i )
or Ci(ζ (1) i ) = 0,
=⇒ Ti(ζ(1)i ) ≥ 8σ2
(∆ (1) i )
2 log(n3m) or Ti(ζ (1) i ) ≥ m.
Note that Ti(t) denotes for the number of pulls of arm Ai at iteration t. We will evaluate arm Ai at most 8σ 2
(∆ (1) i )
2 log(n3m) ∧m times. To pick out the first strongest interaction, we need at most M
evaluations satisfying
M ≤ n∑ i=1
(( 8σ2
(∆ (1) i )
2 log(n3m)
) ∧m ) .
This result can be extend to find all top k strongest interactions concretely,
M ≤ n∑ i=1
(( 8σ2
(∆ (k) i )
2 log(n3m)
) ∧m ) ,
= n∑ i=1
(( 24σ2
(∆ (k) i )
2 log(n 3
√ m) ) ∧m ) ,
≤ n∑ i=1
(( 24σ2
(∆ (k) i )
2 log(nm)
) ∧m ) .
The Algorithm 1 takes O (∑n i=1 log(n) (( σ2 log(nm) ∆2i ) ∧m )) time, where O(log(n)) is the time for maintaining a priority queue(to find the minimal LCB or maximal UCB) in each iteration.
Theorem C.3 (Remark 3.1). If we further assume that ∆i ∼ N (γ, 1), for some γ, and m = cn, for c ∈ [0, 1], then the expected pulls of arms (over randomness in ∆i) satisfies
E[M ] ≤ O(n log(nm) + km)
with probability 1−Θ( 1n2 ) over randomness in Algorithm 1. Thus, the expected running time is well bounded by O(n log(nm) log(n) + km log(n)).
Proof. This follows directly from the Appendix 2 of Bagaria et al. (2018a).

heading: D DATASETS PRE-PROCESSING
text: D.1 SYNTHETIC TEST SUITE
To make our experimental results convincing, we follow the test suite ever used in Tsang et al. (2018a) with the details given in Table 4. Among others, F1 is a widely used test function for interaction detection, which can be generated as described in Sorokina et al. (2008). For all the other functions, the input dimension is set to 10, and x1, . . . , x10 i.i.d.∼ U(−1, 1).
Table 4: Test suite of data-generating functions
F1(x) π x1x2 √ 2x3 − sin−1(x4) + log(x3 + x5)− x9 x10 √ x7 x8 − x2x7
F2(x) π x1x2 √ 2|x3| − sin−1(0.5x4) + log(|x3 + x5|+ 1) +
x9 1 + |x10|
√ x7
1 + |x8| − x2x7
F3(x) exp |x1 − x2|+ |x2x3| − (x23)|x4| + log(x24 + x25 + x27 + x28) + x9 + 1
1 + x210
F4(x) exp |x1 − x2|+ |x2x3| − (x23)|x4| + (x1x4)2 + log(x24 + x25 + x27 + x28) + x9 + 1
1 + x210
F5(x) 1
1 + x21 + x 2 2 + x 2 3
+ √ exp(x4 + x5) + |x6 + x7|+ x8x9x10
F6(x) exp (|x1x2|+ 1)− exp(|x3 + x4|+ 1) + cos(x5 + x6 − x8) + √ x28 + x 2 9 + x 2 10
F7(x) (arctan(x1) + arctan(x2)) 2 + max(x3x4 + x6, 0)−
1
1 + (x4x5x6x7x8)2 +
( |x7|
1 + |x9|
)5 + 10∑ i=1 xi
F8(x) x1x2 + 2 x3+x5+x6 + 2x3+x4+x5+x7 + sin(x7 sin(x8 + x9)) + arccos(0.9x10) F9(x) tanh(x1x2 + x3x4) √ |x5|+ exp(x5 + x6) + log ( (x6x7x8) 2 + 1 ) + x9x10 + 1 1 + |x10| F10(x) sinh (x1 + x2) + arccos (tanh(x3 + x5 + x7)) + cos(x4 + x5) + sec(x7x9)
D.2 REAL DATASETS
All the real datasets are publicly available. We preprocessed the datasets as follows. For Parkinsons data, we remove the column motor UPDRS and predict the total UPDRS. For SkillCraft data, the
target is set to be LeagueIndex. For the Bike sharing data, the feature weather is converted into a onehot representation. Before feeding the data into neural networks, we performed data normalization first for all datasets.
The Drug combination dataset is processed as follows. There are 110 features extracted, among which the first 50 features are the concentration of 50 unique FDA-approved drugs, and the last 60 features are the one-hot encodings for the cell lines.

heading: E ROC-AUC FOR PAIRWISE INTERACTION DETECTION
text: We calculate the ROC-AUC scores for the synthetic datasets, where the ground truth interaction pairs are known. To obtain the ROC-AUC value, two vectors are needed, namely the pairwise interaction score ∈ R45+ and ground truth ∈ {0, 1}45, both of them are of dimension p(p−1) 2 = 45. By setting different thresholds for the pairwise interaction score ranking, different classifiers with False Positive (FP) rate and True Positive (TP) rate can be obtained. The ROC-AUC value is approximated from those (FP, TP) pairs by the trapezoidal rule.

heading: F PAIRWISE INTERACTION DETECTION ON HIGH- DIMENSIONAL DATASETS
text: We created two new high-dimensional datasets by simply combining the 10 synthetic datasets considered in the paper. The new datasets have either 50 (using F1-F5) or 100 features (using F1-F10). The labels of the two datasets are the sum of the original labels in each dataset of input dimension p = 10. The training sample size is increased to 500,000 to mimic big data. We compared 5 different model architectures and our method consistently outperforms NID, see the Table 5 below for the results.
Here we also report the computational complexity of our UCB-based interaction detection method. Due to the randomness of the UCB algorithm, the number of pulls fluctuates. Note that we select top k = 200 interaction pairs for 100-dimensional case, and the total number of pulls is around 40000, while the naive approach needs 100× 100× 99/2 = 495000 pulls. We select top k = 100 interaction pairs for 50-dimensional case, and the total number of pulls is around 21000, while the naive approach needs 100× 50× 49/2 = 122500 pulls. Discussion on High Dimensional Data: For high dimensional data, the Hessian matrix will be too huge to handle. One possible solution is to focus on the important features. We may first take advantage of DeepLIFT(Shrikumar et al., 2017) or SHAP (Lundberg & Lee, 2017) to screen out the most important features, and then apply our interaction detection method. If the number of important features is still too large to handle, we may use ANOVA and F-test to select an affordable number of interaction candidates to pull log(n) times, and explore the remaining pairs with a small probability, say 5%, like the common strategy used in reinforcement learning..

heading: G EXTENSION TO HIGHER-ORDER INTERACTIONS
text: We define three-way interaction as follows, and higher-order interactions can be defined similarly.
Definition G.1 (three-way interaction). A function F : Rp → R is said to exhibit an interaction among three of its variables xi, xj and xk if,
Ex
[ ∂3F (x)
∂xi∂xj∂xk
]2 > 0.
The naive method to detect higher-order interactions is quite similar to FANOVA graph (Muehlenstaedt et al., 2012). We give three examples of detecting higher-order interactions from the synthetic data. We detect the pairwise interaction strength first, then build a weighted graph (the edge weight corresponds to the interaction strength). By setting a proper number of clusters, higher-order interactions can be discovered correctly, see Figure 9.
(a) y = x0x1 + x2x3x4 + x5x6x7 + x8x9
(b) y = e|x0−x1| ∗ x9 + |x2 ∗ x3| − x2|x5|4 + (x6x7)2 + x8 (c) y = x0x1 + |x1 + x2x3|+ x4x5 − (x25)x6 − ex7+x8x9
Algorithm 2 Hierarchical Higher-Order Interaction Detection
Require: The target number of i-way interactions {k(i)}pi=2. 1: Detect top-k(2) pairwise interaction via UCB algorithm. 2: Construct an undirected graph G based on the detected pairwise interactions. 3: Enumerate all the cliques C in the graph G. 4: for i = 3, 4, · · · , p do 5: if C(i) = ∅ then 6: break 7: end if 8: Find C(i)refine = {C ∈ C(i) | every i− 1 complete subgraph of C admits a detected interaction} 9: Detect the top k(i)-strongest i-way interactions in C(i)refine via the UCB algorithm. 10: end for 11: return All the detected interactions with their strengths.
One drawback of the above naive clustering method is that it doesn’t allow for the overlap of nodes, i.e. one variable may appear in several interactions (F8 in the test suite). To address the above problem, we propose the Algorithm 2, which detects the higher-order interactions hierarchically. We first construct an undirected graph G from the detected pairwise interactions. To shrink our search space for higher-order interactions, we restrict ourselves to the set of all cliques C of the graph G, where a clique C ∈ C in a graph is a subset of vertices that are all joined by edges. We use C(i) ⊆ C to denote the collection of the i-cliques (the cliques with i vertices). For example, to detect 3-way interactions, one only needs to search in the C(3), which contains the cliques with three vertices. Furthermore, we can shrink the search space to C(i)refine for i-way interactions, based on the detected (i− 1)-way interactions (see line 8 in Algorithm 2). For example, while detecting 4-way interactions, we put 4-clique {2, 3, 4, 6} into consideration only if all the 3-way interactions {3, 4, 6}, {2, 4, 6}, {2, 3, 6}, {3, 4, 6} exist. We then use the UCB algorithm to verify if the cliques
in C(i)refine admit the interaction relationship. In this way, we detect the i-way interactions from the detected (i− 1)-way interactions information. The 3-way interaction strength can be calculated from Equation 8. However, with the order of interaction increasing, the number of function evaluations for one gradient computation increases geometrically. We do not recommend detecting the interaction whose order is higher than four, which is uneconomical to compute and hard to interpret. We found our proposed Algorithm 2 can successfully detect the higher-order interactions for F8 in the test suite 2.
∂2F (x)
∂xi∂xj∂xk ≈ 1 8h3 [+F (x + h(ei + ej + ek)) + F (x + h(−ei − ej + ek))
+ F (x + h(ei − ej − ek)) + F (x + h(−ei + ej − ek)) − F (x + h(−ei + ej + ek))− F (x + h(ei − ej + ek)) − F (x + h(ei + ej − ek))− F (x + h(−ei − ej − ek))]
(8)

heading: H MORE ON THE PARAMETRIC ACE
text: In section 4, we introduced a parametric ACE model based on the following generalized linear additive model with main effects and pairwise feature interactions, namely,
θ(Y ) = p∑ i=1 βisi (xui ;θsi) + R∑ i=1 βp+iri (xIi ;θri) + . (9)
The consideration of interactions improves the model performance significantly, compare to Generalized Additive Neural Networks (GANN)(Potts, 1999; Agarwal et al., 2020)3, see Table 6.
The idea behind the original nonparametric ACE algorithm (Breiman & Friedman, 1985) is to alternate between an inner process for finding the optimal transformation functions of the inputs and an outer process for finding the optimal transformation function of the output.
Here, we made several modifications. First, all the transformation functions are represented by small subnetworks. Second, we reformulate the above regression equation approximately as
Y ≈ θ−1 (
p∑ i=1 βisi (xui ;θsi) + R∑ i=1 βp+iri (xIi ;θri)
) . (10)
Therefore, the fix-up layer shown in Figure 3 can be understood as the inverse optimal transformation of the output θ−1(·). To be more general, θ−1(·) takes all input transformations as individuals instead of their high-level summary (weighted sum). Our parametric ACE then alternately tunes the fix-up layer conditioned on the current optimal transformation layer (outer iteration) and tunes the optimal transformation layer conditioned on the current fix-up layer until some convergence condition is met.
Another important function of the fix-up layer is to alleviate the negative impacts of wrongly detected pairwise interactions and/or undetected higher-order interactions altogether on the output. For this purpose, we could make this fix-up layer a network of small subnetworks, so that the whole architecture is a network of small subnetworks. Each subnetwork can be regarded as a meta-neuron that is expected to be more competent than any SOTA activation function. With such a nice structure,
2A demo can be found at https://github.com/zhangtj1996/ParaACE. 3GANN is a parametric version of GAM, which only considers the transformations of univariates (single
features).
we have divided the hyper-parameters naturally into blocks and the whole network can hopefully be made adaptive to new tasks more rapidly by tuning just a few blocks.
It is possible for ParaACE to learn “contradictory” models, which present the same function (Lengerich et al., 2020). For example, the main effect might be absorbed into the interaction effects. Thus the direct interpretation through subnetworks may not sound that convincing. Lengerich et al. (2020) proposed an algorithm to purify the interaction effect by calculating the functional ANOVA (Hooker, 2007), based on a piecewise-constant function F̂ . This algorithm is applicable to any function F by first constructing a piecewise-constant approximation F̂ . Particularly, for our ParaACE model, the overall computation will be much easier since we explicitly have the univariate and bi-variate feature transformation. The number of bins required to approximate ParaACE with F̂ reduces significantly, and so does the complexity of the follow-up purifying algorithm.

heading: I EXTENSION TO CLASSIFICATION TASK
text: We generated the datasets for binary classification from the synthetic regression datasets (1000 samples with injected noise). We choose the median of the target as a threshold to separate each synthetic regression the dataset into two classes. The comparison between the baseline OverparaFC and our proposed ParaACE in terms of classification accuracy is shown in Table 7. For real-world datasets, we picked Higgs Boson (Baldi et al., 2014), Spambase (Dua & Graff, 2017), and Diabetes (Turney, 1994).Both the classification accuracy and the compression ratio are reported in Table 8.

heading: J DETECTED PAIRWISE INTERACTIONS FOR SYNTHETIC AND REAL DATA
text: Synthetic data: Figure 10 shows the pairwise interaction strength produced by our proposed method.
Real data: We provide the detected interactions for five real datasets in Figure 11. For the Cal housing dataset, we terminate the UCB algorithm when k = 20 interactions stood out; for the other four datasets, we terminate at k = 50. The green, orange, and blue points represent the UCB, estimated mean, and LCB respectively. For the drug combination dataset, the results are shown in Figure 12.

heading: K EMPIRICAL RESULT FOR SAMPLE EFFICIENCY
text: Figure 13 shows our proposed framework is sample efficient. The data are from the test suite (Table 4), and noise η is injected. Increasing the training sample size, both over-parameterized neural nets and our proposed ParaACE perform better. But ParaACE is less data-hungry since fewer training samples
are needed for ParaACE to achieve a similar performance of the over-parameterized neural nets. The experiments show that ParaACE still performs well, even if with small training samples.

heading: L DETAILS ON THE COMPARISON WITH KD, LTH, AND SYNFLOW
text: L.1 COMPARISON WITH KD
KD is widely used for multi-class classification problems, which extracts knowledge from the ”soft label” by controlling the temperature, but fewer KD methods are there for regression problems. Currently, people mainly use the Teacher bounded regression loss (Chen et al., 2017) and the hint loss to deal with the objective detection problem. Passing unlabeled data to the Teacher model to produce pseudo labels can also help. In these ways, the Student is expected to have a similar performance to the Teacher.
We let an over-parameterized FC (10-5000-900-400-100-30-1) be the Teacher, and a ReLU network with architecture (10-70-70-70-70-30-1) be the Student. The Student is trained with 800 original data and 4000 pseudo data (labeled by the Teacher). The loss function we adopted is L2 + Hint loss, where the Hint loss is used on the layer with 30 neurons.
It shows that KD trained with pseudo data and carefully designed loss function achieves similar performance as the Teacher, while ParaACE trained with original data and simple L2 loss achieves significantly better performance.
L.2 COMPARISON WITH LTH
We used the pruning technique proposed in the lottery tickets hypothesis (LTH) paper (Frankle & Carbin, 2019). We pruned the fully connected neural network after every 20 epochs. Every time 20 % weights were pruned in each layer except for the last layer. We trained 500 epochs, and set the batch size to be 400. We report the best test performance in each round of pruning in Figure 14.
Figure 14 shows that, in most cases pruning the network properly can improve the model performance. But if the network was over-pruned, the performance may drop.
L.3 COMPARISON WITH SYNFLOW
We implemented the single shot SynFlow for regression tasks based on the official code repository https://github.com/ganguli-lab/Synaptic-Flow (Tanaka et al., 2020). The sparsity is set to 10−2.447 for synthetic datasets, thus the neural network is compressed 280 times. For real-world datasets, we
compressed the model by 120 times. And we retrained the neural network after the single-shot prune (post-training). Note that we did not prune the biases here.
We set the hyper parameters for the post-training process as follows. The batch size is set to be 500, and the number of epochs is 100. We use Adam as the optimizer with lr = 0.001, betas = (0.9, 0.99).

heading: M ACCURATE APPROXIMATION FUNCTION BENEFITS INTERACTION DETECTION
text: Suppose the underlying ground truth function is f(x), and our approximation function is gθ(x), we omit the θ in the following. We show that the estimation error for Hessian∣∣∣Ex [ ∂2g(x)∂xi∂xj ]− Ex [ ∂2f(x)∂xi∂xj ]∣∣∣ is bounded by o( ), when |f(x)− g(x)| < for all x. This implies more accurate pre-trained model leads to better interaction detection accuracy. According to the following theorem, deep neural networks are good surrogate functions due to the universal function approximation capability (Hornik, 1991).
Theorem M.1. Assume xi are uniformly drawn from [−1, 1] independently, if there exist an > 0, such that |f(x)− g(x)| ≤ for all x ∈ Rp, then
Ex
[ ∂2g(x)
∂xi∂xj
] − Ex [ ∂2f(x)
∂xi∂xj
] ≤ .
Proof.
Ex
[ ∂2g(x)
∂xi∂xj ] = ∫ ∂2g(x)
∂xi∂xj p(x)dx
=
∫ ∂2g(x)
∂xi∂xj p(xi, xj)p(x\ij)dx
= ∫ p(x\ij) ∫ 1 −1 ∫ 1 −1 ∂2g(x) ∂xi∂xj p(xi, xj)dxidxjdx\ij
[where p(xi, xj) = 1/4]
= 1
4
∫ p(x\ij) ∫ 1 −1 ∫ 1 −1 ∂2g(x) ∂xi∂xj dxidxjdx\ij ,
where ∫ 1 −1 ∫ 1 −1 ∂2g(x) ∂xi∂xj dxidxj
= ∫ 1 −1 ∂g(x\ij , xi = 1, xj)− ∂g(x\ij , xi = −1, xj) ∂xj dxj
= g(x\ij , xi = 1, xj = 1)− g(x\ij , xi = 1, xj = −1) − (g(x\ij , xi = −1, xj = 1)− g(x\ij , xi = −1, xj = −1)).
We can derive Ex [ ∂2f(x) ∂xi∂xj ] in a similar fashion. Thus,
Ex
[ ∂2g(x)
∂xi∂xj
] − Ex [ ∂2f(x)
∂xi∂xj ] = 1
4 Ex\ij [g(x\ij , xi = 1, xj = 1)− f(x\ij , xi = 1, xj = 1) − g(x\ij , xi = −1, xj = 1) + f(x\ij , xi = −1, xj = 1) − g(x\ij , xi = 1, xj = −1) + f(x\ij , xi = 1, xj = −1) + g(x\ij , xi = −1, xj = −1)− f(x\ij , xi = −1, xj = −1)]
≤1 4 Ex\ij [4 ]
= .
Theorem M.2. Assume xi are uniformly drawn from [−1, 1] independently, and further assume∣∣∣ ∂2f(x)∂xi∂xj ∣∣∣ , ∣∣∣ ∂2g(x)∂xi∂xj ∣∣∣ ≤ M,∀(xi, xj). If there exist an > 0, such that |f(x) − g(x)| ≤ for all x ∈ Rp, then
Ex [∣∣∣∣ ∂2g(x)∂xi∂xj ∣∣∣∣2 ] − Ex [∣∣∣∣∂2f(x)∂xi∂xj ∣∣∣∣2 ] ≤ 2M .
Proof.
Ex [∣∣∣∣ ∂2g(x)∂xi∂xj ∣∣∣∣2 ] − Ex [∣∣∣∣∂2f(x)∂xi∂xj ∣∣∣∣2 ]
=Ex [∣∣∣∣ ∂2g(x)∂xi∂xj ∣∣∣∣2 − ∣∣∣∣∂2f(x)∂xi∂xj ∣∣∣∣2 ]
=Ex
[( ∂2g(x)
∂xi∂xj − ∂
2f(x)
∂xi∂xj
)( ∂2g(x)
∂xi∂xj + ∂2f(x) ∂xi∂xj )] ≤ ∣∣∣∣Ex [( ∂2g(x)∂xi∂xj − ∂ 2f(x) ∂xi∂xj )]∣∣∣∣ · ∣∣∣∣Ex [( ∂2g(x)∂xi∂xj + ∂ 2f(x) ∂xi∂xj
)]∣∣∣∣ ≤ · 2M (Using Theorem M.1)
Note that M can be interpreted as the magnitude of the strongest local interaction. For the functions with strong interactions, higher approximation quality is desired.

heading: N ABLATION STUDY FOR THE EFFECTIVENESS OF INTERACTION DETECTION
text: To show how much performance gain we are able to obtain from the interaction detection procedure. We conduct an experiment that input random pairwise interactions to ParaACE to compare with the one with detected interactions on synthetic datasets, see Table 9.

Review summaries:
The authors propose two key ideas. The first is the idea of using the UCB algorithm to identify strong feature interactions in a computationally efficient way; each set of interacting features is an "arm" that could be "pulled", and pulling the arm corresponds to evaluating the strength of the interaction by computing the corresponding entry in the Hessian on a random training example. The finite difference method is used for computing Hessians. The second key idea is that of using the identified pairwise interactions to build a lightweight GAM-like model that they call ParaACE. In experiments, the authors show the UCB approach is effective at identifying feature interactions compared to alternative methods, and they demonstrate that ParaACE offers strong gains in model compression compared to other competing approaches, and often improves performance relative to its overparameterized teacher model.

Identification of feature interactions from black-box models, on a global as well as on a local scale, can lead to a better understanding of the operating modes of the black-box models. In this paper, the authors presented a principled approach to identify global interactions by casting the problem as a multi-arm bandit problem and proposes a solution using UCB algorithm. The authors claim that their method interaction discovery method is free of ad-hoc assumptions with good detection accuracy and stability. Furthermore, the authors showcase the importance of the learned interactions by proposing a new deep learning model based on these interactions and showcase the improvements in model size (thereby competing against pruning methods) as well as in accuracy (thereby competing against generalization methods). 

The paper handles the important area of model that can be interpreted using feature interactions - aligned to theme of explainable deep learning. The authors chose to use the problem of multi-arm bandit, solving it by UCB algorithm with good speed and accuracy. A lightweight and interpretable deep learning model (called ParaACE), built using alternating conditional expectation (ACE) method is the crux of the work. It is shown that the proposed method improves accuracy by 26% and reduces the model size by 100+ times as compared to its Teacher model over various datasets. The paper has extensive supplementary material, as well as shared code.