# Higher-order Neural Additive Models: An Interpretable Machine Learning Model
with Feature Interactions

Minkyu Kim,1 Hyun-Soo Choi,2 Jinho Kim2

###### Abstract

Black-box models, such as deep neural networks, exhibit superior predictive
performances, but understanding their behavior is notoriously difficult. Many
explainable artificial intelligence methods have been proposed to reveal the
decision-making processes of black box models. However, their applications in
high-stakes domains remain limited. Recently proposed neural additive models
(NAM) have achieved state-of-the-art interpretable machine learning. NAM can
provide straightforward interpretations with slight performance sacrifices
compared with multi-layer perceptron. However, NAM can only model
1${}^{\text{st}}$-order feature interactions; thus, it cannot capture the co-
relationships between input features. To overcome this problem, we propose a
novel interpretable machine learning method called higher-order neural
additive models (HONAM) and a feature interaction method for high
interpretability. HONAM can model arbitrary orders of feature interactions.
Therefore, it can provide the high predictive performance and interpretability
that high-stakes domains need. In addition, we propose a novel hidden unit to
effectively learn sharp-shape functions. We conducted experiments using
various real-world datasets to examine the effectiveness of HONAM.
Furthermore, we demonstrate that HONAM can achieve fair AI with a slight
performance sacrifice. The source code for HONAM is publicly
available111https://github.com/gim4855744/HONAM.

## Introduction

Deep neural networks have achieved great success in various fields, such as
computer vision and natural language processing. However, understanding and
interpreting the decision processes of black-box models, such as neural
networks, are notoriously difficult (Agarwal et al. 2020). Therefore, their
application is limited to high-stakes domains, such as healthcare and social
safety. Therefore, explainable artificial intelligence (AI) methods to explain
the behaviors of black box models have been widely studied in recent years.
Local interpretable model-agnostic explanations (LIME) (Ribeiro, Singh, and
Guestrin 2016) is a representative explainable AI method. It attempts to
locally approximate a black-box model using an interpretable model to explain
individual predictions. LIME can explain a local area of the black-box model
if the interpretable model has approximated the black-box model to be
explained well. However, LIME often fails to explain the global behavior of
the black box model. Shapley additive explanation (SHAP) (Lundberg and Lee
2017) is a game-theoretic approach for model-agnostic explanations. SHAP
generates predictions of a black-box model based on possible combinations of
input features and obtains the effect of a feature by exploiting changes in
predictions when the feature is dropped. However, the explanation of SHAP may
not be faithful to the actual predictions of the black-box model. As mentioned
above, explainable AI methods are still limited in their applications to high-
stakes domains because they often fail to explain and may not be faithful to
the actual behaviors of black-box models; thus, we need interpretable AI
methods (Rudin 2019).

The recently proposed neural additive models (NAM) (Agarwal et al. 2020) are
state-of-the-art methods of interpretable machine learning. Since NAM learns a
linear combination of neural networks corresponding to each input feature, it
can model nonlinear relationships between each input feature and the output.
In addition, NAM can provide straightforward interpretations and has a
predictive performance comparable to that of multi-layer perceptron (MLP) and
tree-based methods, such as XGBoost (Chen and Guestrin 2016). However, NAM can
only model 1${}^{\text{st}}$-order feature interactions. Therefore, NAM cannot
capture the co-relationships between input features even though there are many
such relationships in real-world data. More recently, (Chang, Caruana, and
Goldenberg 2022) have proposed NodeGAM and NodeGA${}^{\text{2}}$M which are
neural decision tree-based generalized additive models (GAMs). However,
NodeGAM and NodeGA${}^{\text{2}}$M also merely can learn
1${}^{\text{st}}$-order and 2${}^{\text{nd}}$-order feature interactions,
respectively. If there is an interpretable machine learning method that can
model high-order feature interactions, it will not only provide high-quality
interpretations, but also show better predictive performance.

Methods for explicitly modeling feature interactions using machine learning
have been studied extensively (Rendle 2010; Cheng et al. 2016; Wang et al.
2017; Guo et al. 2017; Lian et al. 2018). Feature interaction methods aim to
model high-order interactions through feature crossing, which iteratively
multiplies input features, and has shown impressive performance in predictive
tasks, such as recommendation (Lian et al. 2018) and regression (Kim, Lee, and
Kim 2020). However, there are several problems with the existing feature
interaction methods. First, explicit feature interaction methods are linear.
Therefore, they cannot learn nonlinear functions, and their expressive powers
are limited. Second, most existing feature interaction modeling studies have
mainly focused on improving predictive performance. Although linear models are
interpretable, their interpretability tends to be neglected in previous
studies.

Since sharp jumps are frequently found in real-world data, making machine-
learning models that can learn sharp jumps is a significant challenge in
interpretable models (Agarwal et al. 2020; Chang et al. 2021). However, neural
networks struggle to capture such sharp jumps (Agarwal et al. 2020) because of
the bias towards smooth functions (Arpit et al. 2017; Rahaman et al. 2019).
(Agarwal et al. 2020) have proposed the exp centered (ExU) unit to address
this problem. However, the ExU unit merely works well under certain
circumstances.

To overcome these problems, we propose a novel interpretable machine learning
method called higher-order neural additive models (HONAM). Since HONAM learns
a linear combination of neural networks corresponding to each input feature
like NAM, it can not only effectively model non-linear relationships between
each input feature and the output, but is also still interpretable. The major
difference between HONAM and NAM is that HONAM can model high-order feature
interactions. It allows HONAM to capture the co-relationships between features
and provides high-quality interpretations and better predictive performance
than NAM. Furthermore, we propose a novel feature interaction modeling method
for high interpretability because we may struggle to interpret the predictions
of HONAM when using the existing feature interaction methods. We conduct
extensive experiments to examine its effectiveness and demonstrate the
interpretability of the HONAM. Experimental results demonstrate that HONAM has
a higher or comparable predictive performance than MLP and XGBoost, even
though it is an interpretable method and can provide high-quality
interpretations. In addition, we propose a novel hidden unit named ExpDive to
learn sharp shape functions, and we show that a trained HONAM can become a
fair AI by removing biases.

Figure 1: The Architecture of HONAM

## Higher-order Neural Additive Models

(a) 1${}^{\text{st}}$-order
features

(b) 2${}^{\text{nd}}$-order features

Figure 2: Feature contributions on the Insurance dataset learned by the
CrossNet

In this section, we introduce HONAM, a novel interpretable machine-learning
method that can model arbitrary orders of feature interactions. First, we
describe NAM (Agarwal et al. 2020) and transform it to model feature
interactions. We then present HONAM with element-wise feature crossing.
Finally, we propose a novel feature interaction method for high
interpretability and apply it to HONAM.

### Transformation of Neural Additive Models

NAM learns a linear combination of MLPs corresponding to each feature.
Accordingly, the computational result of each MLP indicates the direct effect
of each feature on prediction. Owing to this characteristic, the
interpretations of NAM are straightforward. In the original work, NAM was
defined as follows:

$\hat{y}=\sum_{i=1}^{m}f_{i}\left(x_{i}\right)+b,$ (1)

where $x\in\mathbb{R}^{m}$ is the input feature,
$f_{i}:\mathbb{R}\rightarrow\mathbb{R}$ is the MLP corresponding to the
i${}^{\text{th}}$ feature, $b\in\mathbb{R}$ is the output bias, and $m$ is the
number of features. As shown in Equation 1, NAM is an extension of the
generalized additive models (GAM) (Hastie and Tibshirani 1987). The outputs of
all feature networks in NAM are scalars, but the outputs should be vectors to
effectively model the feature interactions. Therefore, we transform NAM as
follows:

$\mathcal{F}(x)=\sum_{i=1}^{m}f_{i}\left(x_{i}\right),$ (2)

where $f_{i}:\mathbb{R}\rightarrow\mathbb{R}^{k}$ is the MLP corresponding to
the i${}^{\text{th}}$ feature, similar to the existing NAM, but the output is
a $k$-dimensional vector. Since
$\mathcal{F}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{k}$ is a function that adds
all representation vectors of the input features, the output of $\mathcal{F}$
can be viewed as the sum of 1${}^{\text{st}}$-order features.

### Modeling High-order Feature Interactions

A simple way to model high-order feature interactions is to enumerate all
feature interactions (Rendle 2010) as follows:

$g_{t}\left(x\right)=\sum_{i_{1}=1}^{m}\sum_{i_{2}=i_{1}+1}^{m}\ldots\sum_{i_{t}=i_{t-1}+1}^{m}\prod_{l=1}^{t}x_{i_{l}},$
(3)

where $t$ is the order of feature interactions we want to model. This method
makes all possible combinations of $t^{\text{th}}$-order features using $m$
features. Simply enumerating all possible interactions is intuitive. However,
this method has $O(km^{t})$ time complexity, which is exponential time; where
$k$ is the embedding size for $x$. Therefore, this method shows very slow
training and inference times.

Another way to model feature interactions is to use multi-layer architecture
(Wang et al. 2017) defined as follows:

$\displaystyle g_{t}\left(x\right)=\left(g_{1}\left(x\right)\odot
g_{t-1}\left(x\right)\right)W_{t},$ (4) $\displaystyle
g_{1}\left(x\right)=xW_{1},$ (5)

where $W_{1}\in\mathbb{R}^{m\times k}$ and $W_{t>1}\in\mathbb{R}^{k\times k}$
are the trainable weights. We refer to this method CrossNet. CrossNet
multiplies $g_{1}(x)$ in every each layer to model high-order feature
interactions. For example, the 1${}^{\text{st}}$ layer $g_{1}(x)=xW_{1}$
indicates the weighted sum of 1${}^{\text{st}}$-order features. The
2${}^{\text{nd}}$ layer $g_{2}(x)=(g_{1}(x)\odot g_{1}(x))W_{2}$ is the
multiplication of two weighted sum of 1${}^{\text{st}}$-order features. In
other words, it indicates the weighted sum of 2${}^{\text{nd}}$-order
features. Its multi-layer structure is similar to MLP, but there is no
activation function. Therefore, CrossNet is a linear model, and $t$-layer
CrossNet can explicitly model $t^{\text{th}}$-order feature interactions. In
addition, CrossNet has $O(tk^{2})$ time complexity, which is significantly
faster than the enumerating method.

However, CrossNet has problems from the interpretability perspective. CrossNet
also models powers of features such as $x_{1}^{2}$, $x_{2}^{4}$, or
$x_{1}x_{2}^{3}$. Such powers of features may confuse us in interpreting
prediction results. For example, Figure 2(a) and Figures 2(b) show the heat
maps of 1${}^{\text{st}}$-order and 2${}^{\text{nd}}$-order features’
contribution values on the Insurance dataset learned by the CrossNet. In the
2${}^{\text{nd}}$-order features, human cannot know the meaning of squared
features such as $(sex=male)^{2}$ and $(smoker=yes)^{2}$. Therefore, powers of
features have no meaning from the interpretability perspective. In addition,
1${}^{\text{st}}$-order feature $(smoker=yes)$ has highly positive value, but
squared feature $(smoker=yes)^{2}$ has negative value. As can be seen in
Figure 2, powers of features can have different contributions with their
1${}^{\text{st}}$-order features. This phenomenon confuses people in
determining the actual effect of $(smoker=yes)$.

3 | California Housing | Insurance | House Prices | Bikeshare | Year
---|---|---|---|---|---
| R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute
2 LR | 0.630 ($\pm$0.001) | 0.441 ($\pm$0.002) | 0.780 ($\pm$0.003) | 0.548 ($\pm$0.014) | 0.830 ($\pm$0.031) | 0.594 ($\pm$0.036) | 0.362 ($\pm$0.012) | 0.238 ($\pm$0.012) | 0.247 ($\pm$0.001) | 0.169 ($\pm$0.001)
CrossNet (t=2) | 0.717 ($\pm$0.001) | 0.525 ($\pm$0.001) | 0.860 ($\pm$0.005) | 0.688 ($\pm$0.020) | 0.860 ($\pm$0.057) | 0.655 ($\pm$0.063) | 0.498 ($\pm$0.008) | 0.316 ($\pm$0.009) | 0.310 ($\pm$0.001) | 0.211 ($\pm$0.001)
XGBoost |  0.825 ($\pm$0.001) |  0.649 ($\pm$0.002) | 0.860 ($\pm$0.009) | 0.735 ($\pm$0.012) | 0.900 ($\pm$0.016) | 0.711 ($\pm$0.011) |  0.948 ($\pm$0.003) |  0.816 ($\pm$0.002) | 0.297 ($\pm$0.000) | 0.201 ($\pm$0.000)
MLP | 0.787 ($\pm$0.005) | 0.614 ($\pm$0.005) | 0.866 ($\pm$0.002) | 0.711 ($\pm$0.004) | 0.856 ($\pm$0.026) | 0.688 ($\pm$0.015) | 0.925 ($\pm$0.006) | 0.762 ($\pm$0.008) |  0.334 ($\pm$0.003) |  0.247 ($\pm$0.001)
EBM | 0.802 ($\pm$0.001) | 0.612 ($\pm$0.002) | 0.879 ($\pm$0.002) | 0.736 ($\pm$0.005) |  0.905 ($\pm$0.006) |  0.726 ($\pm$0.007) | 0.906 ($\pm$0.003) | 0.737 ($\pm$0.004) | 0.282 ($\pm$0.001) | 0.191 ($\pm$0.001)
NAM | 0.744 ($\pm$0.002) | 0.549 ($\pm$0.003) | 0.781 ($\pm$0.003) | 0.544 ($\pm$0.009) | 0.885 ($\pm$0.013) | 0.715 ($\pm$0.009) | 0.694 ($\pm$0.006) | 0.476 ($\pm$0.005) | 0.276 ($\pm$0.001) | 0.185 ($\pm$0.002)
NodeGAM | 0.734 ($\pm$0.004) | 0.540 ($\pm$0.004) | 0.776 ($\pm$0.003) | 0.534 ($\pm$0.004) | 0.897 ($\pm$0.008) | 0.732 ($\pm$0.007) | 0.696 ($\pm$0.005) | 0.476 ($\pm$0.005) | 0.275 ($\pm$0.001) | 0.188 ($\pm$0.001)
NodeGA2M | 0.808 ($\pm$0.004) | 0.620 ($\pm$0.005) | 0.879 ($\pm$0.004) | 0.732 ($\pm$0.010) | 0.887 ($\pm$0.015) | 0.712 ($\pm$0.024) | 0.911 ($\pm$0.004) | 0.740 ($\pm$0.004) | 0.309 ($\pm$0.001) | 0.213 ($\pm$0.001)
HONAM∗ (t=2) | 0.807 ($\pm$0.004) | 0.621 ($\pm$0.004) | 0.880 ($\pm$0.001) | 0.744 ($\pm$0.003) | 0.900 ($\pm$0.018) | 0.725 ($\pm$0.017) | 0.920 ($\pm$0.004) | 0.752 ($\pm$0.008) | 0.320 ($\pm$0.003) | 0.224 ($\pm$0.003)
HONAM (t=2) | 0.810 ($\pm$0.003) | 0.626 ($\pm$0.003) |  0.882 ($\pm$0.002) |  0.742 ($\pm$0.012) | 0.900 ($\pm$0.011) | 0.721 ($\pm$0.009) | 0.925 ($\pm$0.004) | 0.760 ($\pm$0.006) | 0.320 ($\pm$0.002) | 0.226 ($\pm$0.004)
3 |  |  |  |  |  |  |  |  |  |

Table 1: Evaluation on Regression Tasks

3 | FICO | Credit | Support2 | MIMIC III | Click
---|---|---|---|---|---
| AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC
LR | 0.753 ($\pm$0.001) | 0.723 ($\pm$0.002) | 0.978 ($\pm$0.012) | 0.800 ($\pm$0.045) | 0.783 ($\pm$0.009) | 0.577 ($\pm$0.020) | 0.761 ($\pm$0.012) | 0.317 ($\pm$0.031) | 0.617 ($\pm$0.002) | 0.620 ($\pm$0.002)
CrossNet (t=2) | 0.770 ($\pm$0.002) | 0.752 ($\pm$0.004) | 0.964 ($\pm$0.024) | 0.810 ($\pm$0.038) | 0.803 ($\pm$0.009) | 0.612 ($\pm$0.013) | 0.787 ($\pm$0.005) | 0.343 ($\pm$0.022) | 0.619 ($\pm$0.001) | 0.625 ($\pm$0.001)
XGBoost | 0.767 ($\pm$0.003) | 0.744 ($\pm$0.004) | 0.980 ($\pm$0.006) |  0.857 ($\pm$0.035) | 0.800 ($\pm$0.014) | 0.602 ($\pm$0.018) | 0.791 ($\pm$0.009) | 0.342 ($\pm$0.025) |  0.687 ($\pm$0.001) |  0.680 ($\pm$0.001)
MLP | 0.770 ($\pm$0.002) | 0.749 ($\pm$0.003) | 0.980 ($\pm$0.011) | 0.803 ($\pm$0.051) | 0.801 ($\pm$0.008) | 0.606 ($\pm$0.013) | 0.792 ($\pm$0.009) | 0.346 ($\pm$0.025) | 0.627 ($\pm$0.003) | 0.627 ($\pm$0.002)
EBM | 0.700 ($\pm$0.002) | 0.625 ($\pm$0.002) | 0.885 ($\pm$0.023) | 0.699 ($\pm$0.048) | 0.674 ($\pm$0.017) | 0.434 ($\pm$0.019) | 0.563 ($\pm$0.012) | 0.167 ($\pm$0.017) | 0.592 ($\pm$0.001) | 0.551 ($\pm$0.001)
NAM |  0.783 ($\pm$0.003) | 0.760 ($\pm$0.003) | 0.979 ($\pm$0.014) | 0.845 ($\pm$0.036) | 0.815 ($\pm$0.014) | 0.626 ($\pm$0.011) | 0.815 ($\pm$0.005) | 0.380 ($\pm$0.025) | 0.655 ($\pm$0.002) | 0.652 ($\pm$0.002)
NodeGAM | 0.781 ($\pm$0.002) |  0.761 ($\pm$0.002) | 0.980 ($\pm$0.012) | 0.849 ($\pm$0.040) | 0.814 ($\pm$0.012) | 0.626 ($\pm$0.010) | 0.813 ($\pm$0.006) | 0.375 ($\pm$0.017) | 0.643 ($\pm$0.002) | 0.643 ($\pm$0.001)
NodeGA2M | 0.780 ($\pm$0.002) | 0.760 ($\pm$0.003) |  0.982 ($\pm$0.011) | 0.846 ($\pm$0.038) | 0.812 ($\pm$0.012) | 0.624 ($\pm$0.008) | 0.816 ($\pm$0.009) | 0.374 ($\pm$0.023) | 0.641 ($\pm$0.002) | 0.639 ($\pm$0.002)
HONAM∗ (t=2) |  0.783 ($\pm$0.003) |  0.761 ($\pm$0.005) |  0.982 ($\pm$0.009) | 0.842 ($\pm$0.024) | 0.819 ($\pm$0.009) | 0.633 ($\pm$0.011) | 0.825 ($\pm$0.003) | 0.395 ($\pm$0.023) | 0.667 ($\pm$0.004) | 0.663 ($\pm$0.003)
HONAM (t=2) | 0.782 ($\pm$0.002) | 0.760 ($\pm$0.004) | 0.981 ($\pm$0.012) | 0.838 ($\pm$0.026) |  0.823 ($\pm$0.011) |  0.640 ($\pm$0.009) |  0.826 ($\pm$0.005) |  0.399 ($\pm$0.024) | 0.670 ($\pm$0.002) | 0.664 ($\pm$0.003)
3 |  |  |  |  |  |  |  |  |  |

Table 2: Evaluation on Classification Tasks

### Novel Feature Crossing for High Interpretability

We would like to make a novel feature interaction method, which is fast and
powers of features are eliminated. To achieve this goal, we found a recursion
form of the enumerating method.

Theorem 1 We develop a method for recursively eliminating powers of features.
First, we make $t^{\text{th}}$-order features with squared features as
follows:

$\left(\sum_{i=1}^{m}x_{i}\right)\left(\sum_{i_{1}>\ldots>i_{t-1}}^{m}\prod_{l=1}^{t-1}x_{i_{l}}\right).$
(6)

Equation 6 indicates the multiplication of the sum of
$(t-1)^{\text{th}}$-order features and the sum of 1${}^{\text{st}}$-order
features. Therefore, there are $t^{\text{th}}$-order features with squared
features. To remove the squared features, we explicitly subtract them.
However, unexpected features also eliminated in this process. Therefore, we
repeat subtracting powers of features and adding back unexpectedly eliminated
features as follows:

$\displaystyle\begin{split}&\left(\sum_{i=1}^{m}x_{i}\right)\left(\sum_{i_{1}>\ldots>i_{t-1}}^{m}\prod_{l=1}^{t-1}x_{i_{l}}\right)\\\
&-\left(\sum_{i=1}^{m}x_{i}^{2}\right)\left(\sum_{i_{1}>\ldots>i_{t-2}}^{m}\prod_{l=1}^{t-2}x_{i_{l}}\right)\\\
&+\left(\sum_{i=1}^{m}x_{i}^{3}\right)\left(\sum_{i_{1}>\ldots>i_{t-3}}^{m}\prod_{l=1}^{t-3}x_{i_{l}}\right)\\\
&-...\end{split}$ (7)

Now, we can make this iteration to recursion form as follows:

$\displaystyle
g_{t}\left(x\right)=\frac{1}{t}\sum_{i=1}^{t}\left(-1\right)^{i+1}g_{1}\left(x^{i}\right)g_{t-i}\left(x\right),$
(8) $\displaystyle g_{1}\left(x\right)=\sum_{i=1}^{m}x_{i},$ (9)
$\displaystyle g_{0}\left(x\right)=1.$ (10)

Since, same interaction values are multiplied $t$ times in Equation 7, we
divide Equation 8 by $t$. Therefore, Equation 8 is equivalent to Equation 3,
but our method has $O(ktm)$ time complexity, which is linear time when using
dynamic programming. The proof of Theorem 1 can be found in Appendix.

The final prediction process of HONAM is defined as follows:

$\mathcal{G}\left(x\right)=\left(\overset{t}{\underset{i=1}{\mathbin{\|}}}g_{i}\left(x\right)\right)W+b,$
(11) $\hat{y}=\mathcal{G}\left(\mathcal{F}\left(x\right)\right),$ (12)

where $\mathbin{\|}$ denotes the concatenation operator, $W$ is the output
weight, $b$ is the output bias, $\mathcal{G}$ is a feature interaction method
such as our proposed method or CrossNet, and $\mathcal{F}$ is the transformed
NAM defined in Equation 2. Figure 1 shows the architecture of HONAM. We
believe that our feature interaction method is a very useful choice for
interpretable machine learning.

### Relationship with Higher-order Factorization Machines

Although we started this study from a different perspective than HOFM (Blondel
et al. 2016), our feature interaction method is equivalent to the alternative
recursion form of HOFM. Therefore, HOFM can be viewed as a special case of
HONAM. Suppose that each feature network is linear, and the output weight is
an all-one matrix. HONAM is equivalent to HOFM. Furthermore, HONAM overcomes
the limitations of the existing FM/HOFM, feature interaction methods, and NAM.
First, vanilla HOFM in the original work (Rendle 2010) has the time complexity
of $O(km^{t})$, but our feature interactions method has the time complexity of
$O(tkm)$ when using dynamic programming as mentioned above; where $k$ is the
embedding size. Second, since FM/HOFM and other feature interaction methods
are essentially linear models, their expressive powers are limited. HONAM can
learn non-linear functions and has high expressive power because it utilizes
neural networks, and is still interpretable. Third, NAM can only model
1${}^{\text{st}}$-order feature interactions, whereas HONAM can model high-
order feature interactions. Thus, HONAM has higher predictive performance than
NAM and can provide high-quality interpretations.

## Learning Sharp Shape Functions

(a) ExU

(b) ExpDive

Figure 3: Weight shape functions of ExU and ExpDive units

(a) Effects of first-order features

(b) Effects of second-order features

Figure 4: Local interpretations for the FICO dataset

(a) Effects of a first-order feature

(b) Effects of a second-order feature

Figure 5: Global interpretations for the FICO dataset

Since sharp jumps are frequently found in real-world data, making machine-
learning models that can learn sharp jumps is a significant challenge in
interpretable models (Agarwal et al. 2020; Chang et al. 2021). For example,
(Agarwal et al. 2020) found that the average house prices highly jump in San
Francisco and Los Angeles in the California Housing dataset. In addition,
(Chang et al. 2021) found that patients’ mortality risks significantly dropped
at PFratio=332 because they imputed missing values by the mean 332 (healthy
patients commonly do not check PFratio). However, neural networks struggle to
capture such sharp jumps (Agarwal et al. 2020) because of the bias towards
smooth functions (Arpit et al. 2017; Rahaman et al. 2019).

To address this problem, NAM uses the Exp Centered (ExU) unit defined as
follows:

$\sigma\left(\left(x-b\right)\times\exp\left(W\right)\right).$ (13)

where $\sigma$ is the ReLU-n function, $x$ is the input feature, and $W$ and
$b$ are the trainable weight and bias, respectively. The ExU unit can learn
sharp jumps because its outputs are significantly changed with small input
changes. However, the ExU unit works well under certain circumstances. In this
section, we argue the limitations of the ExU unit and propose a novel hidden
unit named the ExpDive unit, which is a more effective way to learn sharp
shape functions.

Theorem 2 Since $exp(W)$ is always positive, the ExU unit does not perform
well when an input feature is negative. Suppose an ExU unit with ReLU
activation, which is the bias term omitted, is as follows:

$ReLU\left(x\times\exp\left(W\right)\right).$ (14)

Equation 14 is always zero when $x$ is negative. In other words, the gradient
for the back-propagation is zero. Therefore, to propagate the gradient when
$x$ is negative, the bias $b$ must be a negative number less than $x$.
Therefore, the ExU unit often fails to learn. The proof of Theorem 2 can be
found in Appendix B.

To overcome this problem, we propose the ExpDive unit. Since the ExpDive unit
exponentially diverges in both the negative and positive ranges, it is not
only suitable for learning sharp shape functions, but also the gradient
vanishing problem does not occur. The ExpDive unit is defined as follows:

$\sigma\left(\left(x-b\right)\times\left(\exp\left(W\right)-\exp\left(-W\right)\right)\right).$
(15)

Figure 3 shows the shape plots of the weights of the ExU and ExpDive units.
Contrary to the ExU unit, the ExpDive unit can propagate gradients if the
weight is negative, even if the input $x$ is negative.

## Experiments

We conduct extensive experiments to demonstrate the effectiveness of HONAM. In
addition, we conduct experiments to learn about sharp jumps and fair AI.
Descriptions of the experimental setup and datasets is provided in Appendix.

### Effectiveness of Modeling Feature Interactions

Although we mainly focus on the interpretability of HONAM in this study, we
intend to show that feature interaction modeling is also effective for
predictive performance. To this end, we compare HONAM with other machine-
learning models using various regression and classification datasets. Table 1
and Table 2 show the performance of the different comparative models on the
regression and classification tasks, respectively. HONAM∗ indicates HONAM with
CrossNet, and HONAM indicates HONAM with the proposed interaction method in
the tables. HONAM significantly outperforms EBM, NAM and NodeGAM in both
regression and classification tasks. It shows that feature interaction
modeling is effective for improving predictive performances. HONAM also
outperforms NodeGA${}^{\text{2}}$M. It demonstrates that the fully neural
network method (HONAM) is more effective than the tree-based method (NodeGAM).
In addition, our HONAM outperforms or has comparable performances to full-
complexity methods, MLP and XGBoost. Some GAM-based methods show better
performance than MLP in several datasets. This may be because the
distributions of each feature were different. Dealing with features with
different distributions using shared trainable parameters may disturb training
(Yan et al. 2022). Because NAM-family methods deal with features using
different trainable parameters, this problem is alleviated. We also conduct an
interaction ablation study in Appendix.

Figure 6: Learned shape functions on the COMPAS dataset

### Interpretations

In HONAM, since each input feature is handled independently by the
corresponding network, the computational result of each feature network
indicates the direct effect of a feature on the output value. Owing to this
characteristic, it is easy to interpret the predictions of HONAM. We can
acquire local interpretations by forward propagation, and by repeating the
process of local interpretation, we can also acquire global interpretations of
HONAM. In this study, we visualize local interpretations of the
1${}^{\text{st}}$-order and 2${}^{\text{nd}}$-order features and visualize
global interpretations to understand the global behavior of HONAM. Note that
interpretations of high-order features, such as the 3${}^{\text{rd}}$-order or
4${}^{\text{th}}$-order, can also be obtained using HONAM.

Figure 4(a) and Figure 4(b) show the visualization of the local
interpretations for the 1${}^{\text{st}}$-order and 2${}^{\text{nd}}$-order
features in the FICO dataset, respectively. The figures show the actual
contribution values of features to the prediction. We can get these values via
the forward propagation of HONAM. For a more detailed explanation, we refer
the representation vectors of each input feature ${x_{1},x_{2},…,x_{m}}$
computed by the feature networks as ${r_{1},r_{2},…,r_{m}}$, where $m$ is the
number of features. To obtain 1st-order interpretation, the contribution value
of feature $i$ is computed by $r_{i}\times w_{0:k}$, where $w$ is the output
matrix and $k$ is the size of the representation vector. Similarly, the
contribution value of 2nd-order feature $(x_{i}*x_{j})$ can be computed as
$(r_{i}*r_{j})\times w_{k:2*k}$. Because the FICO dataset is a binary
classification task, features with positive effects increase the probability
that the model predicts this data sample as positive, whereas features with
negative effects decrease the probability. For a straightforward
understanding, we visualize the interpretation of the 2${}^{\text{nd}}$-order
features using a heat map. In Figure 4(b), red cells indicate positive
features, and blue cells indicate negative features.

In addition, we visualize the global interpretation of HONAM. Figure 5(a) and
Figure 5(b) show the global interpretation of 1${}^{\text{st}}$-order feature
ExternalRiskEstimate and 2${}^{\text{nd}}$-order feature ExternalRistEstimate
$\times$ NumInstallTradesWBalance in the FICO dataset, respectively. In Figure
5(a), the red bars indicate the densities of training samples. In the
1${}^{\text{st}}$-order interpretation, ExternalRistEstimate tends to have a
more positive effect as increases. However, the 2${}^{\text{nd}}$-order
interpretation exhibits a different tendency. In the 2${}^{\text{nd}}$-order
interpretation, even if ExternalRistEstimate has a high value, it has a
negative effect if NumInstallTradesWBalance is small. This is interesting
because this tendency cannot be found in 1${}^{\text{st}}$-order interactions
and can only be found in 2${}^{\text{nd}}$-order interactions. This
demonstrates that HONAM can provide more detailed interpretations than NAM and
why modeling high-order feature interaction is important in interpretable
machine learning.

### HONAM for Fair AI

3 | AUROC | AUPRC
---|---|---
HONAM | 0.832 ($\pm$0.001) | 0.798 ($\pm$0.001)
HONAM (remove sex bias) | 0.829 ($\pm$0.001) | 0.794 ($\pm$0.001)
HONAM (remove race bias) | 0.817 ($\pm$0.002) | 0.787 ($\pm$0.001)
HONAM (remove sex and race biases) | 0.816 ($\pm$0.001) | 0.784 ($\pm$0.002)
3 |  |

Table 3: Predictive performances of biased and unbiased HONAMs

3 African American vs | Overall | Asian | Caucasian | Hispanic | Native American | Other
---|---|---|---|---|---|---
HONAM | 0.349 | 0.327 | 0.398 | 0.293 | 0.358 | 0.174
Unbiased HONAM | 0.524 | 0.383 | 0.547 | 0.458 | 0.823 | 0.502
3 |  |  |  |  |  |

Table 4: Disparate Impact on Race Features

3 | California Housing | Insurance | House Prices
---|---|---|---
| Adj R-squared | Adj R-absolute | Adj R-squared | Adj R-absolute | Adj R-squared | Adj R-absolute
2 HONAM${}_{\text{ExU}}$ | 0.782 ($\pm$0.011) | 0.594 ($\pm$0.015) |  0.878 ($\pm$0.004) |  0.735 ($\pm$0.018) | 0.847 ($\pm$0.023) | 0.613 ($\pm$0.037)
HONAM${}_{\text{ExpDive}}$ |  0.797 ($\pm$0.012) |  0.612 ($\pm$0.015) |  0.878 ($\pm$0.004) | 0.727 ($\pm$0.016) |  0.850 ($\pm$0.026) |  0.630 ($\pm$0.022)
3 |  |  |  |  |  |

Table 5: Predictive performance comparison of different units on regression
tasks

3 | FICO | Credit | MIMIC IV
---|---|---|---
| AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC
2 HONAM${}_{\text{ExU}}$ | 0.776 ($\pm$0.006) | 0.756 ($\pm$0.009) | 0.973 ($\pm$0.005) | 0.842 ($\pm$0.019) |  0.823 ($\pm$0.001) |  0.074 ($\pm$0.002)
HONAM${}_{\text{ExpDive}}$ |  0.780 ($\pm$0.004) |  0.762 ($\pm$0.006) |  0.975 ($\pm$0.006) |  0.849 ($\pm$0.014) |  0.823 ($\pm$0.001) |  0.074 ($\pm$0.001)
3 |  |  |  |  |  |

Table 6: Predictive performance comparison of different units on
classification tasks

Recently, fair AI has also attracted attention as a crucial machine learning
research topic. Real-world datasets may have unexpected biases, and many
researchers are working hard to eliminate such biases in machine-learning
models. Fortunately, owing to the modularity of NAM (Agarwal et al. 2020),
removing biases from trained HONAM is easy. We show that a HONAM trained on
the COMPAS recidivism dataset is biased and introduce a method for eliminating
biases.

Figure 6 shows the visualizations of the feature shape plots learned by HONAM
for the COMPAS dataset. As shown in the figure, the trained HONAM is highly
biased towards race features. These biases are undesirable and may be
ethically problematic. Since HONAM handles each feature independently, it is
easy to remove these biases. For example, we can drop out all feature networks
corresponding to race features to eliminate racial bias from the trained
HONAM. Consequently, the effect of the race field is set to 0. Therefore, we
can make a fair HONAM by sacrificing a slight predictive performance through
this strategy, as shown in Table 3, which shows the predictive performance of
biased and unbiased HONAMs. In addition, Table 4 shows the disparate impacts
on race features for biased and unbiased HONAMs. The disparate impacts are
significantly improved on all race features.

(a) ExU

(b) ExpDive

Figure 7: Learned shape functions of different units on the synthetic
classification task

### Comparison of Different Units

In this section, we demonstrate that the proposed ExpDive unit is more
effective in learning sharp shape functions than the ExU unit using the
synthetic classification dataset. Subsequently, we compare ExU, and ExpDive
units to determine whether these units work well on real-world datasets.

Figure 7 shows the visualization of the feature shape plots learned by each
unit for the synthetic classification dataset. Consistent with our theory, the
ExU unit was not properly trained in the negative range. In addition, we can
observe that the ExpDive unit trained the sharpest function on the dataset.
The process of creating the synthetic dataset is described in Appendix, and
more information on the experiment is also provided in Appendix.

Table 5 and Table 6 show the predictive performance of HONAMs with different
units. The ExpDive unit outperforms the ExU unit. It demonstrates that our
ExpDive unit is more effective than the ExU unit. We anticipate that the
ExpDive unit will work well in fields where sharp jumps frequently occur, such
as stock price predictions.

## Conclusion

In this study, we proposed HONAM, a novel interpretable machine learning
method that can model arbitrary orders of feature interactions. To develop the
HONAM, we redefined the existing NAM and applied CrossNet to model high-order
feature interactions. Furthermore, we proposed a novel feature interaction
method for high interpretability. We also proposed the ExpDive unit, which
overcomes the limitations of the ExU unit, to effectively learn sharp shape
functions. We conducted extensive experiments to examine the effectiveness of
HONAM. The experimental results demonstrated that feature interaction modeling
is effective for enhancing predictive performance, and HONAM outperforms or
has comparable performance to MLP and XGBoost. We visualized both the local
and global interpretations of HONAM for the 1${}^{\text{st}}$-order and
2${}^{\text{nd}}$-order features. These interpretations have shown that it is
important to model high-order feature interactions in terms of
interpretability, and HONAM can capture high-order co-relationships that NAM
cannot. In addition, we showed that it is possible to create a fair HONAM by
eliminating the biases of HONAM. HONAM has achieved high-performance
interpretable machine-learning models required by high-stakes domains. We
expect HONAM to be widely used in various domains, such as healthcare and
mortality prediction.

## References

  * Agarwal et al. (2020) Agarwal, R.; Melnick, L.; Frosst, N.; Zhang, X.; Lengerich, B.; Caruana, R.; and Hinton, G. 2020.  Neural Additive Models: Interpretable Machine Learning with Neural Nets.  In _Proceedings of the 34th International Conference on Neural Information Processing Systems_.
  * Arık and Pfister (2021) Arık, S. O.; and Pfister, T. 2021.  TabNet: Attentive Interpretable Tabular Learning.  In _AAAI_.
  * Arpit et al. (2017) Arpit, D.; Jastrzebski, S.; Ballas, N.; Krueger, D.; Bengio, E.; Kanwal, M. S.; Maharaj, T.; Fischer, A.; Courville, A.; Bengio, Y.; et al. 2017.  A Closer Look at Memorization in Deep Networks.  In _International conference on machine learning_.
  * Blondel et al. (2016) Blondel, M.; Fujino, A.; Ueda, N.; and Ishihata, M. 2016.  Higher-Order Factorization Machines.  In _Proceedings of the 29th International Conference on Neural Information Processing Systems_.
  * Caruana (1997) Caruana, R. 1997.  Multitask Learning.  _Machine learning_ , 28: 41–75.
  * Caruana et al. (2015) Caruana, R.; Lou, Y.; Gehrke, J.; Koch, P.; Sturm, M.; and Elhadad, N. 2015.  Intelligible Models for HealthCare: Predicting Pneumonia Risk and Hospital 30-day Readmission.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_.
  * Chang, Caruana, and Goldenberg (2022) Chang, C.-H.; Caruana, R.; and Goldenberg, A. 2022.  NODE-GAM: Neural Generalized Additive Model for Interpretable Deep Learning.  In _International Conference on Learning Representations_.
  * Chang et al. (2021) Chang, C.-H.; Tan, S.; Lengerich, B.; Goldenberg, A.; and Caruana, R. 2021.  How Interpretable and Trustworthy are GAMs?  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_.
  * Chen and Guestrin (2016) Chen, T.; and Guestrin, C. 2016.  XGBoost: A Scalable Tree Boosting System.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.
  * Cheng et al. (2016) Cheng, H.-T.; Koc, L.; Harmsen, J.; Shaked, T.; Chandra, T.; Aradhye, H.; Anderson, G.; Corrado, G.; Chai, W.; Ispir, M.; et al. 2016.  Wide & Deep Learning for Recommender Systems.  In _Proceedings of the 1st Workshop on Deep Learning for Recommender Systems_.
  * Cheng, Shen, and Huang (2020) Cheng, W.; Shen, Y.; and Huang, L. 2020.  Adaptive Factorization Network: Learning Adaptive-Order Feature Interactions.  In _Proceedings of the AAAI Conference on Artificial Intelligence_.
  * Dal Pozzolo (2015) Dal Pozzolo, A. 2015.  Adaptive Machine Learning for Credit Card Fraud Detection.  _PhD Thesis, Université libre de Bruxelles_.
  * Guo et al. (2017) Guo, H.; Tang, R.; Ye, Y.; Li, Z.; and He, X. 2017.  DeepFM: A Factorization-Machine Based Neural Network for CTR Prediction.  In _Proceedings of the 26th International Joint Conference on Artificial Intelligence_.
  * Hastie and Tibshirani (1987) Hastie, T.; and Tibshirani, R. 1987.  Generalized Additive Models: Some Applications.  _Journal of the American Statistical Association_ , 82: 371–386.
  * Johnson et al. (2021) Johnson, A.; Bulgarelli, L.; Pollard, T.; Horng, S.; Celi, L. A.; and Mark, R. 2021\.  MIMIC-IV.  _PhysioNet_.  Https://doi.org/10.13026/s6n6-xd98.
  * Johnson et al. (2016) Johnson, A. E.; Pollard, T. J.; Shen, L.; Lehman, L.-w. H.; Feng, M.; Ghassemi, M.; Moody, B.; Szolovits, P.; Anthony Celi, L.; and Mark, R. G. 2016.  MIMIC-III, A Freely Accessible Critical Care Database.  _Scientific data_ , 3(1): 1–9.
  * Kim, Choi, and Kim (2022) Kim, M.; Choi, H.-S.; and Kim, J. 2022.  Explicit Feature Interaction-aware Graph Neural Networks.  _arXiv preprint arXiv:2204.03225_.
  * Kim, Lee, and Kim (2020) Kim, M.; Lee, S.; and Kim, J. 2020.  Combining Multiple Implicit-Explicit Interactions for Regression Analysis.  In _2020 IEEE International Conference on Big Data (Big Data)_.
  * Lian et al. (2018) Lian, J.; Zhou, X.; Zhang, F.; Chen, Z.; Xie, X.; and Sun, G. 2018.  xDeepFM: Combining Explicit and Implicit Feature Interactions for Recommender Systems.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_.
  * Lou, Caruana, and Gehrke (2012) Lou, Y.; Caruana, R.; and Gehrke, J. 2012.  Intelligible Models for Classification and Regression.  In _Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.
  * Lou et al. (2013) Lou, Y.; Caruana, R.; Gehrke, J.; and Hooker, G. 2013.  Accurate Intelligible Models with Pairwise Interactions.  In _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.
  * Lundberg and Lee (2017) Lundberg, S. M.; and Lee, S.-I. 2017.  A Unified Approach to Interpreting Model Predictions.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_.
  * Nori et al. (2019) Nori, H.; Jenkins, S.; Koch, P.; and Caruana, R. 2019.  InterpretML: A Unified Framework for Machine Learning Interpretability.  _arXiv preprint arXiv:1909.09223_.
  * Popov, Morozov, and Babenko (2019) Popov, S.; Morozov, S.; and Babenko, A. 2019.  Neural Oblivious Decision Ensembles for Deep Learning on Tabular Data.  _arXiv preprint arXiv:1909.06312_.
  * Potts (1999) Potts, W. J. 1999.  Generalized Additive Neural Networks.  In _Proceedings of the fifth ACM SIGKDD international conference on Knowledge discovery and data mining_.
  * Rahaman et al. (2019) Rahaman, N.; Baratin, A.; Arpit, D.; Draxler, F.; Lin, M.; Hamprecht, F.; Bengio, Y.; and Courville, A. 2019.  On the Spectral Bias of Neural Networks.  In _International Conference on Machine Learning_.
  * Rendle (2010) Rendle, S. 2010.  Factorization Machines.  In _2010 IEEE International Conference on Data Mining_.
  * Ribeiro, Singh, and Guestrin (2016) Ribeiro, M. T.; Singh, S.; and Guestrin, C. 2016.  ”Why Should I Trust You?” Explaining the Predictions of Any Classifier.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_.
  * Rudin (2019) Rudin, C. 2019.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ , 1: 206–215.
  * Wang et al. (2017) Wang, R.; Fu, B.; Fu, G.; and Wang, M. 2017.  Deep & Cross Network for Ad Click Predictions.  In _Proceedings of the ADKDD’17_.
  * Xiao et al. (2017) Xiao, J.; Ye, H.; He, X.; Zhang, H.; Wu, F.; and Chua, T.-S. 2017.  Attentional Factorization Machines: Learning the Weight of Feature Interactions via Attention Networks.  In _Proceedings of the 26th International Joint Conference on Artificial Intelligence_.
  * Yan et al. (2022) Yan, B.; Wang, P.; Zhang, K.; Li, F.; Xu, J.; and Zheng, B. 2022.  APG: Adaptive Parameter Generation Network for Click-Through Rate Prediction.  _arXiv preprint arXiv:2203.16218_.

##  Appendix A Proof of Theorem 1

We prove Theorem 1 using mathematical induction.

Definition A.1. The sum of 1${}^{\text{st}}$-order feature interactions in
which powers of features are completely removed can be defined as follows:

$\displaystyle\sum_{i_{1}=1}^{m}\sum_{i_{2}=i_{1}+1}^{m}\ldots\sum_{i_{t}=t_{t-1}+1}^{m}\prod_{l=1}^{t}x_{i_{l}}$
$\displaystyle=:\sum_{i_{1}>\ldots>i_{t}}^{m}\prod_{l=1}^{t}x_{i_{l}}$ .

Therefore, the proof must satisfy Definition A.1.

Lemma A.2. When the interaction order $t=2$, Theorem 1 is as follows:

$\displaystyle g_{2}\left(x\right)$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{2}\left(-1\right)^{i+1}g_{1}\left(x^{i}\right)g_{t-i}\left(x\right)$
$\displaystyle=\frac{1}{2}\left\\{g_{1}\left(x\right)g_{1}\left(x\right)-g_{1}\left(x^{2}\right)g_{0}\left(x\right)\right\\}$
$\displaystyle=\frac{1}{2}\left\\{\left(\sum_{i=1}^{m}x_{i}\right)\left(\sum_{i=1}^{m}x_{i}\right)-\sum_{i=1}^{m}x^{2}_{i}\right\\}$
$\displaystyle=\sum_{i_{1}=1}^{m}\sum_{i_{2}=i}^{m}x_{i_{1}}x_{i_{2}}.$

Therefore, Theorem 1 holds when $t=2$ because Lemma A.1 satisfies Definition
A.1.

Assumption A.3. Then, we assume that Theorem 1 is true when $t=k-1$ as
follows:

$\displaystyle g_{k-1}\left(x\right)$
$\displaystyle=\frac{1}{k-1}\sum_{i=1}^{k-1}\left(-1\right)^{i+1}g_{1}\left(x^{i}\right)g_{k-i-1}\left(x\right)$
$\displaystyle=\sum_{i_{1}>\ldots>i_{k-1}}^{m}\prod_{l=1}^{k-1}x_{i_{l}}.$

Proof A.4. The sum of $k^{\text{th}}$-order interactions that include powers
of features is defined as follows:

$\displaystyle\left(\sum_{i=1}^{m}x_{i}\right)\left(\sum_{i_{1}>\ldots>i_{k-1}}^{m}\prod_{l=1}^{k-1}x_{i_{l}}\right)$
$\displaystyle=k\sum_{i_{1}>\ldots>i_{k}}^{m}\prod_{l=1}^{k}x_{i_{l}}+\sum_{i=1}^{m}x_{i}^{2}\underset{i\neq
j_{*}}{\sum_{j_{1}>\ldots>j_{k-1}}^{m}}\prod_{l=1}^{k-1}x_{j_{l}}$
$\displaystyle=k\sum_{i_{1}>\ldots>i_{k}}^{m}\prod_{l=1}^{k}x_{i_{l}}+\left(\sum_{i=1}^{m}x_{i}^{2}\right)\left(\sum_{j_{1}>\ldots>j_{k-1}}^{m}\prod_{l=1}^{k-1}x_{j_{l}}\right)$
$\displaystyle\hskip 83.0pt-\sum_{i=1}^{m}x_{i}^{3}\underset{i\neq
j_{*}}{\sum_{j_{1}>\ldots>j_{k-2}}^{m}}\prod_{l=1}^{k-2}x_{j_{l}}$
$\displaystyle=k\sum_{i_{1}>\ldots>i_{k}}^{m}\prod_{l=1}^{k}x_{i_{l}}+\left(\sum_{i=1}^{m}x_{i}^{2}\right)\left(\sum_{j_{1}>\ldots>j_{k-1}}^{m}\prod_{l=1}^{k-1}x_{j_{l}}\right)$
$\displaystyle\hskip
83.0pt-\left(\sum_{i=1}^{m}x_{i}^{3}\right)\left(\sum_{j_{1}>\ldots>j_{k-2}}^{m}\prod_{l=1}^{k-2}x_{j_{l}}\right)+\ldots$
$\displaystyle=k\sum_{i_{1}>\ldots>i_{k}}^{m}\prod_{l=1}^{k}x_{i_{l}}$
$\displaystyle\hskip
50.0pt+\sum_{i=1}^{k-1}\left(-1\right)^{i+1}\left(\sum_{j=1}^{m}x_{j}^{i}\right)\left(\sum_{j_{1}>\ldots>j_{k-i}}^{m}\prod_{l=1}^{k-i}x_{j_{l}}\right).$

Therefore, eliminating powers of features from the above equation can be
defined as follows:

$\displaystyle\left(\sum_{i=1}^{m}x_{i}\right)\left(\sum_{i_{1}>\ldots>i_{k-1}}^{m}\prod_{l=1}^{k-1}x_{i_{l}}\right)$
$\displaystyle\hskip
50.0pt-\sum_{i=1}^{k-1}\left(-1\right)^{i+1}\left(\sum_{j=1}^{m}x_{j}^{i}\right)\left(\sum_{j_{1}>\ldots>j_{k-i}}^{m}\prod_{l=1}^{k-i}x_{j_{l}}\right)$
$\displaystyle=\sum_{i=1}^{k}\left(-1\right)^{i+1}\left(\sum_{j=1}^{m}x_{j}^{i}\right)\left(\sum_{j_{1}>\ldots>j_{k-i}}^{m}\prod_{l=1}^{k-i}x_{i_{l}}\right)$
$\displaystyle=k\sum_{i_{1}>\ldots>i_{k}}^{m}\prod_{l=1}^{k}x_{i_{l}}.$

Subsequently, we remove duplicated interactions by dividing $k$ as follows:

$\displaystyle\sum_{i_{1}>\ldots>i_{k}}^{t}\prod_{l=1}^{k}x_{i_{l}}$
$\displaystyle\hskip
25.0pt=\frac{1}{k}\sum_{i=1}^{k}\left(-1\right)^{i+1}\left(\sum_{j=1}^{m}x_{j}^{i}\right)\left(\sum_{j_{1}>\ldots>j_{k-i}}^{m}\prod_{l=1}^{k-i}x_{i_{l}}\right).$

Finally, we can obtain the recursion form of interactions by using Lemma A.2
and Assumption A.3 as follows:

$\displaystyle\sum_{i_{1}>\ldots>i_{k}}^{t}\prod_{l=1}^{k}x_{i_{l}}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}\left(-1\right)^{i+1}\left(\sum_{j=1}^{m}x_{j}^{i}\right)\left(\sum_{j_{1}>\ldots>j_{k-i}}^{m}\prod_{l=1}^{k-i}x_{i_{l}}\right)$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}\left(-1\right)^{i+1}g_{1}(x^{i})g_{k-i}(x).$

Consequently, Theorem 1 is true because it satisfies Definition A.1.

$\square$

##  Appendix B Proof of Theorem 2

ExU units often fail to learn when using the ReLU activation, as shown below.
The $w$ with an exponential function always greater than zero as follows:

$\forall w\in\mathbb{R}:\exp\left(w\right)>0.$

Therefore, the state of an ExU unit with the ReLU activation can be expressed
as follows:

$ReLU\left(x\times\exp\left(w\right)\right)=\begin{cases}x\times\exp\left(w\right)\hskip
10.0ptif\;x>0\\\ 0\hskip 55.0ptotherwise\end{cases}.$

Consequently, the ExU units cannot learn the negative range of $x$ because the
back-propagation gradients are zero as
$ReLU\left(x\times\exp\left(w\right)\right)$ is always zero when $x\leq 0$.

In addition, we can suppose activation functions other than the ReLU. ExU
units are also limited to other activation functions, as shown below. The $w$
with an exponential function always greater than zero is as follows:

$\forall w\in\mathbb{R}:\exp\left(w\right)>0.$

Therefore, the state of a linear transformation of an ExU unit is as follows:

$\begin{cases}x\times\exp\left(w\right)\geq 0\hskip 10.0ptif\;x\geq 0\\\
x\times\exp\left(w\right)<0\hskip 10.0ptif\;x<0\end{cases}.$

When using ExU units, positive input $x$ always has a positive effect, and
negative input $x$ always has a negative effect.

$\square$

##  Appendix C Datasets

3 | # samples | # features | positive rate | task
---|---|---|---|---
2 California Housing | 20,640 | 8 | - | regression
Insurance | 1,338 | 6 | - | regression
House Prices | 1,460 | 60 | - | regression
Bikeshare | 17,389 | 16 | - | regression
Year | 515,345 | 90 | - | regression
FICO | 10,459 | 23 | 0.478 | classification
Credit | 284,807 | 29 | 0.002 | classification
Support2 | 9,105 | 29 | 0.259 | classification
MIMIC III | 27,348 | 57 | 0.098 | classification
MIMIC IV | 451,608 | 8 | 0.018 | classification
COMPAS | 11,742 | 4 | 0.437 | classification
Click | 1M | 11 | 0.500 | classification
3 |  |  |  |

Table 7: Datasets statistics

We use the following 12 public datasets for our experiments:

  * •

The California Housing dataset contains the 1990 California census
information. The goal of this dataset is to predict the median values of
houses in given California districts.
https://www.kaggle.com/datasets/camnugent/california-housing-prices

  * •

The Insurance dataset contains personal health information. The goal of this
dataset is to predict individual medical costs billed by health insurance.
https://www.kaggle.com/datasets/mirichoi0218/insurance

  * •

The House Prices dataset contains housing information such as location and the
total number of rooms. The goal of this dataset is to predict the sale prices
of each house. https://www.kaggle.com/competitions/house-prices-advanced-
regression-techniques

  * •

The Bikeshare dataset contains the counts of rental bikes with the
corresponding weather and seasonal features.
https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset

  * •

The Year dataset contains the features of songs ranging from 1922 to 2011. The
goal of this dataset is to predict the release year of a song.
https://archive.ics.uci.edu/ml/datasets/yearpredictionmsd

  * •

The FICO dataset contains information from customers who requested credit
lines. The goal of this dataset is to predict whether consumers who opened
credit accounts have been past due for more than 90 days in the last 24
months. https://community.fico.com/s/explainable-machine-learning-challenge

  * •

The Credit dataset contains completely de-identified features. The goal of
this dataset is to detect credit fraud. Please see the original work (Dal
Pozzolo 2015).

  * •

The Support2 dataset contains patients information from study to understand
prognoses preferences outcomes and risks of treatment.
https://biostat.app.vumc.org/wiki/Main/SupportDesc

  * •

The MIMIC III and MIMIC IV datasets are large-scale datasets containing a
variety of hospital information, such as hospitalization and prescription
information. Although there are many sub-tasks in this dataset, we focus on
mortality prediction (Johnson et al. 2016, 2021).

  * •

The COMPAS dataset contains personal and criminal information. The goal of
this dataset is to predict whether the criminal risks are low or high. We set
the target value to 0 if text_score is Low and 1 if text_score is Medium or
High as same as the original work. https://github.com/propublica/compas-
analysis

  * •

The Click dataset contains advertising information. The goal of this dataset
is to predict whether a user clicks an ad or not. We extracted positive and
negative samples 500K as same as the previous work, respectively (Popov,
Morozov, and Babenko 2019). https://www.kaggle.com/c/kddcup2012-track2

### Preprocessing

For categorical features, we do ordinary encoding to reduce memory usage. For
continuous features, we do standard scaling with mean and standard deviation.
Then we do quantile transformation for all features. When fitting quantile
transformation we add small gaussian nose. Adding noise is crucial to have
mean and standard deviation close to 0 and 1, respectively (Chang, Caruana,
and Goldenberg 2022).

##  Appendix D Experimental Setup

We conducted experiments to demonstrate the effectiveness of the HONAM. The
specifications of the experimental machine is Intel i7-8700 CPU, NVIDIA
GeForce RTX 3090 TI GPU, and 64 GB RAM.

### Dataset Split

For the California Housing, Insurance, House Prices, Bikeshare, FICO, Credit,
Support2, MIMIC, and COMPAS datasets, we randomly split
training/validation/test sets into 60%/20%/20% ratio with 5 random seeds. For
Year and Click datasets, we use training/validation/test sets provided by the
previous work (Popov, Morozov, and Babenko 2019).

### Hyperparameters

We use linear/logistic regression, CrossNet (Wang et al. 2017), XGBoost (Chen
and Guestrin 2016), MLP, EBM (Nori et al. 2019), NAM (Agarwal et al. 2020),
NodeGAM (Chang, Caruana, and Goldenberg 2022), NodeGA${}^{\text{2}}$M, and
HONAM as experimental methods. All models, except XGBoost and EBM, is
implemented using PyTorch. For MLP, NAM, and HONAM, the number of hidden
layers is set to 3, the numbers of hidden units are set to [32, 64, 32], and
LeakyReLU is used as the activation function. The number of units in CrossNet
is set to 32. All PyTorch models are trained for 1000 epochs with the learning
rate 0.001, and the batch size is used as the number of data samples divided
by 100 to reduce the training time. The best validation model is selectd for
the test. For NodeGAM and NodeGA${}^{\text{2}}$M, we use the default
hyperparameters provided by (Chang, Caruana, and Goldenberg 2022). We use the
open-source version of XGBoost and EBM. For XGBoost, the number of boosting
rounds is set to 1000 to ensure convergence, and eta is set to 0.3. For EBM,
the number of boosting round is set to 20000 to ensure convergence, inner bags
and outer bags are set to 8, and learning rate is set to 0.01 All models are
trained with 5 random seeds, and we report mean scores and standard
deviations.

### Evaluation Metrics

The mean squared error and root mean squared error are representative
evaluation metrics for regression. Still, because they are not scaled, it
cannot be easy to evaluate whether the model really works well. Therefore, we
used the R-squared score, which is a scaled evaluation metric. However,
R-squared score relies only on mean squared error. To evaluate the regression
models from a different perspective, we suggest R-absolute score, a novel
scaled regression metric that relies on the mean absolute error. The
R-absolute score is defined as follows:

$R\text{-}absolute=\frac{\sum_{i=1}^{N}\left|y_{i}-\hat{y}_{i}\right|}{\sum_{i=1}^{N}\left|y_{i}-\bar{y}\right|},$
(16)

where $N$ is the number of data samples. R-squared and R-absolute scores are
used as evaluation metrics for regression tasks, and AUROC and AUPRC are used
as evaluation metrics for classification tasks.

##  Appendix E Feature Interacion Ablation Study

3 | California Housing | Insurance | House Prices | Bikeshare | Year
---|---|---|---|---|---
| R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute | R-squared | R-absolute
2 CrossNet (t=2) | 0.717 ($\pm$0.001) | 0.525 ($\pm$0.001) | 0.860 ($\pm$0.005) | 0.688 ($\pm$0.020) | 0.885 ($\pm$0.011) | 0.683 ($\pm$0.003) | 0.498 ($\pm$0.008) | 0.316 ($\pm$0.009) | 0.310 ($\pm$0.001) | 0.211 ($\pm$0.001)
CrossNet (t=3) | 0.739 ($\pm$0.004) | 0.549 ($\pm$0.004) | 0.854 ($\pm$0.007) | 0.691 ($\pm$0.020) | 0.903 ($\pm$0.007) | 0.720 ($\pm$0.007) | 0.590 ($\pm$0.017) | 0.403 ($\pm$0.009) | 0.324 ($\pm$0.001) | 0.228 ($\pm$0.002)
CrossNet (t=4) | 0.741 ($\pm$0.007) | 0.562 ($\pm$0.006) | 0.855 ($\pm$0.010) | 0.690 ($\pm$0.007) | 0.903 ($\pm$0.007) | 0.718 ($\pm$0.008) | 0.645 ($\pm$0.013) | 0.453 ($\pm$0.011) | 0.327 ($\pm$0.002) | 0.236 ($\pm$0.002)
HONAM∗ (t=2) | 0.807 ($\pm$0.004) | 0.621 ($\pm$0.004) | 0.880 ($\pm$0.001) |  0.744 ($\pm$0.003) | 0.900 ($\pm$0.018) | 0.725 ($\pm$0.017) | 0.920 ($\pm$0.004) | 0.752 ($\pm$0.008) | 0.320 ($\pm$0.003) | 0.224 ($\pm$0.003)
HONAM∗ (t=3) |  0.810 ($\pm$0.008) |  0.629 ($\pm$0.009) |  0.882 ($\pm$0.001) | 0.743 ($\pm$0.005) | 0.905 ($\pm$0.020) | 0.730 ($\pm$0.022) | 0.945 ($\pm$0.005) | 0.806 ($\pm$0.010) | 0.329 ($\pm$0.004) | 0.235 ($\pm$0.006)
HONAM∗ (t=4) | 0.804 ($\pm$0.005) | 0.624 ($\pm$0.004) | 0.881 ($\pm$0.002) | 0.741 ($\pm$0.006) | 0.904 ($\pm$0.025) | 0.728 ($\pm$0.034) | 0.949 ($\pm$0.003) | 0.816 ($\pm$0.004) |  0.331 ($\pm$0.005) |  0.243 ($\pm$0.008)
HONAM (t=2) |  0.810 ($\pm$0.003) | 0.626 ($\pm$0.003) |  0.882 ($\pm$0.002) | 0.742 ($\pm$0.012) | 0.900 ($\pm$0.011) | 0.721 ($\pm$0.009) | 0.925 ($\pm$0.004) | 0.760 ($\pm$0.006) | 0.320 ($\pm$0.002) | 0.226 ($\pm$0.004)
HONAM (t=3) | 0.804 ($\pm$0.004) | 0.623 ($\pm$0.001) | 0.881 ($\pm$0.002) | 0.741 ($\pm$0.004) | 0.904 ($\pm$0.016) | 0.728 ($\pm$0.004) |  0.952 ($\pm$0.003) | 0.817 ($\pm$0.004) | 0.327 ($\pm$0.005) | 0.237 ($\pm$0.006)
HONAM (t=4) | 0.805 ($\pm$0.003) | 0.619 ($\pm$0.004) | 0.881 ($\pm$0.002) | 0.740 ($\pm$0.005) |  0.909 ($\pm$0.015) |  0.735 ($\pm$0.011) |  0.952 ($\pm$0.003) |  0.819 ($\pm$0.004) | 0.326 ($\pm$0.002) | 0.237 ($\pm$0.005)
3 |  |  |  |  |  |  |  |  |  |

Table 8: Interaction Ablation Study on Regression Tasks

3 | FICO | Credit | Support2 | MIMIC III | Click
---|---|---|---|---|---
| AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC | AUROC | AUPRC
CrossNet (t=2) | 0.770 ($\pm$0.002) | 0.752 ($\pm$0.004) | 0.964 ($\pm$0.024) | 0.810 ($\pm$0.038) | 0.803 ($\pm$0.009) | 0.612 ($\pm$0.013) | 0.787 ($\pm$0.005) | 0.343 ($\pm$0.022) | 0.619 ($\pm$0.001) | 0.625 ($\pm$0.001)
CrossNet (t=3) | 0.772 ($\pm$0.001) | 0.754 ($\pm$0.002) | 0.969 ($\pm$0.016) | 0.824 ($\pm$0.042) | 0.794 ($\pm$0.013) | 0.601 ($\pm$0.019) | 0.768 ($\pm$0.010) | 0.341 ($\pm$0.022) | 0.616 ($\pm$0.003) | 0.616 ($\pm$0.003)
CrossNet (t=4) | 0.772 ($\pm$0.002) | 0.755 ($\pm$0.003) | 0.965 ($\pm$0.024) | 0.838 ($\pm$0.033) | 0.798 ($\pm$0.007) | 0.605 ($\pm$0.010) | 0.775 ($\pm$0.009) | 0.342 ($\pm$0.027) | 0.617 ($\pm$0.003) | 0.619 ($\pm$0.003)
HONAM∗ (t=2) |  0.783 ($\pm$0.003) |  0.761 ($\pm$0.005) | 0.982 ($\pm$0.009) | 0.842 ($\pm$0.024) | 0.819 ($\pm$0.009) | 0.633 ($\pm$0.011) | 0.825 ($\pm$0.003) | 0.395 ($\pm$0.023) | 0.667 ($\pm$0.004) | 0.663 ($\pm$0.003)
HONAM∗ (t=3) | 0.781 ($\pm$0.002) | 0.759 ($\pm$0.003) |  0.987 ($\pm$0.006) |  0.853 ($\pm$0.031) | 0.819 ($\pm$0.008) | 0.632 ($\pm$0.014) | 0.824 ($\pm$0.002) | 0.388 ($\pm$0.023) | 0.667 ($\pm$0.003) | 0.656 ($\pm$0.003)
HONAM∗ (t=4) | 0.780 ($\pm$0.001) | 0.758 ($\pm$0.002) | 0.982 ($\pm$0.017) | 0.851 ($\pm$0.061) | 0.819 ($\pm$0.010) | 0.633 ($\pm$0.012) | 0.822 ($\pm$0.005) | 0.386 ($\pm$0.022) | 0.664 ($\pm$0.006) | 0.657 ($\pm$0.005)
HONAM (t=2) | 0.782 ($\pm$0.002) | 0.760 ($\pm$0.004) | 0.981 ($\pm$0.012) | 0.838 ($\pm$0.026) |  0.823 ($\pm$0.011) |  0.640 ($\pm$0.009) |  0.826 ($\pm$0.005) |  0.399 ($\pm$0.024) |  0.670 ($\pm$0.002) |  0.664 ($\pm$0.003)
HONAM (t=3) | 0.782 ($\pm$0.002) | 0.760 ($\pm$0.004) | 0.981 ($\pm$0.009) | 0.843 ($\pm$0.031) | 0.820 ($\pm$0.010) | 0.637 ($\pm$0.014) |  0.826 ($\pm$0.004) | 0.397 ($\pm$0.022) | 0.669 ($\pm$0.002) | 0.663 ($\pm$0.002)
HONAM (t=4) | 0.780 ($\pm$0.002) | 0.758 ($\pm$0.004) | 0.983 ($\pm$0.009) | 0.840 ($\pm$0.035) | 0.819 ($\pm$0.010) | 0.636 ($\pm$0.013) | 0.824 ($\pm$0.004) | 0.390 ($\pm$0.019) | 0.668 ($\pm$0.005) | 0.662 ($\pm$0.005)
3 |  |  |  |  |  |  |  |  |  |

Table 9: Interaction Ablation Study on Classification Tasks

We conduct feature interaction ablation study. Table 8 and Table 9 show the
predictive performances of CrossNet and HONAM. HONAM∗ indicates HONAM with
CrossNet, and HONAM indicates HONAM with the proposed feature interaction
method. We use 2${}^{\text{nd}}$-, 3${}^{\text{rd}}$-, and
4${}^{\text{th}}$-order CrossNet and HONAM for the ablation study. The
predictive performances are improved as the interaction order is increased.
However, significant performance improvements do not appear in some datasets.
Perhaps it is because they are vulnerable to overfitting, or there are no
feature interactions in the dataset. Nevertheless, we believe that high-order
interactions are important from interpretability perspective because they
allow the discovery of co-relationships between different features.

##  Appendix F Learning Sharp Jumps

### Synthetic Dataset Generation

The synthetic classification dataset was generated as follows: We randomly
sampled 100 points in [-1, 1] for the input $x$. We then randomly sampled $p$
in [0.1, 0.9] and generated 100 labels (0 or 1) from a binomial distribution
with probability $p$ for each $x$. The synthetic classification dataset
contained 10,000 samples. The shape function of the training model should fit
the log odds of the probability $p$. Additionally, we generated a synthetic
regression dataset for an additional experiment. The synthetic regression
dataset was generated as follows: We randomly sampled 100 points in [-1, 1]
for input $x$, and the corresponding target $y$ was also sampled in [-1, 1].
The synthetic regression dataset contained 100 samples.

### Learned Shape Plots on Synthetic Dataset

(a) ExU

(b) ExpDive

Figure 8: Learned shape functions of different units with single feature
network on synthetic regression task

(a) ExU

(b) ExpDive

Figure 9: Learned shape functions of different units with single feature
network on synthetic classification task

(a) ExU

(b) ExpDive

Figure 10: Learned shape functions of different units with multiple feature
networks on synthetic regression task

(a) ExU

(b) ExpDive

Figure 11: Learned shape functions of different units with multiple feature
networks on synthetic classification task

Figure 8 and Figure 9 show the shape plots learned by ExU, and ExpDive units
on synthetic regression and classification datasets, respectively. Consistent
with Theorem 2, the ExU unit cannot learn well in the range of negative $x$ on
both datasets. Our ExpDive unit learned the sharpest function on both
datasets. In addition, we found that the joint learning models can learn
extremely sharp shape functions. Figure 10 and 11 show the shape plots learned
by the joint learning model with 100 ExU, and ExpDive units, respectively.
Since the joint learning models have high capacity and sub-networks cooperate
to predict a single data point, they can learn extremely sharp shape
functions. Therefore, the ExU unit works well in the negative range of $x$.

##  Appendix G Interpreting HONAM

### Extracting Feature Shape Plots

We can acquire direct feature contribution values via forward propagation of
HONAM. For example, we refer the representation vectors of each input feature
${x_{1},x_{2},…,x_{m}}$ computed by the feature networks as
${r_{1},r_{2},…,r_{m}}$, where $m$ is the number of features. To obtain 1st-
order interpretation, the contribution value of feature $i$ is computed by
$r_{i}\times w_{0:k}$, where $w$ is the output matrix and $k$ is the size of
the representation vector. Similarly, the contribution value of 2nd-order
feature $(x_{i}*x_{j})$ can be computed as $(r_{i}*r_{j})\times w_{k:2*k}$.
The purification, as in NodeGAM, is not needed because, 1) in our HONAM, a
feature has only one corresponding representation vector, 2) this
representation vector is directly connected to the output layer, and 3) our
feature interaction method does not learn multi-terms polynomials. For
example, our interaction method only learns interactions between different
features such as $x_{1}\times x_{2}$ or $x_{1}\times x_{3}\times x_{7}$, it
does not learn multi-terms polynomials such as $x_{1}x_{2}+x_{2}$ or
$x_{2}x_{3}+x_{2}$.

### Visualization of Feature Contributions

Figure 12 and 13 show the shape plots of the 1${}^{\text{st}}$-order and
2${}^{\text{nd}}$-order features on the CA Housing dataset. Figure 14 and 15
show the shape plots of the 1${}^{\text{st}}$-order and
2${}^{\text{nd}}$-order features on the FICO Housing dataset.

We can get some insights from the figures. First, HONAM can learn soft (or
smooth) nonlinear functions because it uses neural networks. Because the
conventional GAM can only learn linear functions, its expressive power is
limited. In addition, because tree-based models, such as boosting trees, can
only learn discretized functions, they require a massive number of trees to
learn smooth-looking functions. Therefore, it is difficult to extend various
settings, such as multi-label or multi-class learning (Agarwal et al. 2020),
and their explanations are difficult to understand. Second, HONAM can learn
sharp jumps. In the California Housing dataset, HONAM learned a few sharp
jumps. We believe that our ExpDive unit is more effective in a field where
sharp jumps frequently occur, such as in stock price predictions. Third, we
can view the relationships between features via HONAM. In other words, HONAM
provides more detailed interpretations than NAM and other interpretable
methods.

##  Appendix H Discussion

In this section, we would like to discuss the limitations of interpretable
machine learning, NAM-family methods, and HONAM.

Interpretable machine learning methods allow us to understand their decision-
making processes. However, it is impossible to determine why certain features
are important and why certain features have positive or negative influences.
For example, in Figure 4(a), ExternalRistEstimate has the greatest influence.
In other words, the model determined that this feature was the most important;
however, why? Human assistance is required to understand this. Therefore,
interpretable machine learning still requires human effort, and we take care
when interpreting the results of a model.

The NAM family of methods requires different feature networks corresponding to
each feature. This strategy creates an interpretable model but is time-
consuming. Fast training and inference are important in real-world systems.
Therefore, we need an interpretable method that is faster than the NAM.

HONAM uses all feature combinations to model interactions. However, modeling
all interactions slows down the execution time and can reduce the predictive
performance because some interactions may be noisy (Cheng, Shen, and Huang
2020). If HONAM can adaptively learn feature interactions, it may have higher
predictive performance, interpretability, and faster execution time.

Figure 12: Learned shape plots of 1${}^{\text{st}}$-order features on the CA
Housing dataset

Figure 13: Learned shape plots of 2${}^{\text{nd}}$-order features on the CA
Housing dataset

Figure 14: Learned shape plots of 1${}^{\text{st}}$-order features on the FICO
dataset

Figure 15: Learned shape plots of 2${}^{\text{nd}}$-order features on the FICO
dataset

Algorithm 1 Forward propagation of Feature Network

Input:

input features $x$

Parameter:

number of layers $l$

trainable weights $w_{1},w_{2},\ldots,w_{l}$

trainable biases $b_{1},b_{2},\ldots,b_{l}$

1:for $i\leftarrow 1$ to $l$ do

2: $x\leftarrow xw_{i}+b_{i}$ $\triangleright$ linear transformation of $x$

3: $x\leftarrow\sigma\left(x\right)$ $\triangleright$ activation function

4:end for

5:return $x$

Algorithm 2 Forward propagation of HONAM

Input:

input features $x$

Parameter:

number of features $m$

order of interaction $t$

feature networks corresponding to each feature $f_{1},f_{2},...f_{m}$

trainable weight $w$

trainable bias $b$

1:for $i\leftarrow 1$ to $m$ do

2: $repr_{i}\leftarrow f\left(x_{i}\right)$ $\triangleright$ run feature
network, see Algorithm 1

3:end for

4:$fi_{0},pfi_{0}\leftarrow 1,1$ $\triangleright$ init interactions

5:$fi_{1},pfi_{1}\leftarrow\sum_{i=1}^{m}repr_{i},\sum_{i=1}^{m}repr_{i}$
$\triangleright$ sum of 1${}^{\text{st}}$-order interactions

6:for $i\leftarrow 2$ to $t$ do $\triangleright$ high-order interaction
modeling

7: $sign\leftarrow 1$

8: $pfi_{i}\leftarrow\sum_{k=1}^{m}repr^{i}_{k}$

9: $fi_{i}\leftarrow 0$

10: for $j\leftarrow 1$ to $i$ do

11: $fi_{i}\leftarrow fi_{i}+sign\times pfi_{j}\times fi_{i-j}$

12: $sign=-sign$

13: end for

14: $fi_{i}\leftarrow fi_{i}/i$

15:end for

16:$pred\leftarrow fi_{1:m}w+b$ $\triangleright$ prediction

17:return $pred$