# Predicting the cascading dynamics in complex networks via the bimodal
failure size distribution

Chongxin Zhong  School of Systems Science, Beijing Normal University, Beijing,
P.R. China  Yanmeng Xing  School of Systems Science, Beijing Normal
University, Beijing, P.R. China  Ying Fan  School of Systems Science, Beijing
Normal University, Beijing, P.R. China  An Zeng  School of Systems Science,
Beijing Normal University, Beijing, P.R. China<EMAIL_ADDRESS>
###### Abstract

Cascading failure as a systematic risk occurs in a wide range of real-world
networks. Cascade size distribution is a basic and crucial characteristic of
systemic cascade behaviors. Recent research works have revealed that the
distribution of cascade sizes is a bimodal form indicating the existence of
either very small cascades or large ones. In this paper, we aim to understand
the properties and formation of such bimodal distribution of cascade sizes in
complex networks, and further predict the final cascade size. We first find
that the bimodal distribution of cascade sizes is ubiquitous in both synthetic
and real networks. Moreover, the large cascade sizes distributed in the right
peak of bimodal distribution are resulted from either the failure of nodes
with high load at the first step of the cascade or multiple rounds of cascades
triggered by the initial failure. Accordingly, we propose a hybrid load metric
(HLM), which combines the load of the initial broken node and the load of
failed nodes triggered by the initial failure, to predict the final size of
cascading failures. Finally, we validate the effectiveness of HLM by computing
the accuracy of identifying the cascades belonging to the right and left peaks
of the bimodal distribution. The results show that HLM is a better predictor
than commonly used network centrality metrics in both synthetic and real-world
networks.

Keywords:  cascading failure, bimodal distribution, complex network

##  1 Introduction

Cascading failure is a common phenomenon in many complex systems [1, 2, 3].
The vast majority of large-scale cascading failures can be traced back to a
single-element failure. With the catastrophic effects[4, 5], a failure (due to
a random breakdown or intentional attack) of systems can trigger further
failures which may lead to the whole system being dysfunctional and
collapsing. The cascades not only take place in infrastructure like electrical
power grids and transportation networks [6, 7], but also present in social and
economic systems [2, 8]. Applying network science to uncover complicated
dynamic behaviors emerging from these real complex systems is an important
research direction [9]. Stacks of studies have already discussed various
important aspects of cascading dynamics in complex networks, including
modeling the process of failure propagation[4, 10, 11, 12, 13], vulnerability
and robustness of networks [14, 15, 16, 17], mitigation and recovery
strategies for cascading failures [18, 19, 20], and the nature of phase
transition induced by attacks[21, 22].

One of the most paradigmatic models depicting propagation mechanism of
cascading failure is the Motter-Lai model[4], which shows the non-local faults
diffusion driven by global load redistribution. Non-local cascades have been
observed in real systems like blackouts or air-traffic disruptions, all of
which can cause global damage to systems[23]. In particular, cascade size
distribution is a basic and crucial characteristic of systemic cascade
behaviors[24]. From the empirical data, it has been shown that the blackout
size distribution approximately follows a power law indicating that the power
grid is in a state of self-organized criticality[25]. Many cascade models have
been established to reproduce power-law distributions of cascade sizes[10, 6,
26]. However, recent studies have suggested that the cascade sizes in complex
networks follow a bimodal distribution with either very small cascades or
large ones, which is a typical appearance that usually shows up in the abrupt
phase transition[27, 28, 29, 22, 30]. The abrupt transition is so sudden and
devastating that it is harder to be predicted than a continuous transition.
Recently, some researchers have put their attention on the bimodal
distribution of cascade sizes[31, 30].

Previous research works have revealed that both power-law and bimodal
distributions of cascading failures exist in various cascade models. The
presence of power-law and bimodal distributions have a dependence on network
topology, model parameters, and rules of cascading dynamics[31, 30]. However,
current research works focus more on power-law distribution and few research
works on bimodal distribution of cascades. These works only show the bimodal
distribution phenomenon observed in various cascade models, however, the
properties and formation of the bimodal distribution still remain to be
explored. Understanding the formation of the bimodal distribution of cascading
sizes would help us not only to predict the damage caused by the cascading
dynamics, but also to design methods to mitigate catastrophic cascades.
Therefore, it is necessary to get an exhaustive comprehension on the rules
behind bimodal cascading sizes and further aid in prediction of the final
cascade size.

In this paper, we investigate in detail the formation of the bimodal
distribution of cascades, aiming to distinguish that whether a node belongs to
the category in which its breakdown completely disrupts the whole network or
leaves it almost intact. To this end, we study the characteristics of nodes
among the networks triggering different peaks of the bimodal distribution of
cascade sizes after a single initial failure. We find that the final cascade
size has close ties with the maximal load of nodes that failed at the first
step of the cascade as well as the number of sequential rounds of cascades
triggered by the initial failure. Based on this, we propose a hybrid load
metric (HLM), in which the load of the initial failed node and the maximal
load of First Failures (the failed nodes at the first step of a cascade) are
considered, to predict the final cascade size. We then employ AUC, defined by
the area under the receiver operating characteristic curve, to measure the
accuracy of HLM in identifying a cascade belonging to the left or right peak
of the bimodal distribution. By validating on synthetic and real-world
networks, we show that HLM outperforms the commonly used network centrality
metric, betweenness, in predicting final cascade sizes.

##  2 Cascading failure model

Real networks sometimes undergo undesirable attacks or errors of nodes
(edges), which can trigger further failures in other components. There is an
established fact that cascading failures can range from only sporadic
components failing to catastrophic collapse, as shown in Fig. 1. In IEEE 118
Bus[32], the initial failure of node $i$ and node $j$ with similar load
trigger completely different cascade sizes. The former causes nearly 40% of
nodes to fail, while the latter only leads to a few nodes being broken. Recent
studies have indicated that the distribution of avalanche sizes is bimodal in
which some nodes only trigger a limited scope of cascading failures, while
others would cause fatal global cascades. In this paper, we are mainly
concerned about that how the bimodal distribution of cascades forms and how to
distinguish whether a node is capable to bring a widespread network crash that
results in a cascade size belonging to the right peak of the bimodal
distribution.

Figure 1: Illustration of differences in cascade sizes triggered by two
individual nodes in a power grid. The network is IEEE 118 Bus
$\left(N=118,\langle k\rangle=1.58\right)$, in which nodes represent buses and
links represent transmission lines. The initial attacked node is marked in
red. The red nodes in (a) and (b) are with a similar load. The yellow nodes
are subsequent failures triggered by removing the red node and the links
connecting to it. The Motter-Lai model is employed to simulate the cascading
dynamics here.

First of all, we briefly introduce the model of cascading failures that we
employed. We consider a simple overload model, which is proposed by Motter and
Lai[4]. In the Motter-Lai model, node $i$ has a certain load $L_{i}$ given by
betweenness, which is the number of shortest paths through this node.
$\left(1+\alpha\right)$ times the initial load of node $i$ ($L_{i}$) is
defined as its capacity $C_{i}$ , which is the maximum load that node $i$ can
carry.

$\displaystyle C_{i}=\left(1+\alpha\right)L_{i},\qquad i=1,2,3...N,$ (1)

where $\alpha$ is a tunable parameter, called the tolerant parameter. $N$ is
the number of nodes in the network. In this paper, we fix $\alpha=0.1$, a
relatively small value, so that the cascade sizes of nodes are not so large
and the cascading results are largely affected by the properties of the
initially failed node. In a network, once a node is attacked, all edges
connected to this node are also deleted. Therefore, some shortest paths may be
diverted, resulting in a change of the load of remaining nodes. If the load of
remaining nodes exceed their capacities, these nodes will also be removed.
Repeat above process until no node fails in the network, as the end of a
cascade. The number of times that this process repeats is donated as $T$,
which indicates the number of rounds of cascades triggered by the initial
failure. To simplify the description, the first round of cascades triggered by
the initial failure is called the first step of a cascade, and so forth. This
model is paradigmatic of non-local propagation of failures and widely used in
research works to study various faults propagation.

##  3 Results

We first investigate the kinds of networks where the bimodal distribution of
cascade sizes can emerge. In particular, we consider a vital property of
networks, the average degree $\langle k\rangle$, which tremendously affects
the ability of networks to withstand failures. As for information and disease
transmission, it is easier for networks with a higher average degree to have
wider propagation.[24]. However, cascading failure dynamics is different. For
example, a network with a higher $\langle k\rangle$ , in which more nodes can
share the redistributed load of the failed nodes, is less inclined to trigger
a large cascade. Therefore, we first study how the average degree of the
network impacts the distribution of cascade sizes. The relative size of a
cascade is quantified by the ratio of failed nodes after the initial attack in
the network. We stimulate the Motter-Lai model on Erdős–Rényi (ER) random
networks, Small-world (SW) networks and Barabasi-Albert (BA) networks,
respectively. In each realization, we attack one node (delete this node and
all edges connecting to it) and stimulate the Motter-Lai model to calculate
final cascade sizes in the network with a fixed $\langle k\rangle$. The size
distribution of cascading failures is recorded, until we have attacked all
nodes of networks separately. The results are shown in Fig. 2. The redder the
color, the higher frequency of such cascade sizes within the corresponding
interval in networks with corresponding $\langle k\rangle$. We can clearly see
that there are bimodal distributions of cascades in all synthetic networks,
shown by two redder regions longitudinally at a fixed $\langle k\rangle$. For
the ER networks and the BA networks, the bimodal distribution of cascades is
more likely to occur in networks with a lower average degree. In networks with
too small or too large $\langle k\rangle$, cascade size distributions are
unimodal and cascade sizes are mainly at a low level close to 0. This
phenomenon is resulted from that a network with too small $\langle k\rangle$
is prone to split the network into some separate connected components that can
prevent the propagation of failures and those with too large $\langle
k\rangle$ have more links to share the redistributed load of failed nodes[30].
The SW networks, however, exhibit different characteristics that the bimodal
distribution always exists, even though the average degree $\langle k\rangle$
is large. This indicates that the small-world property might bring more
cascading risks for networks.

Figure 2: The distribution of cascade sizes in synthetic networks with
different average degree $\mathbf{\langle k\rangle}$. (a), (b) and (c) show
the results in ER networks, SW networks and BA networks respectively. The
cascade is triggered by a single initial broken node in the network. We get
the cascade size of each node and then take binning frequency statistics,
which is shown by the heat map. The redder the color is, the higher frequency
of such cascade sizes is within the corresponding interval in networks with a
specific average degree. Two redder regions longitudinally at a fixed $\langle
k\rangle$ indicate that there is a bimodal distribution of cascade sizes in
networks with the corresponding $\langle k\rangle$. The Motter-Lai model is
employed to simulate the cascading dynamics. The results are computed from the
10 independent networks with a fixed $\langle k\rangle$ and $N=500$. Figure 3:
The spreading distance $\langle d\rangle$ of failures under different cascade
sizes. The networks are: (a) ER networks $\left(N=500,\langle
k\rangle=2\right)$, (b) SW networks $\left(N=500,\langle k\rangle=4\right)$,
(c) BA networks $\left(N=500,\langle k\rangle=4\right)$. The spreading
distance $\langle d\rangle$ shown by yellow scatters is defined as the average
topological distance between the initial attacked node and subsequent failed
nodes. Here, subsequent failed nodes (failures) refer to all failures
triggered by the initial single failure. And the $\langle d\rangle$-random
shown by orange scatters is the result of randomized $\langle d\rangle$ that
is obtained by recalculating the mean distance between randomly selected
failed nodes (preserving the failure size) and the initial broken node. The
insert in (c) is the distribution of cascade sizes in the double logarithm
coordinate. The solid green line indicates the inflection point of $\langle
d\rangle$. The above results are computed in 100 independent networks with the
same $\langle k\rangle$. Figure 4: The relationship between First Failures and
final cascade failures. Networks are : (a, d) ER networks $\left(N=500,\langle
k\rangle=2\right)$, (b, e) SW networks $\left(N=500,\langle
k\rangle=4\right)$, (c, f) BA networks $\left(N=500,\langle
k\rangle=4\right)$. First Failures is defined by the failures (failed nodes)
at the first step of a cascade and those in the second and later steps is
called Second Failures. (a), (b) and (c) show the distribution of First
Failures sizes and Second failures sizes in corresponding synthetic networks.
(d), (e) and (f) show the relationship between First Failures sizes ($S_{f}$)
and final cascading failures sizes ($S_{all}$) by the intensity scatter. The
color of each point represents the density of nodes within a rectangular area
centered on it. The Kendall’s tau rank correlation coefficient $\tau$ is
applied to estimate how the First Failures size is correlated to the final
cascade size. The above results are computed in 100 independent networks with
the same $\langle k\rangle$.

The bimodal distribution of cascades has been observed in three different
paradigmatic synthetic networks. It is natural to consider what are the
differences between the cascades in different peaks of the bimodal
distribution. To answer this, we further analyze the relationship between the
spreading distance $\langle d\rangle$ of failures and cascade sizes to
understand the impact of the bimodal distribution on the propagation range.
The spreading distance $\langle d\rangle$ is defined as the average
topological distance between the initial attacked node and subsequent
failures, which can imply the nonlocality of failure propagation[23]. A low
$\langle d\rangle$ means that the cascading failures are near the initial
breakdown (locality), otherwise are far from the initial one (nonlocality).
The bimodal distribution of cascades, in which the left peak separates from
the right peak by a gap where cascades rarely occur at these sizes, is shown
by a bar chart in different synthetic networks in Fig. 3. The horizontal axis
represents cascade size and the blue bars represent the corresponding
probability, which is the proportion of a fixed cascade size among cascade
sizes of 50000 nodes from 100 independent networks with the same average
degree $\langle k\rangle$. We can see that the bimodality is relatively weak
and cascade sizes are also lower in BA networks, which reflects the robustness
of BA networks against cascading failures. For spreading distance $\langle
d\rangle$, shown by yellow scatters, it grows as the cascade size increases
and then roughly keeps the same with the randomized $\langle d\rangle$
($\langle d\rangle$-random), which is obtained by recalculating the mean
distance between randomly selected failed nodes (preserving the failure size)
and the initial broken node. As the reference of real $\langle d\rangle$, the
randomized $\langle d\rangle$ represents a random counterpart where the
failures are cascaded randomly and globally. Obviously, there is an inflection
point of spreading distance $\langle d\rangle$ around the gap between two
peaks in these networks. To the left of the inflection point, $\langle
d\rangle$ increases as the cascade size grows and is significantly lower than
the randomized $\langle d\rangle$. However, to the right of the inflection
point, the average distance $\langle d\rangle$ is close to the randomized
$\langle d\rangle$, which indicates that such cascade sizes in the right peak
could cause the same range of damage as the random counterpart. As the results
shown above, the cascade sizes in the right peak are nonlocal during the
propagation of failures, which can destroy the network globally and randomly.
Note that the average distance is slightly lower than the randomized one on
the right of the inflection point in ER networks. This implies that the
failure propagation in homogeneous networks tends to be more local than those
in heterogeneous networks, when compared to a random counterpart.

Figure 5: The dependence of the cascade size on the maximal load of First
Failures (MLF) and cascading rounds $T$ triggered by the initial failure in
different synthetic networks. Networks are : (a) ER networks
$\left(N=500,\langle k\rangle=2\right)$, (b) SW networks $\left(N=500,\langle
k\rangle=4\right)$, (c) BA networks $\left(N=500,\langle k\rangle=4\right)$.
(d) shows the characteristics of MLF and $T$ in these networks. The blue bar
is the ratio of nodes in the right peak with higher MLF than the median MLF of
the whole data ($\rm{MLF_{r}>MLF_{1/2}}$). The orange bar has a similar
meaning that the ratio of nodes in the right peak with higher $T$ than the
median $T$ of the whole data. The green bar shows the ratio of nodes with
higher MLF or $T$ in the right peak than the corresponding median of the whole
data. The above results are computed in 100 independent networks with the same
$\langle k\rangle$.

The spreading distance $\langle d\rangle$ of failures in the two peaks shows
nonlocality characteristics of cascades. The nonlocality of the cascade
suggests that the damage is distributed globally, which is hard to predict and
recover. Indeed, figuring out the formation of bimodal distribution is the key
step in predicting and preventing large-scale cascades. Next, we aim to
investigate how the bimodal distribution forms. Firstly, we consider
intuitively whether the failure size at the first step of a cascade determines
the final cascade size. We donate the failures (failed nodes) at the first
step of a cascade as First Failures and those at the second and later steps
are called Second Failures here. In Figs. 4 (a), (b) and (c), we show the
distributions of First Failures sizes and Second Failures sizes in ER
networks, SW networks and BA networks respectively. Interestingly, the
distribution of First Failures sizes is unimodal and narrow. The result is
totally different from the original bimodal distribution, which implies that
First Failures cannot represent the final cascading failures. However, the
distribution of Second Failures sizes shows weak bimodality, consistent with
the original bimodal distribution. In addition, we quantify the dependence of
the final cascade sizes on First Failures by the Kendall’s tau rank
correlation coefficient[33] $\tau$. In Figs. 4(a), (b), we do not observe
obvious linear relation between final cascade sizes and First Failures sizes,
and the Kendall’s tau rank correlation coefficients are also small. In Fig.
4(c), the correlation coefficient is slightly higher. This may be because of
high heterogeneity in BA networks, considering the initial failed nodes with
high heterogeneity would play decisive roles in their cascading outcomes. To
sum up, the final bimodal distribution of cascade sizes is mainly caused by
the Second Failures and has no significant correlation with First Failure in
ER networks and SW networks. This is significantly different from the
propagation of disease and information, where the wider spread probably
emerges from more infected agents at the first step of propagation[34].

Figure 6: The AUC value of HLM under different parameter $\lambda$. The networks are: (a) ER networks $\left(N=500,\langle k\rangle=2\right)$, (b) SW networks $\left(N=500,\langle k\rangle=4\right)$, (c) BA networks $\left(N=500,\langle k\rangle=4\right)$. The results are obtained by averaging 10 independent realizations. Table 1: Structural properties and AUC value of the different metrics in synthetic and real networks. Structural properties include network size $(N)$ and average degree ($\langle k\rangle$). Metrics include the load of the initial attacked node $(L)$, size of First Failures ($S_{1}$), the spreading distance ($d_{1}$) at the first step that is defined by the average distance between the initial attacked node and those of First Failures, the maximal load of First Failures (MLF) and HLM with an optimal $\lambda$ achieving the largest AUC value ($\rm{HLM^{*}}$). Networks | $\mathbf{N}$ | $\mathbf{\langle k\rangle}$ | $\mathbf{\langle d_{1}\rangle}$ | $\mathbf{S_{f}}$ | $\mathbf{L}$ | MLF | $\mathbf{\rm}{HLM^{*}}$
---|---|---|---|---|---|---|---
BA | 500 | 4 | 0.8303 | 0.8375 | 0.8483 | 0.8924 | 0.9022
ER | 500 | 2 | 0.7891 | 0.7833 | 0.8157 | 0.8405 | 0.8546
SW | 500 | 4 | 0.5747 | 0.5597 | 0.7646 | 0.8062 | 0.8339
Bus 1494 | 1494 | 2.89 | 0.7763 | 0.7967 | 0.8145 | 0.8453 | 0.8459
Euro-road | 1039 | 2.51 | 0.5000 | 0.7185 | 0.6816 | 0.7949 | 0.7949
Bus 494 | 494 | 2.37 | 0.7800 | 0.7752 | 0.8134 | 0.8750 | 0.8785
Econ-poli | 2343 | 2.28 | 0.7772 | 0.8053 | 0.8455 | 0.9038 | 0.9039
Bio-yeast | 1458 | 2.67 | 0.6798 | 0.8140 | 0.8263 | 0.8263 | 0.8821
Figure 7: The distribution of First Failures sizes and Second Failures sizes
in real-world networks. The networks are:(a) Bus 494 $\left(N=494,\langle
k\rangle=2.37\right)$, (b) Euro-road $\left(N=1039,\langle
k\rangle=2.51\right)$, (c) Econ-poli $\left(N=2343,\langle
k\rangle=2.28\right)$, (d) Bio-yeast $\left(N=1458,\langle
k\rangle=2.67\right)$.

Based above analysis, the size of First Failures does not determine the final
cascading size to a large extent, which implies that it cannot be used to
predict the final cascade size. Now, we turn our attention to the
characteristics of First Failures to see whether they influence the final
cascade size. The load of nodes plays an essential role, which not only
determines the redistributed flow, but also the capacity of nodes in the
cascading dynamics. Therefore, we further explore the dependence of the
cascade size on the maximal load of First Failures (MLF) in different
artificial networks. The number of rounds of cascades $T$ triggered by the
initial failure as a vital variable of cascading dynamics is also considered
here. The results are shown in Fig. 5. Each cascade is triggered by a single
node. MLF and $T$ are recorded after the cascade ends. Therefore, every node
will be with a cascade size, an MLF value and a $T$ value. Beforehand, we have
observed that the distribution of cascade sizes is bimodal in these networks
in Fig. 3. The boundary between the left peak and right peak is defined by the
cascade sizes with the lowest probability in the gap of two peaks here. In
Figs. 5 (a), (b) and (c), it is obvious that the large cascade sizes
distributed in the right peak are associated with higher MLF and $T$ than
those in the left peak. To characterize this phenomenon more precisely, we
further measure the values of MLF and $T$ of the cascades in the right peak
compared with the whole data. We first calculate the median MLF value and $T$
value of all nodes in networks ($\rm{MLF_{1/2}}$, $T_{1/2}$). Then, we
calculate the ratio of nodes in the right peak with higher MLF values than
$\rm{MLF_{1/2}}$. The higher such ratio is, the more nodes would be with
higher MLF values in the right peak. Similarly, the ratio of nodes in the
right peak with higher $T$ values than $T_{1/2}$ is recorded. Moreover, the
union of them that the percentage of nodes in the right peak with higher MLF
values or $T$ values than their corresponding median of the whole data is also
considered. The results are shown in Fig. 5 (d). We can see that the ratios
are over 0.5 and the union of them is more than 0.9 in these networks, which
indicates that most large cascades in the right peak tend to be with either
higher MLF values or higher $T$ values. A higher MLF value means that the
nodes with high load failed at the first step of the cascade and a higher $T$
means that more rounds of cascades triggered by the initial failure. In other
words, most large cascades in the right peak is resulted from either explosive
breakdown in only a few steps or accumulation of multiple rounds of cascades
triggered by the initial failure.

The higher MLF and $T$ in the right peak where the larger cascade sizes are
distributed, can provide us a new view of predicting the sizes of cascading
failures. The MLF value only contains the information from the first step of
the cascade, although $T$ is derived from the entire process of cascading
dynamics. In the previous studies, the load of the initial attacked node is
always considered to predict the final cascade size. However, there are some
scenarios ignored, in which the final damage is not strongly dependent on the
initial failure but closely related to the detailed cascade process. For
instance, a node with a low load failed might trigger the failure of a core
node, which would cause a large cascade. Accordingly, we propose a hybrid load
metric (HLM), combining the load of initial failure and MLF with a tunable
parameter, to predict the final cascade size after an initial attack. HLM
value of node $i$ is defined as

$\displaystyle\centering{\rm{HLM}}_{i}=\lambda{L}_{i}+\left(1-\lambda\right){\rm{MLF}_{i}},\@add@centering$
(2)

where $L_{i}$ is the load of the initial attacked node $i$. $\rm{{MLF}_{i}}$
is the maximal load of First Failures triggered by the initial attacked node
$i$. $\lambda$ is a tunable parameter controlling their weight. Actually, when
$\lambda=0$, HLM is equivalent to MLF. When $\lambda=1$, HLM returns to the
commonly used centrality index, namely the betweenness (load) of the initial
failed node. To validate the effectiveness of HLM, we employ Area Under Curve
(AUC)[35] to evaluate its effectiveness for predicting the final cascade size.
We randomly choose a node in the left and right peaks respectively and then
compare their HLM values at each time. In $n$ independent comparisons, if
there are $n^{\prime}$ times the node in the right peak having a higher HLM
and $n^{\prime\prime}$ times they have the same HLM value, the AUC is
expressed as

$\displaystyle\centering{\rm{AUC}}=\frac{n^{\prime}+0.5n^{\prime\prime}}{n}.\@add@centering$
(3)

The AUC value is about 0.5 for the random prediction. The AUC value can be
interpreted as the probability that a randomly chosen node in the right peak
is given a higher HLM value than that in the left peak. A higher AUC value
implies that the metric has a better performance on prediction or
classification. We first investigate how the parameter $\lambda$ affects the
AUC value in synthetic networks, which is shown in Fig. 6. We can see that
there is an optimal $\lambda$ achieving the largest AUC in three synthetic
networks and it is apparently higher than the AUC value when $\lambda=0$ or 1.
In addition, AUC with $\lambda=0$ is higher than that with $\lambda=1$, which
means that MLF has a better performance than $L$ in predicting the final
cascade size. Hence, the MLF value can better symbolize the final cascade size
compared to the load of initial broken node $L$, which is commonly used. The
detailed results are reported in Table. 1. For comparison, we also test other
general metrics including the load of the initial attacked node $L$, the size
of First Failures $S_{f}$, and the spreading distance $\langle d_{1}\rangle$
at the first step of cascades which is the average distance between the
initial attacked node and those of First Failures. We can see that MLF and HLM
apparently outperform other metrics, and HLM is the best one. In fact, HLM
combining the load of initial failure and information of subsequent failed
nodes, makes up for the deficiency of only focusing on the load of the initial
failed node by taking the cascading dynamics at the early stage of cascades
into consideration.

To further support the effectiveness of our method, we apply it to real
networks including power grids (Bus 1494, Bus 494)[36], a road network (Euro-
road)[37], an economic network (Econ-poli)[36] and a biological network (Bio-
yeast) [38]. Their structural properties are shown in Table. 1. Likewise,
there is a bimodal distribution of cascade sizes in each of these real
networks, part of which are shown in Fig. 7. First failures and Second
Failures also obey the same features as in synthetic networks that the
distribution of First Failures sizes is unimodal and the distribution of
Second Failures sizes shows weak bimodality. The AUC of the metrics in real
networks is displayed in Table. 1. MLF and HLM also outperform other metrics,
which is consistent with those observed in synthetic networks.

##  4 Discussion

The cascade size is a big concern when researchers discuss cascading failures.
Numerous research works have focused on the power-law distribution of cascade
sizes, which implies that the system is driven by a self-organized state.
However, recent studies suggest that the distribution of cascade sizes is
bimodal in which the cascade size triggered by one node is either very small
or very large. Indeed, both power-law and bimodal distributions of cascading
failures exist in complex networks, which depend on network topology, model
parameters and rules of cascading dynamics. Previous studies have intensively
discussed the power-law cascade size distribution. However, the bimodal
distribution is rarely involved. Therefore, it is necessary to study the rules
and characteristics behind the bimodal distribution, which can help to
understand deeply cascading dynamics and aid prediction of the final cascade
size.

In this article, we focus on the formation and prediction of bimodal
distribution of cascade sizes. What we do first is to investigate the features
of nodes in different peaks of the bimodal distribution. The results show that
the large cascade sizes distributed in the right peak are resulted from either
the failure of higher-load nodes at the first step or triggering multiple
rounds of cascades after the initial failure. In other words, large-scale
cascading failures are the results of either an explosive breakdown in only a
few steps or cumulative multiple rounds of cascades triggered by an initial
attack. In this context, we propose MLF (the maximal load of First Failures)
and a hybrid load metric (HLM) combining the load of the initial attacked node
and MLF with a tunable parameter $\lambda$. Then, employ them to distinguish
whether the breakdown of a node completely disrupts the whole network or
leaves it almost intact, which respectively correspond two peaks of the
bimodal distribution. By comparing the AUC value, both MLF and HLM show better
performance compared to common indicators, such as the load of the initial
attacked node and the size of First Failures. Moreover, real-world networks
also show similar results.

The bimodal distribution of cascade sizes is an interesting phenomenon that
deserves further effort to investigate the mechanism behind it. Our work
reveals that the statistical properties of bimodal cascade size distribution
and propose an effective metric HLM to estimate the size of cascading failures
triggered by a single node. Despite its effectiveness, there are several
limitations of it. For example, HLM can achieve only a rough distinction of
cascades between two peaks of bimodal distribution, not a highly accurate
prediction. More accurate prediction methods are asked for further
exploration. Moreover, the information of the first step in cascade progress
is also involved in HLM, which means we cannot predict by HLM before the
cascades occur. In fact, the prediction of the failure spreading is rarely
studied, because cascades are the results of global interactions among nodes
which would make the prediction harder. Our object is to predict the final
cascade at the early stage of cascades. Finally, our research only considers a
simple case where the cascades are triggered by a single node, and does not
involve the cascades due to massive initial failures, which is also an
important research direction of cascading dynamics asked for further
investigation.

## References

  * [1] Albert, R., Jeong, H. & Barabási, A.-L.  Error and attack tolerance of complex networks.  _nature_ 406, 378–382 (2000).
  * [2] Haldane, A. G. & May, R. M.  Systemic risk in banking ecosystems.  _Nature_ 469, 351–355 (2011).
  * [3] Fairley, P.  The unruly power grid.  _IEEE Spectrum_ 41, 22–27 (2004).
  * [4] Motter, A. E. & Lai, Y.-C.  Cascade-based attacks on complex networks.  _Physical Review E_ 66, 065102 (2002).
  * [5] Pourbeik, P., Kundur, P. S. & Taylor, C. W.  The anatomy of a power grid blackout-root causes and dynamics of recent major blackouts.  _IEEE Power and Energy Magazine_ 4, 22–29 (2006).
  * [6] Yang, Y., Nishikawa, T. & Motter, A. E.  Small vulnerable sets determine large network cascades in power grids.  _Science_ 358, eaan3184 (2017).
  * [7] Zhao, J., Li, D., Sanhedrai, H., Cohen, R. & Havlin, S.  Spatio-temporal propagation of cascading overload failures in spatially embedded networks.  _Nature communications_ 7, 1–6 (2016).
  * [8] Contreras, M. G. A. & Fagiolo, G.  Propagation of economic shocks in input-output networks: A cross-country analysis.  _Physical Review E_ 90, 062812 (2014).
  * [9] Boccaletti, S., Latora, V., Moreno, Y., Chavez, M. & Hwang, D.-U.  Complex networks: Structure and dynamics.  _Physics reports_ 424, 175–308 (2006).
  * [10] Crucitti, P., Latora, V. & Marchiori, M.  Model for cascading failures in complex networks.  _Physical Review E_ 69, 045104 (2004).
  * [11] Wang, J.-W. & Rong, L.-L.  Cascade-based attack vulnerability on the us power grid.  _Safety science_ 47, 1332–1336 (2009).
  * [12] Buldyrev, S. V., Parshani, R., Paul, G., Stanley, H. E. & Havlin, S.  Catastrophic cascade of failures in interdependent networks.  _Nature_ 464, 1025–1028 (2010).
  * [13] Schäfer, B., Witthaut, D., Timme, M. & Latora, V.  Dynamically induced cascading failures in power grids.  _Nature communications_ 9, 1–13 (2018).
  * [14] Heide, D., Schäfer, M. & Greiner, M.  Robustness of networks against fluctuation-induced cascading failures.  _Physical Review E_ 77, 056103 (2008).
  * [15] Wang, W.-X. & Chen, G.  Universal robustness characteristic of weighted networks against cascading failure.  _Physical Review E_ 77, 026101 (2008).
  * [16] Huang, X., Gao, J., Buldyrev, S. V., Havlin, S. & Stanley, H. E.  Robustness of interdependent networks under targeted attack.  _Physical Review E_ 83, 065101 (2011).
  * [17] Zeng, A. & Liu, W.  Enhancing network robustness against malicious attacks.  _Physical Review E_ 85, 066130 (2012).
  * [18] Motter, A. E.  Cascade control and defense in complex networks.  _Physical Review Letters_ 93, 098701 (2004).
  * [19] Hu, F., Yeung, C. H., Yang, S., Wang, W. & Zeng, A.  Recovery of infrastructure networks after localised attacks.  _Scientific reports_ 6, 1–10 (2016).
  * [20] Smolyak, A., Levy, O., Vodenska, I., Buldyrev, S. & Havlin, S.  Mitigation of cascading failures in complex networks.  _Scientific reports_ 10, 1–12 (2020).
  * [21] Kornbluth, Y. _et al._ Network overload due to massive attacks.  _Physical Review E_ 97, 052309 (2018).
  * [22] Artime, O. & De Domenico, M.  Abrupt transition due to non-local cascade propagation in multiplex systems.  _New Journal of Physics_ 22, 093035 (2020).
  * [23] Witthaut, D. & Timme, M.  Nonlocal effects and countermeasures in cascading failures.  _Physical Review E_ 92, 032809 (2015).
  * [24] Watts, D. J.  A simple model of global cascades on random networks.  _Proceedings of the National Academy of Sciences_ 99, 5766–5771 (2002).
  * [25] Carreras, B. A., Newman, D. E. & Dobson, I.  North american blackout time series statistics and implications for blackout risk.  _IEEE Transactions on Power Systems_ 31, 4406–4414 (2016).
  * [26] Nesti, T., Sloothaak, F. & Zwart, B.  Emergence of scale-free blackout sizes in power grids.  _Physical Review Letters_ 125, 058301 (2020).
  * [27] Pahwa, S., Scoglio, C. & Scala, A.  Abruptness of cascade failures in power grids.  _Scientific reports_ 4, 1–9 (2014).
  * [28] Spiewak, R., Soltan, S., Forman, Y., Buldyrev, S. V. & Zussman, G.  A study of cascading failures in real and synthetic power grid topologies.  _Network Science_ 6, 448–468 (2018).
  * [29] Burkholz, R., Herrmann, H. J. & Schweitzer, F.  Explicit size distributions of failure cascades redefine systemic risk on finite networks.  _Scientific reports_ 8, 1–8 (2018).
  * [30] Kornbluth, Y., Cwilich, G., Buldyrev, S. V., Soltan, S. & Zussman, G.  Distribution of blackouts in the power grid and the motter and lai model.  _Physical Review E_ 103, 032309 (2021).
  * [31] Valdez, L. D. _et al._ Cascading failures in complex networks.  _Journal of Complex Networks_ 8, cnaa013 (2020).
  * [32] Christie, R.  Power systems test case archive. university of washington, department of electrical engineering (1993).
  * [33] Kendall, M. G.  A new measure of rank correlation.  _Biometrika_ 30, 81–93 (1938).
  * [34] Pei, S., Muchnik, L., Andrade, J. S., Jr., Zheng, Z. & Makse, H. A.  Searching for superspreaders of information in real-world social media.  _Scientific Reports_ 4, 5547 (2014).
  * [35] Hanley, J. A. & Mcneil, B. J.  The meaning and use of the area under a receiver operating characteristic (roc) curve.  _Radiology_ 143, 29 (1982).
  * [36] Rossi, R. A. & Ahmed, N. K.  The network data repository with interactive graph analytics and visualization.  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ (2015).
  * [37] Kunegis, J.  Konect: The koblenz network collection.  In _Proceedings of the 22nd International Conference on World Wide Web_ , WWW ’13 Companion, 1343–1350, DOI: 10.1145/2487788.2488173 (Association for Computing Machinery, New York, NY, USA, 2013).
  * [38] Jeong, H., Mason, S. P., Barabási, A.-L. & Oltvai, Z. N.  Lethality and centrality in protein networks.  _Nature_ 411, 41–42 (2001).

## Acknowledgements

This work was supported by National Natural Science Foundation of China
(71731002).

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.


# Explainable Recommendation Systems by Generalized Additive Models with
Manifest and Latent Interactions

Yifeng Guo1, Yu Su2, Zebin Yang1, and Aijun Zhang1,
1Department of Statistics and Actuarial Science, The University of Hong Kong,
Pokfulam Road, Hong Kong SAR, China
2AI Lab, Suoxinda Holdings Limited, Shenzhen, Guangdong, China Corresponding
author. Email<EMAIL_ADDRESS>
###### Abstract

In recent years, the field of recommendation systems has attracted increasing
attention to developing predictive models that provide explanations of why an
item is recommended to a user. The explanations can be either obtained by
post-hoc diagnostics after fitting a relatively complex model or embedded into
an intrinsically interpretable model. In this paper, we propose the
explainable recommendation systems based on a generalized additive model with
manifest and latent interactions (GAMMLI). This model architecture is
intrinsically interpretable, as it additively consists of the user and item
main effects, the manifest user-item interactions based on observed features,
and the latent interaction effects from residuals. Unlike conventional
collaborative filtering methods, the group effect of users and items are
considered in GAMMLI. It is beneficial for enhancing the model
interpretability, and can also facilitate the cold-start recommendation
problem. A new Python package GAMMLI is developed for efficient model training
and visualized interpretation of the results. By numerical experiments based
on simulation data and real-world cases, the proposed method is shown to have
advantages in both predictive performance and explainable recommendation.

Keywords: Generalized additive model; Interaction effects; Matrix
factorization; Model interpretability; Personalized recommendation; Cold start
problem.

##  1 Introduction

In the last decades, the recommendation system (RecSys) that predicts user
preferences on multiple items has been widely adopted in many fields,
including banking, the Internet, E-commerce, and social media. Traditional
recommendation systems mainly fall into two categories: content-based and
collaborative filtering (CF)-based recommendation. The former originates from
information retrieval and information filtering Adomavicius and Tuzhilin,
(2005), and it models user and item profiles with assorted available content
information. Two examples of content-based recommendation are the term
frequency/inverse document frequency (TF-IDF) that specifies keyword weights
(Salton,, 1989) and the Fab system that recommends Web pages to users
(Balabanović and Shoham,, 1997). The CF-based recommendation systems are
developed for matching users with similar tastes and preferences. User-based
CF for the GroupLens news recommendation system is one of the earliest CF
approaches by Resnick et al., (1994) and later the item-based CF method is
discussed by Sarwar et al., (2001). Popularized by Netflix competition, the
latent factor model based on matrix factorization (MF; Koren et al.,, 2009)
became a widely known CF method with relatively high prediction performance.
The RecSys research community has extended MF techniques in various aspects,
including combining multiple matrices with network relationships among the
variables (Zhang et al.,, 2011), integrating MF with topic models (Wang et
al.,, 2015) and integrating side information using matrix co-factorization
technique (Pujahari and Sisodia,, 2020). Both content-based and CF-based
methods have pros and cons. In particular, the CF techniques do not model the
user or item features as carefully as the content-based methods. Combining
these two types of methods leads to the hybrid recommendation system. Among
others, the most successful hybrid is the factorization machine (FM; Rendle,,
2010, 2012) that extends matrix factorization for generic modeling of
features. In recent years, the RecSys technique has been enhanced through deep
learning in various ways, including wide and deep network (Cheng et al.,,
2016), DeepFM (Guo et al.,, 2017), neural collaborative filtering (He et al.,,
2017) and multi-criteria collaborative filtering (Nassar et al.,, 2020) which
all utilize the complex neural networks to richly express the user and item
features and their dummy IDs.

The state of the art recommendation systems includes MF, FM, and DeepFM
reviewed above, which have superb predictive performance. However, they merely
make recommendations, but can hardly tell why an item is recommended to a user
in terms of user or item features. In this sense, the MF, FM, and DeepFM
methods are usually referred to as black-box models, which are known to be
lacking transparency, persuasiveness, effectiveness, trustworthiness, and user
satisfaction. There exist two categories of explainable recommendations in the
literature. One is the post-hoc diagnostics that generates explanations after
the recommendation has been made; see e.g., the explanation mining by Peake
and Wang, (2018). The other category is intrinsic embedding that tries to
embed explicit interpretability into the recommendation models. Examples
include the explicit factor model based on phrase-level sentiment analysis
(Zhang et al.,, 2014) and the neural attentional regression model with review-
level explanations by Chen et al., (2018). Both examples are supplemented by
user reviews on each item through text analysis. Such user-item response data
together with textual review collection can be also modeled by the tensor
analysis framework; see e.g., Wang et al., (2018). A more comprehensive review
of such sort of explainable recommendations can be referred to Zhang and Chen,
(2020).

In this paper, we develop a new kind of intrinsically interpretable RecSys
from the statistical modeling point of view, in particular, based on the idea
of the generalized additive model (GAM; Hastie and Tibshirani,, 1990).
Consider a special case of RecSys problem with multiple users but only one
item (either label or rating), then it corresponds to the standard machine
learning problem. For data collected with user features
$\bm{x}\in\mathbb{R}^{p}$ and a single item response $y\in\mathbb{R}$, the GAM
takes the following form

$g(\mathbb{E}[y])=\mu+h_{1}(x_{1})+\cdots+h_{p}(x_{p}),$ (1)

where $\mu$ is the intercept, $h(\cdot)$ is a smooth nonparametric function,
and $g$ is a pre-specified link function. In the fast-emerging area of
interpretable machine learning, GAM is known to process a good balance between
prediction performance and inherent interpretability. Such a statistically
interpretable model is a special case of additive index model (AIM) that has
been successfully re-deployed into the explainable neural network (xNN)
architectures (Vaughan et al.,, 2018; Yang et al.,, 2020) through univariate-
input subnetwork parametrization of $h(\cdot)$ in (1). It has been further
extended to the GAMI-Net recently proposed by Yang et al., (2020) to include
structured pairwise interaction effects through bivariate-input subnetworks.
Meanwhile, it draws our attention that GAM plays an essential role in the
recent work of Agarwal et al., (2020) under the rephrased name of neural
additive models, which corresponds to a special case of xNN and GAMI-Net.

The RecSys problem has multiple item responses per user, and the data
collection includes user features $\bm{x}\in\mathbb{R}^{p}$, item features
$\bm{z}\in\mathbb{R}^{q}$, and usually incomplete (user, item)-responses. It
is our goal to find an appropriate way of applying GAM in (1) for modeling
such RecSys data of user-item features and responses. It is straightforward to
model the main effects based on either user or item features. For the RecSys
problem, the user-item cross interactions play critical roles in predicting
the (user, item)-responses. We propose a two-way approach to model such user-
item interactions, including an explicit way of modeling the manifest
interactions induced from the observed user and item features, and an implicit
way of modeling the latent interactions induced from other unobserved
features. This leads to generalized additive models with manifest and latent
interactions (GAMMLI) for explainable recommendation systems. Moreover, for
pursuing enhanced model interpretability, the latent interactions are modeled
with user and item group structures. The details about the GAMMLI methodology
and how to interpret GAMMLI results are presented in Section 2.

The computational methods for GAMMLI estimation are presented in Section 3,
where we discuss the stage-wise network training algorithm, latent interaction
estimation, and hyperparameter tuning. In Section 4, we conduct numerical
experiments based on both synthetic data and three real-world cases.It is
shown that the proposed method has not only high prediction accuracy but also
superior model interpretability. Finally, Section 5 concludes the paper with
directions for future work.

##  2 GAMMLI Methodology

Assume a RecSys data consists of $m$ users (indexed by $i=1,\ldots,m$) and $n$
items (indexed by $j=1,\ldots,n$). The observed data includes user features
$\bm{x}_{i}\in\mathbb{R}^{p}$, item features $\bm{z}_{j}\in\mathbb{R}^{q}$,
and (user, item)-responses by $y_{ij}$. Based on (1), we propose a novel
explainable RecSys framework which provides not only the intuitive
interpretation of observed features but also the heuristic interpretation of
unobserved features. GAMMLI consists of three components, i.e., main effects,
manifest interactions, and latent interactions.

  * •

Main user and item effects: the main effects can be modeled by (1) based on
user features $\bm{x}$ and item features $\bm{z}$.

  * •

Manifest user-item interactions: the manifest interaction effects
$I^{(m)}_{ij}$ are induced by observed user features $\bm{x}$ and item
features $\bm{z}$. They can be modeled explicitly by
$I^{(m)}_{ij}=\sum_{a,b}h(x_{ia},z_{jb})$ for selective $a\in\\{1,\ldots,p\\}$
and $b\in\\{1,\ldots,q\\}$.

  * •

Latent user-item interactions: the latent interaction effects $I^{(l)}_{ij}$
are induced by other unobserved user or item features. They can be modeled
implicitly through low-rank approximation
$I^{(l)}_{ij}\approx\sum_{c=1}^{r}u_{ic}v_{jc}$ for a relatively small integer
$r$.

Putting these effects together, we have the GAMMLI model of the form

$g(\mathbb{E}[y_{ij}])=\mu+\sum_{a=1}^{p}h_{a}^{(u)}(x_{ia})+\sum_{b=1}^{q}h_{b}^{(v)}(z_{jb})+\sum_{a,b}h^{(u,v)}_{a,b}(x_{ia},z_{jb})+\sum_{c=1}^{r}u_{ic}v_{jc},$
(2)

where the link function $g$ is specified according to different types of
(user, item)-responses.

Figure 1: The network architecture of GAMMLI. In the left part, the network
fits $p$ user main effects and $q$ item main effects based on explicitly
observed features. In the middle part, the user-item manifest interactions are
modeled by bivariate network inputs. The right part presents the group-
constrained low-rank approximation for latent user-item interactions. These
three modules are additively combined to form the final response.

The model architecture of GAMMLI is shown in Fig. 1. In particular, GAMMLI is
fitted in a sequential way with three stages. In the first and second stages,
the explicitly observed user and item features are used to model main effects
and manifest user-item interactions. In the third stage, the residuals from
the first two stages are used to fit the latent user-item interactions. More
details about such a sequential strategy are discussed in Section 3.

###  2.1 Main Effects and Manifest Interactions

Consider the three $h(\cdot)$ terms in (2), which are usually fitted by
nonparametric methods like splines or trees. Recently, (Yang et al.,, 2020)
proposed a GAMI-Net based on additive subnetworks with better accuracy and
interpretability than multiple existing machine learning methods. Such neural
network parametrization is incorporated in this paper for modeling smooth
$h(\cdot)$ functions in (2). Each main effect function, either
$h_{a}^{(u)}(x_{a})$ of a user feature, or $h_{b}^{(v)}(z_{b})$ of an item
feature, can be represented by a univariate-input subnetwork. Each user-item
manifest interaction $h^{(u,v)}_{a,b}(x_{a},z_{b})$ can be represented by a
bivariate-input subnetwork. For the sake of identifiability, each main effect
or manifest interaction is assumed to have mean zero. The categorical features
can be modeled by bias nodes, where each bias node captures the intercept
effect of a corresponding dummy variable. More details about the subnetwork
configuration can be referred to Yang et al., (2020).

Manifest user-item interactions are fitted to the residuals from the main
effect prediction. Each manifest interaction subnetwork consists of two input
nodes, one for a user feature and the other for an item feature. Note that if
we directly use GAMI-Net (Yang et al.,, 2020), all the user-user, user-item,
and item-item interactions would be taken into account. However, in this
paper, we restrict to the user-item cross manifest interactions only, since
they are more interesting to model the interactions among user-item responses.
On the other hand, we consider all two-way user-item manifest interactions
without imposing weak heredity constraint or interaction filtering.

###  2.2 Latent Interaction Effects

Latent user-item interactions $I^{(l)}_{ij}$ are used to model the unobserved
features after fitting the main effect and manifest interactions. Denote
$\bm{M}_{m\times n}$ as the residual matrix with missing values. The observed
entries are indexed by $\Omega=\\{(i,j):\bm{M}_{ij}\mbox{ is not missing}\\}$.
It can be approximated by low-rank matrix factorization as
$I^{(l)}_{ij}\approx\sum_{c=1}^{r}u_{ic}v_{jc}$ for each user $i$ and each
item $j$. Here, the column vectors $\\{\bm{u}_{c},\bm{v}_{c}\\}$ for
$c=1,\ldots,r$ are called user latent factors and item latent factors,
respectively. Using matrix notations, we can write them as
$\bm{U}=[\bm{u}_{1},\ldots,\bm{u}_{r}]$ of size $m\times r$ and
$\bm{V}=[\bm{v}_{1},\ldots,\bm{v}_{r}]$ of size $n\times r$.

Unlike the existing matrix factorization methods, we take the group effect
into consideration. The users and items are simultaneously clustered based on
explicitly observed features with a biclustering technique. The group
constraint can increase the quality of recommendation in heterogeneous groups.
Denote by $S_{k}^{(u)}$ for $k=1,\ldots,K$ the $K$ disjoint clusters for the
users, and by $S_{l}^{(v)}$ for $l=1,\ldots,L$ the $L$ disjoint clusters of
the items. For each user-item pair $(i,j)$, let $k(i)$ be the corresponding
cluster index for user $i$ and let $l(j)$ be the corresponding cluster index
for item $j$. Then, we can find the cluster centroids

$\widetilde{\bm{U}}_{i\ast}=\frac{1}{\big{|}S_{k(i)}^{(u)}\big{|}}\sum_{i^{\prime}\in
S_{k(i)}^{(u)}}\bm{U}_{i^{\prime}\ast}\ \mbox{ and }\
\widetilde{\bm{V}}_{j\ast}=\frac{1}{\big{|}S_{l(j)}^{(v)}\big{|}}\sum_{j^{\prime}\in
S_{l(j)}^{(v)}}\bm{V}_{j^{\prime}\ast},\quad\forall(i,j),$ (3)

so as to form the cluster centroid matrices $\widetilde{\bm{U}}$ and
$\widetilde{\bm{V}}$. Then, we suggest to estimate the latent factors by the
following regularized least squares,

$\min_{\bm{U},\bm{V}}\|P_{\Omega}(\bm{M}-\bm{U}\bm{V}^{T})\|_{F}^{2}+\lambda\left(\|\bm{U}-\widetilde{\bm{U}}\|_{F}^{2}+\|\bm{V}-\widetilde{\bm{V}}\|_{F}^{2}\right),$
(4)

where $P_{\Omega}$ is the projection of a matrix that preserves the entries
indexed by $\Omega$, the subscript $F$ denotes the Frobenius norm of a matrix,
and $\lambda\geq 0$ is the regularization strength.

The problem (4) can be solved by alternating least squares similar to the
SoftImpute algorithm (Mazumder et al.,, 2010; Hastie et al.,, 2015). Given the
pilot estimates for $\hat{\bm{U}},\hat{\bm{V}}$, and pre-specified clusters,
we may obtain the cluster centroid matrices by (3), then alternatively
estimate $\bm{V}$ by

$\min_{\bm{V}}\|P_{\Omega}(\bm{M}-\hat{\bm{U}}\bm{V}^{T})\|_{F}^{2}+\lambda\|\bm{V}-\widetilde{\bm{V}}\|_{F}^{2},$

and estimate $\bm{U}$ with updated $\hat{\bm{V}}$ by

$\min_{\bm{U}}\|P_{\Omega}(\bm{M}-\bm{U}\hat{\bm{V}}^{T})\|_{F}^{2}+\lambda\|\bm{U}-\widetilde{\bm{U}}\|_{F}^{2},$

each corresponding to a multi-response ridge regression with closed-form
solutions. A detailed procedure based on soft-imputed singular value
decomposition and alternative least squares is presented in Algorithm 1.

Algorithm 1 Latent Interactions-ALS

1:Initialize residual matrix $\bm{M}$,
$\hat{\bm{U}}=\bm{U}^{\ast}\bm{\Sigma},\hat{\bm{V}}=\bm{V}^{\ast}\bm{\Sigma}$
where $\bm{U}^{\ast}$ is orthonormal matrix, $\bm{\Sigma}=\bm{I}_{r}$ is an
identity matrix and $\bm{V}^{\ast}=0$. Based on the pre-specified clusters,
calculate cluster centroid matrices $\widetilde{\bm{U}},\widetilde{\bm{V}}$.

2:Given
$\hat{\bm{U}}=\bm{U}^{\ast}\bm{\Sigma},\hat{\bm{V}}=\bm{V}^{\ast}\bm{\Sigma},\widetilde{\bm{U}},\widetilde{\bm{V}}$,
solve:
$\min_{\bm{V}}\|P_{\Omega}(\bm{M}-\hat{\bm{U}}\bm{V}^{T})\|_{F}^{2}+\lambda\|\bm{V}-\widetilde{\bm{V}}\|_{F}^{2},$

3:Let $\ddot{\bm{V}}=\bm{V}-\widetilde{\bm{V}}$ and
$\bm{M}^{\ast}=P_{\Omega}(\bm{M}-\hat{\bm{U}}\hat{\bm{V}}^{T})+\hat{\bm{U}}(\hat{\bm{V}}-\widetilde{\bm{V}})^{T}$.
We update $\hat{\bm{V}}$ by solving:
$\min_{\ddot{\bm{V}}}\|\bm{M}^{\ast}-\hat{\bm{U}}\ddot{\bm{V}}^{T}\|_{F}^{2}+\lambda\|\ddot{\bm{V}}\|_{F}^{2},$
with solution
$\ddot{\bm{V}}^{T}=(\bm{\Sigma}^{2}+\lambda\bm{I})^{-1}\bm{\Sigma}{\bm{U}^{\ast}}^{T}\bm{M}^{\ast}.$

4:$\hat{\bm{V}}=\ddot{\bm{V}}+\widetilde{\bm{V}}$. Obtain SVD result for
$\hat{\bm{V}}\bm{\Sigma}$ and update by the following steps:

  1. i

update $\bm{V}^{\ast}\leftarrow$ left-singular vectors;

  2. ii

update $\bm{\Sigma}\leftarrow$ singular values;

  3. iii

update $\hat{\bm{V}}\leftarrow\bm{V}^{\ast}\bm{\Sigma}$ and cluster centroid
matrix $\widetilde{\bm{V}}$.

5:By symmetry, solve for $\bm{U}$ with solution
$\ddot{\bm{U}}=\bm{M}^{\ast}\bm{V}^{\ast}\bm{\Sigma}(\bm{\Sigma}^{2}+\lambda\bm{I})^{-1}.$

6:$\hat{\bm{U}}=\ddot{\bm{U}}+\widetilde{\bm{U}}$. Obtain SVD result for
$\hat{\bm{U}}\bm{\Sigma}$ and update by the following steps:

  1. i

update $\bm{U}^{\ast}\leftarrow$ left-singular vectors;

  2. ii

update $\bm{\Sigma}\leftarrow$ singular values;

  3. iii

update $\hat{\bm{U}}\leftarrow\bm{U}^{\ast}\bm{\Sigma}$ and cluster centroid
matrix $\widetilde{\bm{U}}$.

7:Repeat 2-6 until convergence.

8:return reconstruction matrix $\hat{\bm{M}}=\hat{\bm{U}}\hat{\bm{V}}^{T}$.

The group constraints are incorporated into (4) to enhance the
interpretability of the proposed GAMMLI method for explainable RecSys, with
particular benefits in solving cold-start problems. The cold-start problem
targets brand new users or items without existing information about user-item
responses, and they widely exist in the real RecSys applications. Cold-start
problems are hard to be dealt with by collaborative filtering methods since
there are no records to find similar patterns with other users or items. By
the proposed GAMMLI framework, we can take advantage of the auxiliary features
associated with new users or items. Such features can be used not only for
computing main effects and manifest interactions but also for obtaining the
specific cluster assignments. The cluster centroids can be used as the latent
interactions for the cold-start cases. Therefore in GAMMLI, the clustering
procedure bridges the observed features with latent factors and provides a
heuristic interpretation of latent interactions.

###  2.3 GAMMLI Results Interpretation

It is straightforward to interpret the GAMMLI results by additive effect
decomposition (2). We may characterize the contribution of each main effect
and manifest interaction, as well as the contribution by the latent
interactions. Besides, we can visualize the results.

The importance ratio (IR) is a way to quantify the importance by the
percentage of variation explained. For the modeled main effects and manifest
interactions, their variations can be empirically evaluated by the sum of
squared marginal predictions. For the modeled low-rank latent interactions, we
can simply evaluate their sum of squares, too. Converting them into the
percentage scale, we define the IR of each modeled effect as follows,

$\displaystyle{\rm IR}(\hat{h}_{a}^{(u)})$
$\displaystyle=\frac{1}{|\Omega|-1}\sum_{(i,j)\in\Omega}\hat{h}_{a}^{(u)}(x_{ia})^{2}/T\equiv
D(\hat{h}_{a}^{(u)})/T,$ (5) $\displaystyle{\rm IR}(\hat{h}_{b}^{(v)})$
$\displaystyle=\frac{1}{|\Omega|-1}\sum_{(i,j)\in\Omega}\hat{h}_{b}^{(v)}(z_{jb})^{2}/T\equiv
D(\hat{h}_{b}^{(v)})/T,$ $\displaystyle{\rm IR}(\hat{h}^{(u,v)}_{a,b})$
$\displaystyle=\frac{1}{|\Omega|-1}\sum_{(i,j)\in\Omega}\hat{h}^{(u,v)}_{a,b}(x_{ia},z_{jb})^{2}/T\equiv
D(\hat{h}^{(u,v)}_{a,b})/T,$ $\displaystyle{\rm
IR}(\\{\hat{\bm{U}},\hat{\bm{V}}\\})$
$\displaystyle=\frac{1}{|\Omega|-1}\sum_{(i,j)\in\Omega}\Big{(}\hat{\bm{U}}_{i\ast}\hat{\bm{V}}_{j\ast}^{T}-\bar{M}\Big{)}^{2}/T\equiv
D(\\{\hat{\bm{U}},\hat{\bm{V}}\\})/T,$

where $|\Omega|$ denotes the total number of observed entries in the training
data, $\bar{M}$ denotes the mean of estimated latent interactions and
$T=\sum_{a=1}^{p}D(\hat{h}_{a}^{(u)})+\sum_{b=1}^{q}D(\hat{h}_{b}^{(v)})+\sum_{a,b}D(\hat{h}^{(u,v)}_{a,b})+D(\\{\hat{\bm{U}},\hat{\bm{V}}\\})$
represents the total variation.

The partial relationship of the main effects and manifest interactions can be
visualized by univariate line plots and bivariate heatmaps, similar to Yang et
al., (2020). For the latent interactions subject to group structures, we can
also use the heatmap to visualize the user-item cross effects. Examples of
these plots are provided in Section 4.

##  3 Computational Methods

The training of GAMMLI is conducted sequentially in three stages, including
the estimations of main effects, manifest interactions, and latent
interactions.

Step 1: In the main effects training, all main effect subnetworks are
estimated simultaneously, and the user-item interaction effects are fixed to
zero. The trainable parameters in the network are updated by mini-batch
gradient descent with adaptive learning rates determined by the Adam
optimizer. The main effects training step will stop when the maximum training
epochs are reached or the given early stopping criterion is met. Each main
effect is then centered to have a mean zero. The contribution of each main
effect can be quantified by their variations as shown in (5). By selecting
top-$s_{1}$ main effect subnetworks, nonsignificant main effects are
accordingly removed. Note the optimal value of $s_{1}$ is automatically
selected in the modeling process subject to the validation loss.

Step 2: After the main effects training finished, the subnetworks of pairwise
user-item interactions are trained using the prediction residuals of the first
stage. All user-item manifest interactions are taken into account due to the
hierarchical principle. Each estimated manifest interaction is centered to
have mean zero, for which the offset is added to the bias node in the output
layer. The manifest interactions training will stop similarly as in the first
stage. The manifest interaction effects are also pruned by selecting the
top-$s_{2}$ interaction effects according to $D(\hat{h}^{(u,v)}_{a,b})$ in (5)
for selective $a\in\\{1,\ldots,p\\}$ and $b\in\\{1,\ldots,q\\}$. Also note
there is no need to manually specify $s_{2}$, as it is automatically selected
using the validation loss. Then, the main effects and manifest interactions
will be jointly fine-tuned to improve performance. More details can be
referred to Yang et al., (2020).

Step 3: Before the latent interactions are estimated, all users and items will
be assigned into $K$ and $L$ clusters, respectively. Various clustering
algorithms have been proposed in the machine learning literature. And in fact,
any clustering algorithm can be employed here. For simplicity, we just
incorporate the popularly used K-means clustering. The clustering results,
i.e., $S_{k}^{(u)}$ ($k=1,\ldots,K$) for users and $S_{l}^{(v)}$
($l=1,\ldots,L$) for items will be used for the next step modeling.

Step 4: The prediction residuals of the first/second stage are used to fit
latent interaction effects as shown in Algorithm 1, which alternatively solves
$\hat{\bm{U}}$ and $\hat{\bm{V}}$ by solving a multi-response ridge
regression. For regression tasks, residuals are calculated by
$y-\hat{y_{1}}-\hat{y_{2}}$. In classification scenarios, residuals are
calculated as is in gradient boosting machine (GBM;Friedman, (2001)). That is,
we first transform the binary response from $\\{0,1\\}$ to $\\{-1,1\\}$. The
residual is then computed by
$\tilde{y}=2y/(1+\exp(2y_{i}(\hat{y_{1}}+\hat{y_{2}})))$, where $y$ is true
response, and $\hat{y_{1}},\hat{y_{2}}$ are estimated responses in the first
and second stages, respectively. Finally, The training of GAMMLI is summarized
in Algorithm 2.

Algorithm 2 GAMMLI Training

1:Train and select the top-$s_{1}$ main effect subnetworks.

2:Train and select the top-$s_{2}$ manifest interaction subnetworks.

3:Cluster (K-means) user and item features into $K$ and $L$ clusters,
respectively.

4:Train latent interaction effects via Algorithm 1.

Hyperparameter Tuning. There exist several hyperparameters in GAMMLI, which
are critical for the final prediction performance. The GAMMLI is configured
with the following empirical settings. The subnetwork structure is configured
to [20, 10] with hyperbolic tangent nodes. For the neural network training,
the initial learning rate is set to 0.001, and the number of tuning epochs is
fixed to 200. The epochs of training main effects, manifest interactions, and
latent interactions are all set to 1000. A 20% validation set is split for
early stopping in each stage, and the early stopping threshold is set to 100
epochs. For interaction training, the rank of latent interactions $r$ is set
to 3.

Three important hyperparameters are identified in the latent interaction
estimation stage, including users group number $K$, items group number $L$,
and regularization strength $\lambda$. To pursue better performance, a coarse-
to-fine grid-search method is employed to optimize $(K,L,\lambda)$.
Specifically, the initial search space is set as $K$ (2 to 30), $L$ (2 to 30),
and $\lambda$ (0 to 50). The search space is evenly partitioned such that each
hyperparameter has 5 grid points. We enumerate all the possible 125
hyperparameter combinations and select the best performing one. Then, the
second-iteration search is conducted sequentially by repeatedly centering on
the best performing grid point, reduce the search space, and scatter high-
granularity grid points. Such a procedure will be repeated and will be stopped
after the fifth iteration. After hyperparameter optimization, a suitable
combination of $(K,L,\lambda)$ can be obtained.

##  4 Numerical Experiments

In this section, we evaluate our model via two simulation scenarios and three
real datasets in offline settings. We begin by introducing the experimental
setup and then report and analyze the experimental results.

Four competitive benchmark models are introduced for comparison, including
eXtreme Gradient Boosting (XGBoost), singular value decomposition (SVD),
factorization machine (FM), and DeepFM. Four compared models cover different
types of recommendation systems, i.e., machine learning model (XGBoost),
collaborative filtering model (SVD), and complicated monolithic method (FM &
DeepFM).

Each dataset is split into training, validation, and test sets. In SVD, we set
the rank to 5. For both FM and DeepFM, we set different loss types for
different tasks, i.e., mean squared error (MSE) for regression tasks and Log
Loss for classification tasks. The number of training epochs is set
differently for different tasks. In XGBoost, the maximum tree depth is
selected from {3,4,5,6,7,8}. Other hyperparameters of benchmark models are set
to their default values. All experiments are repeated 10 times, where we use
different random train test split for modeling fitting and evaluation. For
each dataset, we evaluate the test set root mean squares error (RMSE) and mean
absolute error (MAE) for regression tasks and the area under the ROC curve
(AUC) and Log Loss for binary classification tasks.

All the experiments are implemented using the Python environment. The proposed
GAMMLI is implemented using the TensorFlow 2.0 platform, with the open-sourced
package uploaded to the Github111https://github.com/SelfExplainML/GAMMLI. The
SVD is implemented using surprise package222http://surpriselib.com/. FM and
DeepFM are implemented by the TensorFlow
package333https://github.com/ChenglongChen/tensorflow-DeepFM and XGBoost is
implemented by the popular xgboost package444https://github.com/dmlc/xgboost.

###  4.1 Simulation Study

A synthetic data is first generated, where the numbers of users and items are
set to $m=1000$ and $n=1000$, respectively. Each user or item has 5 numerical
features and a latent vector with rank 3. The user/item features and latent
vectors are generated from isotropic Gaussian blobs. The features are scaled
within $[0,1]$ and latent vectors are scaled within $[-1,1]$. The actual
numbers of user and item groups are both set to 10. To mimic the real data
case where only a small percentage of responses is collected, the response
matrix is generated as a sum of feature effect and latent interactions with
size $m\times n$ and a missing rate $\omega=0.9$. A noise term generated from
the standard normal distribution is also included when generating the response
matrix. Both regression and classification scenarios will be included in the
simulation study.

In the regression scenario, the feature effects of user and item are generated
as follows:

$\displaystyle y=5x_{1}+5z_{1}^{2}+0.5e^{-4(z_{2}+x_{3})+4}+5sin(2\pi
x_{2}z_{3})+3\sum_{k=1}^{3}u_{k}v_{k}^{T}+\epsilon,$ (6)

where $x,z$ represent user and item features and
$\sum_{k=1}^{3}u_{k}v_{k}^{T}$ denotes the sum of latent interactions. This
case exactly follows the model formulation of GAMMLI. It can be observed that
each user or item has two noisy features. As for the classification scenario,
the formulation follows the regression setting but the response is binarized
as $y=\mathbb{I}_{\\{y>0.5\\}}$, where $\mathbb{I}$ is an indicator function.

Fig. 2 shows the training and validation losses of the simulation study. It
can be observed that the losses drop significantly as manifest interactions
and latent interactions are added to the model. Each component of GAMMLI
contributes greatly to model accuracy.

(a) Regression (b) Classification

Figure 2: Training and validation losses of the simulation study.

Fig. 3 draws the global interpretation of main effects, manifest interactions
of the simulation study (regression scenario). Note that in the original
formulation, it is assumed to have 6 active main effects $\\{x_{1},z_{1}\\}$
and 2 active interaction effects $\\{(x_{2},z_{3}),(x_{3},z_{2})\\}$. We can
equivalently separate marginal effects from the interactions so that the
active main effects also include $\\{x_{2},x_{3},z_{2},z_{3}\\}$. It can be
seen that GAMMLI captures all the important feature effects.

Fig. 4 presents the global interpretation of latent interaction effects with
the IR of 14.63%. The user group number $K$ (the y-axis) and item group number
$L$ (the x-axis) are set to 10. From Fig. 4, we can find that user groups 7
and 9 are most similar, as they have similar margin effects and similar
preferences for different item groups. Moreover, user group 5 has the
strongest preference of item group 8 among all group interactions.

Except for the global interpretation, GAMMLI can also be locally explained,
leading to a transparent personal recommendation system. Given a user and an
item, the model outputs the final decision and why we recommend this item to
this user. The additive components of main effects, manifest interactions, and
latent interactions are provided. We can also get the local interpretation of
the simulation study (regression scenario) from the prediction diagnosis of
one sample point in Fig. 5. This information will be useful to understand the
recommendation results.

(a) Ground truth (b) GAMMLI

Figure 3: Main effects and manifest interactions of the simulation regression
setting. Figure 4: Global interpretation of latent interactions (simulation
regression).

Figure 5: Local interpretation for one randomly picked sample (simulation
regression).

Table 1 reports the average test accuracy of our model and several competitive
black-box recommendation system models under both regression and
classification scenarios. For each scenario, the best model is highlighted in
bold. It can be observed that our proposed GAMMLI performs the best under both
regression and classification scenario. While maintaining superior prediction
accuracy, we achieve much better interpretability, which leads our model more
transparent and understandable.

Table 1: Testing accuracy comparison under simulation study. Method | Regression | Classification
---|---|---
MAE | RMSE | AUC | Log Loss
GAMMLI |  0.3613$\pm$0.0998 |  0.5405$\pm$0.1357 |  0.9760$\pm$0.0022 |  0.1972$\pm$0.0091
FM | 0.5239$\pm$0.0813 | 0.7208$\pm$0.1094 | 0.9389$\pm$0.0115 | 0.3599$\pm$0.0393
DeepFM | 0.4449$\pm$0.0194 | 0.6149$\pm$0.0292 | 0.9437$\pm$0.0177 | 0.3396$\pm$0.0566
XGBoost | 0.4574$\pm$0.0040 | 0.6360$\pm$0.0058 | 0.9223$\pm$0.0003 | 0.4121$\pm$0.0024
SVD | 1.6215$\pm$0.0000 | 2.0697$\pm$0.0000 | 0.7795$\pm$0.0000 | 0.6610$\pm$0.0000

###  4.2 Real Data Applications

We further consider three real-world datasets, MovieLens-1m dataset, the
Santander bank dataset, and the University Recommendation dataset.

Movielens 555https://grouplens.org/datasets/movielens/1m/ is a movie rating
dataset that has been widely used to evaluate recommendation system
algorithms. We used the version Movielens-1m, i.e., it contains one million
ratings. Several features are taken into consideration, including gender, age,
occupation, and genres. There are 6040 users and 3706 movies in this dataset,
and the sparsity is 95.5%.

Santander bank 666https://www.kaggle.com/c/santander-product-recommendation is
an open Kaggle dataset. The dataset provides 1.5 years of customer behavior
data from Santander bank to predict what new products the customers will
purchase on 2016-06-28. In the Santander dataset, 22 bank products and 24
profile features of users are taken into consideration. The 22 bank products
are further summarized into three categories, which form a new item feature,
item property. To remove the effect of time and make it convenient for
modeling, we only collect samples where users have purchase behavior on
2016-06-28 as positive samples. We set the records that users without
purchased items as the initial negative samples. To distinguish the real
negative samples and the unknown samples, we first regard the negative samples
as missing values and use the SoftImpute algorithm to fill the missing values.
After imputing the missing values, we set a proper threshold to define
reliable negative samples (Zhang and Lee,, 2005).

University Recommendation
777https://www.kaggle.com/nitishabharathi/university-recommendation is an open
Kaggle dataset. It contains over 14000 students’ profiles in computer science-
related fields with admits/rejects to 54 different universities in the USA.
Universities’ rankings and public status are manually collected from US News
888https://www.usnews.com/best-colleges as item features.

Table 2: Testing accuracy comparison under real data. Method | Movielens | Santander | University Recommendation
---|---|---|---
MAE | RMSE | AUC | Log Loss | AUC | Log Loss
GAMMLI |  0.7095$\pm$0.0047 |  0.9187$\pm$0.0038 |  0.9954$\pm$0.0006 |  0.0463$\pm$0.0026 |  0.7486$\pm$0.0064 |  0.6304$\pm$0.0093
FM |  0.6755$\pm$0.0024 |  0.8584$\pm$0.0013 |  0.9901$\pm$0.0016 |  0.1025$\pm$0.0052 |  0.7046$\pm$0.0420 |  0.7105$\pm$0.1155
DeepFM |  0.6748$\pm$0.0022 |  0.8578$\pm$0.0008 |  0.9899$\pm$0.0020 |  0.1030$\pm$0.0066 |  0.7207$\pm$0.0288 |  0.6801$\pm$0.0619
XGBoost |  0.8634$\pm$0.0000 |  1.0622$\pm$0.0000 |  0.9947$\pm$0.0002 |  0.0507$\pm$0.0010 |  0.7780$\pm$0.0047 |  0.5647$\pm$0.0048
SVD |  0.9041$\pm$0.0000 |  1.1373$\pm$0.0000 |  0.9867$\pm$0.0000 |  0.0851$\pm$0.0000 |  0.5110$\pm$0.0000 |  0.6924$\pm$0.0000

As shown in Table 2, GAMMLI shows comparative predictive performance to that
of other benchmark models. It can be observed that GAMMLI outperforms other
benchmark models in the Santander dataset and keeps a small disparity with the
most accurate benchmark models in Movielens and University Recommendation.
Moreover, other benchmark models are not that stable among the three datasets.
GAMMLI preserves relatively advantageous performance, which also shows its
robustness. In practice, GAMMLI is more likely to have better predictive
performance when the observed features have strong signals. While maintaining
relatively superior prediction accuracy, a clear interpretation of GAMMLI is
generated. We take the University Recommendation dataset as an example to
demonstrate the interpretability of GAMMLI.

The global interpretation of the main effects of the University Recommendation
dataset is shown in Fig. 6. From the univariate plots, we can get four
significant main effects in University Recommendation which are ranking
(University US ranking), cgpa, PublicOrNot, major with IRs equal to 26.4%,
9.2%, 6.1%, 3.4%, respectively. In addition to the univariate plots, the
density plots of the corresponding variables are also provided. The most
significant item feature is ranking. It can be observed that students would
have a higher chance to be admitted if they apply to universities with lower
rankings. For those top universities, though many students submit their
applications, most of them got rejected. Another significant item feature is
PublicOrNot. The admission rate of private universities is around 60%, while
public universities’ admission rate is only around 50%. The most significant
user feature is cgpa. The high score of cumulative GPA has a monotonically
positive impact on admission decisions.

Regarding manifest interactions, it can be observed that two manifest
interactions are captured, which are greV vs. ranking and major vs. ranking
with IRs equal to 9.5% and 3.8%. As the main effects are removed due to the
marginal clarity constraint, the estimated pairwise interactions can be
individually illustrated. The most significant manifest interactions is greV
vs. ranking. For the students with GRE verbal higher than around 155, the
students may be more likely to be admitted by the Top 50 universities. In
practice, for students major in computer science-related subjects, most of
them are skilled in mathematics and get high scores in the GRE quantitative
test. A high GRE verbal score might be equivalent to a high GRE score.

Figure 6: Global interpretation of main effects and manifest interactions
(University Recommendation).

Fig. 7 shows the global interpretation of latent interaction effects are
significantly important with the IR equals to 41.48%. The user group number
$K=11$ and item group number $L=3$ are chosen by the validation set. It can be
observed that user groups 6 and 10 are most similar among all user groups for
they have similar sidebar colors and similar preferences for different item
groups. The strongest preference is obtained between user group 0 and item
group 0 among all latent interactions. We use radar plot Fig. 8 to portray the
centroids of similar groups in terms of their profiles. From fig. 8, we can
find that the difference between user group 6 and user group 10 mainly
attributes to ugCollege which refers to as applicants’ undergraduate college.
Note ugCollege is not an important feature from the analysis of main and
manifest interactions.

Figure 7: Global interpretation of latent interactions (University
Recommendation). Figure 8: Radar plot for user groups 6 and 10 (University
Recommendation).

###  4.3 Cold Start Problem

When a user or item is input into the system, it is automatically judged
whether it is new. As for a regular user or item, the system will assign the
existing latent vectors to it. If it is a new one, only the user/item features
are given, and there will be no existing latent vectors. Therefore, we can
find its nearest feature group and assign it with the group’s centroid of
latent factors. Finally, the results of all three stages are combined to
provide recommendation results without any historical information. We test the
performance of the cold start scenario by evaluating the test samples related
to the new users or items under both simulation and University Recommendation
datasets. The new users or items denote the ones that do not appear in the
training dataset but appear in the test dataset. Table 3 shows that the
performance of GAMMLI exceeds benchmark models by a great margin.

Table 3: Cold start testing accuracy comparison. Method | Simulation | University Recommendation
---|---|---
MAE | RMSE | AUC | Log Loss
GAMMLI |  1.3196$\pm$0.4966 |  1.5130$\pm$0.4029 |  0.7301$\pm$0.0185 |  0.5727$\pm$0.0647
FM | 1.9265$\pm$0.6417 | 2.3353$\pm$0.7770 | 0.6536$\pm$0.0382 | 0.6856$\pm$0.2656
DeepFM | 2.1997$\pm$0.5650 | 2.5896$\pm$0.4705 | 0.6654$\pm$0.0190 | 0.6113$\pm$0.0415
XGBoost | 1.3400$\pm$0.2624 | 1.5200$\pm$0.2568 | 0.7169$\pm$0.0173 | 0.5964$\pm$0.0133
SVD | 3.6663$\pm$0.0000 | 4.1352$\pm$0.0000 | 0.6173$\pm$0.0000 | 0.7059$\pm$0.0000

##  5 Conclusion

In this paper, a novel intrinsically explainable recommendation system called
GAMMLI is proposed based on the generalized additive model with manifest
interactions and latent interactions. It captures the main effects and
manifest interaction of observed features and discovers the latent
interactions of unobserved features. All effects can be interpreted by
visualization methods, including 1-D line plots and 2-D heatmaps. The
experiment results show that GAMMLI has advantages in both predictive
performance and interpretability. Meanwhile, the experiment results support
that GAMMLI can tackle cold start problems.

## References

  * Adomavicius and Tuzhilin, (2005) Adomavicius, G. and Tuzhilin, A. (2005).  Toward the next generation of recommender systems: A survey of the state-of-the-art and possible extensions.  IEEE Transactions on Knowledge and Data Engineering, 17(6):734–749.
  * Agarwal et al., (2020) Agarwal, R., Frosst, N., Zhang, X., Caruana, R., and Hinton, G. E. (2020).  Neural additive models: Interpretable machine learning with neural nets.  arXiv preprint arXiv:2004.13912.
  * Balabanović and Shoham, (1997) Balabanović, M. and Shoham, Y. (1997).  Fab: content-based, collaborative recommendation.  Communications of the ACM, 40(3):66–72.
  * Chen et al., (2018) Chen, C., Zhang, M., Liu, Y., and Ma, S. (2018).  Neural attentional rating regression with review-level explanations.  In Proceedings of the 2018 World Wide Web Conference, pages 1583–1592.
  * 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, pages 7–10.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: a gradient boosting machine.  Annals of Statistics, pages 1189–1232.
  * 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.  arXiv preprint arXiv:1703.04247.
  * Hastie et al., (2015) Hastie, T., Mazumder, R., Lee, J. D., and Zadeh, R. (2015).  Matrix completion and low-rank svd via fast alternating least squares.  Journal of Machine Learning Research, 16(1):3367–3402.
  * Hastie and Tibshirani, (1990) Hastie, T. J. and Tibshirani, R. J. (1990).  Generalized additive models, volume 43.  CRC press.
  * He et al., (2017) He, X., Liao, L., Zhang, H., Nie, L., Hu, X., and Chua, T.-S. (2017).  Neural collaborative filtering.  In Proceedings of the 26th International Conference on World Wide Web, pages 173–182.
  * Koren et al., (2009) Koren, Y., Bell, R., and Volinsky, C. (2009).  Matrix factorization techniques for recommender systems.  Computer, 42(8):30–37.
  * Mazumder et al., (2010) Mazumder, R., Hastie, T., and Tibshirani, R. (2010).  Spectral regularization algorithms for learning large incomplete matrices.  Journal of Machine Learning Research, 11:2287–2322.
  * Nassar et al., (2020) Nassar, N., Jafar, A., and Rahhal, Y. (2020).  A novel deep multi-criteria collaborative filtering model for recommendation system.  Knowledge-Based Systems, 187:104811.
  * Peake and Wang, (2018) Peake, G. and Wang, J. (2018).  Explanation mining: Post hoc interpretability of latent factor models for recommendation systems.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2060–2069.
  * Pujahari and Sisodia, (2020) Pujahari, A. and Sisodia, D. S. (2020).  Pair-wise preference relation based probabilistic matrix factorization for collaborative filtering in recommender system.  Knowledge-Based Systems, page 105798.
  * Rendle, (2010) Rendle, S. (2010).  Factorization machines.  In 2010 IEEE International Conference on Data Mining, pages 995–1000. IEEE.
  * Rendle, (2012) Rendle, S. (2012).  Factorization machines with libfm.  ACM Transactions on Intelligent Systems and Technology, 3(3):1–22.
  * Resnick et al., (1994) Resnick, P., Iacovou, N., Suchak, M., Bergstrom, P., and Riedl, J. (1994).  Grouplens: an open architecture for collaborative filtering of netnews.  In Proceedings of the 1994 ACM conference on Computer Supported Cooperative Work, pages 175–186.
  * Salton, (1989) Salton, G. (1989).  Automatic text processing: The transformation, analysis, and retrieval of.  Reading: Addison-Wesley, 169.
  * Sarwar et al., (2001) Sarwar, B., Karypis, G., Konstan, J., and Riedl, J. (2001).  Item-based collaborative filtering recommendation algorithms.  In Proceedings of the 10th International Conference on World Wide Web, pages 285–295.
  * Vaughan et al., (2018) Vaughan, J., Sudjianto, A., Brahimi, E., Chen, J., and Nair, V. N. (2018).  Explainable neural networks based on additive index models.  The RMA Journal, pages 40–49.
  * Wang et al., (2015) Wang, H., Wang, N., and Yeung, D.-Y. (2015).  Collaborative deep learning for recommender systems.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1235–1244.
  * Wang et al., (2018) Wang, N., Wang, H., Jia, Y., and Yin, Y. (2018).  Explainable recommendation via multi-task learning in opinionated text data.  In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 165–174.
  * Yang et al., (2020) Yang, Z., Zhang, A., and Sudjianto, A. (2020).  Enhancing explainability of neural networks through architecture constraints.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–12.
  * Yang et al., (2020) Yang, Z., Zhang, A., and Sudjianto, A. (2020).  GAMI-Net: An explainable neural network based on generalized additive models with structured interactions.  arXiv preprint arXiv:2003.07132.
  * Zhang and Lee, (2005) Zhang, D. and Lee, W. S. (2005).  A simple probabilistic approach to learning from positive and unlabeled examples.  In Proceedings of the 5th Annual UK Workshop on Computational Intelligence, pages 83–87.
  * Zhang et al., (2011) Zhang, S., Li, Q., Liu, J., and Zhou, X. J. (2011).  A novel computational framework for simultaneous integration of multiple types of genomic data to identify microrna-gene regulatory modules.  Bioinformatics, 27(13):i401–i409.
  * Zhang and Chen, (2020) Zhang, Y. and Chen, X. (2020).  Explainable recommendation: A survey and new perspectives.  Foundations and Trends® in Information Retrieval, 14(1):1–101.
  * Zhang et al., (2014) Zhang, Y., Lai, G., Zhang, M., Zhang, Y., Liu, Y., and Ma, S. (2014).  Explicit factor models for explainable recommendation based on phrase-level sentiment analysis.  In Proceedings of the 37th International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 83–92.