\section{Introduction}\label{sec:introduction}

Deep Reinforcement Learning (DRL) has been playing a significant role and achieving remarkable success in a variety of challenging problems, such as chess games \cite{alphagozero}, real-time video games \cite{alphastar, shao2018starcraft}, robotics control \cite{levine2016end}, and image classification \cite{zhao2016vehicle}. As an extension, Multi-Agent Reinforcement Learning (MARL) has also received more and more attention in many scenarios where a stand-alone agent fails in accomplishing complicated tasks due to the lack of cooperation. The existence of multiple agents poses some common issues, such as non-static environment \cite{hernandez2017survey}, partially observability \cite{maddpg, vdn}, dimension explosion \cite{fql}, credit assignment \cite{coma, qmix}, and so on. In recent research \cite{szer2004improving}, it has been demonstrated that through internal communication, agents are able to share local information and pursue the same goal, which is important to address the nonstationarity and partially observability in multi-agent environment. Of particular interest is the distinction between two lines of research, that is hand-crafted communication protocols \cite{zhang2013coordinating} and learnable communication protocols \cite{dial, commnet}. Especially the advent of MARL allows to learn the protocols in an end-to-end way. Unfortunately, communication networks in the real world have the limited bandwidth. If there are a large number of agents and they send messages excessively, the network can be easily blocked, delaying message transmission and impairing cooperative effects.

Some research in MARL has been proposed to learn communication with limited bandwidth  \cite{atoc, schednet, pruning, imac}. But some key problems are still not well studied. Existing methods focus more on the reduction of sending behaviours but pay less attention on the explicit definition of the network bandwidth. Hence their conditions under which agents decide whether to send messages are not directly designed to fulfill bandwidth limitation.

Motivated by that, this paper proposes a new Event-Triggered Communication Network (ETCNet) to realize efficient communication in MARL faced with limited bandwidth. First of all, the limited bandwidth is translated into a penalty threshold mathematically, which is further put into an optimization problem as constraints. Then the event-triggered concept is realized in the architecture, and each sending behaviour is determined by an event-triggered module with a gating policy. The open of the gating policy poses a penalty for the occupation of bandwidth, but it enhances multi-agent cooperation because of sharing information. Therefore the synthesis of gating policy is put into a constrained Markov decision process (MDP) optimization, with the multi-agent performance as the objective and the limited bandwidth as the constraint. After introducing the Lagrange multiplier, reinforcement learning adaptively finds the optimal solution in a trial-and-error manner.

To verify the effectiveness, two typical multi-agent tasks, including cooperative navigation and predator-prey, are simulated.
We compare our method with other MARL methods that also consider the limited-bandwidth constraint, including SchedNet \cite{schednet}, Gated-ACML \cite{pruning}, and Message-dropout \cite{messagedropout}. After comparison, our ETCNet is significant in reducing the bandwidth consumption and preserves the whole system with marginal impact.

\section{Related Work}\label{sec:relatedWork}

Learning communication protocols of multi-agent systems has attracted considerable attention in literature \cite{accnet, bicnet}. Existing research directions include the message content \cite{a3c2}, the robustness of communication \cite{noisy}, metrics of emergent communication \cite{measure}, the attention mechanism of learning compacted messages \cite{peng2018attention, geng2019attention, daacmp}, memory-driven communication \cite{memorydriven}, and parameter sharing \cite{parashare}. But they pay little attention to the restriction of limited bandwidth in communication network.

To efficiently utilize finite communication resources, some recent MARL methods make agents learn to choose what, when and with whom to communicate. IC3Net \cite{ic3net} extends the work of CommNet \cite{commnet} by means of Long Short-Term Memory (LSTM) and the gating mechanism.
Gated-ACML \cite{pruning} and ATOC \cite{atoc} both evaluate the importance of communication by comparing the Q-difference between sending messages and not. If the difference is
greater than a threshold, agents consider the message is valuable and choose to communicate.
SchedNet \cite{schednet} leverages weight generators to choose top-$k$ agents with apparently more valuable observations to participate in the communication group, and broadcasts their messages to the others.
The purpose of the above methods is to reduce the bandwidth consumption but there is no mathematical definition of bandwidth constraints.
IMAC \cite{imac} argues that explicit mathematical relations exist between the entropy of messages and the bandwidth, and introduces the mutual information to approximate message entropy. By restricting the mutual information to an upper bound, the problem becomes a constrained optimization that aims to learn the efficient message generators. However, the system still transmits messages at each moment, causing the waste of communication resources if the messages at consecutive moments have similar or even the totally same content.

Another drawback of above mentioned work is that agents decide whether to communicate is just based on the current observation.
Event-triggered control is an important concept in the field of control theory to reduce the update of control signals in networked control systems \cite{et2012, zhu2016event, zhang2017event}. The occupation of communication network to send signals is conditioned on the difference of a predefined energy function between the current observation and a previous one, so the bandwidth usage is restricted.
We extend the event-triggered concept to multi-agent communication and learn communication protocols to decide when sending messages.

\section{Preliminary on DEC-POMDP with Communication}\label{sec:Background}

We consider MARL in the framework of Decentralized Partially Observable Decision Process (DEC-POMDP), which is described as a tuple $\left\langle \mathcal{S}, \boldsymbol{\mathcal{A}},P,\boldsymbol{R}, \boldsymbol{\mathcal{O}}, Z, N,\gamma \right\rangle$, where $N$ is the number of agents; $\mathcal{S}$ denotes the state space of the problem;
$\boldsymbol{\mathcal{O}} = \left\{ \mathcal{O}_i \right\}_{i=1,2,...,N}$ represents the sets of observations for each agent;
$\boldsymbol{\mathcal{A}} = \left\{ \mathcal{A}_i \right\}_{i=1,2,...,N}$ denotes the sets of actions.
$Z(s,i): \mathcal{S} \rightarrow O_i$ is the observation function that determines the private observation, and
the agent $i$ receives a private observation by $o_i = Z(s,i)$.
$P(s'|s, \boldsymbol{a}) : \mathcal{S} \times \boldsymbol{\mathcal{A}} \times \mathcal{S} \rightarrow [0,1]$ represents the state transition function, where $\boldsymbol{a}=[a_1,a_2,...,a_N]$ is the joint action.
$\boldsymbol{R} = \left\{ \mathcal{R}_i \right\}_{i=1,2,...,N}: \mathcal{S} \times \boldsymbol{\mathcal{A}} \rightarrow \mathbb{R}^N$ represents the set of reward functions.
$\gamma \in [0,1]$ denotes the discount factor.
Each agent aims to learn a policy $\pi_i(a_i|o_i): \mathcal{O}_i \rightarrow \mathcal{A}_i$ that maximizes the expected discounted return $\mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t r_{i,t} \right], r_{i,t} \sim \mathcal{R}_i(s_t,\boldsymbol{a}_t)$.

Sharing observations improves the performance of the whole multi-agent system, and makes each agent learn the better policy. In such case, the policy is written by $\pi_i(a_i|o_i, \boldsymbol{o}_{-i}): \mathcal{O}_i \times \boldsymbol{\mathcal{O}}_{-i} \rightarrow \mathcal{A}_i$, where $\boldsymbol{o}_{-i} = \left[ o_{1}, ... o_{i-1}, o_{i+1}, ... o_{N}  \right] \in \boldsymbol{\mathcal{O}}_{-i}$ is the joint observation of other agents except $i$. If observations are high-dimensional, they have to be encoded to low-dimensional representations to reduce data transmission.
The policy with communication is denoted by $\pi_i(a_i|o_i, \boldsymbol{m}_{-i}): \mathcal{O}_i \times \boldsymbol{\mathcal{M}}_{-i} \rightarrow \mathcal{A}_i$, where $\boldsymbol{m}_{-i} = \left[ m_{1}, ..., m_{i-1}, m_{i+1}, ..., m_{N}  \right] \in \boldsymbol{\mathcal{M}}_{-i}$ denotes the messages that agent $i$ receives from its teammates.

In this way, we extend DEC-POMDP to a communicative one to enhance coopration. The process of communication consists of encoding, transmission, and decoding. The encoding process maps the message to a bitstream, which is transmitted through a communication channel. Decoding is the inverse operation of encoding to recover the message.
For ease of analysis, we assume the message $m$ satisfies a certain distribution $M$ and the entropy is denoted by $H(M)$. The communication network has a bandwidth $B$.
According to Shannon's source coding theorem \cite{Shannon1948}, in order to encode the message without the risk of information loss, the average number of bits $N_b$ must satisfy $ N_b \geq H(M)$.
The maximum data rate $R_{max}$ (bits per second) \cite{Freeman2004} in a noiseless channel has $R_{max} = 2B \log_2 K$, where $K$ is the number of discrete levels in the signal.
IMAC \cite{imac} combines the above two requirements together and argues that the relationship between bandwidth and message entropy has
$2B \log_2 K = R_{max} \geq n N_b \geq n H(M)$.

\section{Event-Triggered Communication Network}\label{sec:ETCNet}

Now we formally present our ETCNet. First, we detail the architecture to show its advantage of saving bandwidth and maintaining multi-agent cooperation. Second, the limited bandwidth is converted to the penalties of sending behaviours. By adding the new constraint to the multi-agent cooperative objective, the gating policy is learned by reinforcement learning to solve a constrained optimization problem.

\subsection{Architecture}\label{sec:architecture}

Figure \ref{fig:architecture} presents the architecture of ETCNet in multi-agent settings.
The execution process of each agent consists of three phases: ($i$) observation encoding; ($ii$) message gating and sending; ($iii$) message receiving and decision making.
The gating module is designed in an event-triggered way such that the sending behaviour happens only if it is necessary. At the message-receiving and decision-making stage, if an agent opens its transmission gate, other agents will receive its current message to decide cooperative actions. Otherwise, they will use the lastly received message, memorized by a zero-oreder holder (ZOH) module, to continue to cooperate.

\begin{figure}[!ht]
   \centering
   \begin{center}
   \scriptsize
     \includegraphics*[width=3.3in]{architecture.pdf}\\
   \caption{Framework of the proposed ETCNet. ENC$_i$ is the encoding module; ET$_i$ is the event-triggered gating module; ACT$_i$ is the agent policy module; ZOH represents zero-order holder.}
   \label{fig:architecture}
   \end{center}
\end{figure}

Before digging into the detailed design, we list some key notations as follows: consider $N$ agents, for agent $i$ at time $t$, its observation is denoted by $o_{i,t}$; its current message is $m_{i,t}=e_{i}\left( o_{i,t} \right)$, where $e_{i}(\cdot)$ is the encoding function; its gating action is denoted by $g_{i,t} \sim \mu_{i}(\cdot)$ , where $\mu_{i}(\cdot)$ represents the gating policy function; its action executed on the environment is denoted by $a_{i,t} \sim \pi_{i}(\cdot)$ , where $\pi_{i}(\cdot)$ represents the agent policy function.
The gating action samples from $\{0,1 \}$, where $1$ represents the event is triggered and the transmission is open. Otherwise, $0$ is not.

We specify $U_{i,t} = \left[ t^i_0, ..., t^i_r, ... \right]$ to denote the set of event-triggered time points $t^i_r$ at the current $t$, shown as Figure \ref{fig:U}.
The more explicit expression for the above mentioned variables and functions is presented.
The gating policy is denoted by $g_{i,t} \sim \mu_{i} ( m_{i,t}, m_{i,\hat{t}^i} )$, where $m_{i,\hat{t}^i}$ represents the message at the lastly triggering moment, memorized by ZOH, and $\hat{t}^i = \underset{\kappa \in U_{i,t-1}}{\arg \min } \left \{ t - \kappa \right \}$). Note that $U_{i,t-1}$ could not be updated to $U_{i,t}$ before agent $i$ makes gating decisions at $t$.
An agent chooses to send messages only when it considers the change of two inputs will facilitate the cooperation.
The agent action follows $a_{i,t} \sim \pi_{i}( o_{i,t}, \tilde{\boldsymbol{m}}_{-i,t})$, where $\tilde{\boldsymbol{m}}_{-i,t} = \left[ m_{1,\tilde{t}^1}, ..., m_{i-1,\tilde{t}^{i-1}}, m_{i+1,\tilde{t}^{i+1}}, ..., m_{N,\tilde{t}^N} \right]$, where $\tilde{t}^j= \underset{\kappa \in U_{j,t}}{\arg \min } \left \{ t - \kappa \right \}$.
In addition to own observations, the policy uses received messages (if there are) or memorized messages from others to realize cooperation.

\begin{figure}[!ht]
   \centering
   \begin{center}
   \scriptsize
     \includegraphics*[scale=0.75]{U.pdf}\\
   \caption{Illustration of $U_{i,t}$ whose elements are event-triggered time points of agent $i$.}
   \label{fig:U}
   \end{center}
\end{figure}

Compared with existing work, the biggest difference in architecture is that our ETCNet not only uses the current observation to define the sending condition, but also relies on the lastly sent message.
Beyond that, if no message is received, the agent uses the memorized message rather than the zero vector (e.g. Gated-ACML \cite{pruning}, ATOC \cite{atoc}, and SchedNet \cite{schednet}) to prevent the lose of information and preserve the cooperation performance.

\subsection{Limited-bandwidth Constraint and Penalty Threshold}\label{sec:definition}

According to the Preliminary, the maximum symbols per second on a limited-bandwidth channel satisfy
\begin{equation} \label{eq:n1}
\begin{gathered}
n  \leq \frac{2B \log_2 K}{H(M)}.
\end{gathered}
\end{equation}
However, the distribution and entropy of message $M$ is generally unknown, and all we can get are its statistic properties like mean and variance. The principle of maximum entropy  \cite{principle} proves that the Gaussian distribution has the maximum entropy compared with other probability distributions with the same mean and variance.
We can take the entropy of a Gaussian distribution $H(X)=\log (2\pi e \sigma^2), X \sim N(\mu, \sigma)$ as an upper bound of $H(M)$, where $\mu$ and $\sigma$ are the mean and variance of $M$.
Substituting it back to (\ref{eq:n1}) yields
\begin{equation} \label{eq:n2}
\begin{gathered}
n \leq \frac{2B \log_2 K}{H(X)} = \frac{4B \log_2 K}{\log (2\pi e \sigma^2)}.
\end{gathered}
\end{equation}

Suppose the gating policy has a probability of $p$ sending messages at each step, and a message has a length of $L$ symbols. The system sampling frequency is $F$. For a number of $N$ agents, the average number of symbols on the channel is equal to $N(N-1)LFp $ and should be no greater than $n$.
Then we are able to deduce an upper bound of probability that each agent is allowed to send messages at each step
\begin{equation} \label{eq:p2}
\begin{gathered}
p \leq p_{sup} = clip \left( \frac{4B \log_2 K}{\log (2\pi e \sigma^2) N(N-1)LF}, 0, 1 \right).
\end{gathered}
\end{equation}

Since sending messages or not corresponds to the open or close of the gating policy, so we can describe the occupation of bandwidth as a sum of penalties over the time horizon:
\begin{equation} \label{eq:C}
\begin{gathered}
C = \mathbb{E}\left[ \sum_{t=0}^{\infty} \gamma^t \mathbb I(g_{i,t}=1) \right]  \leq \frac{p_{sup}}{1-\gamma}=C_{sup}
\end{gathered}
\end{equation}
where $C_{sup}$ indicates the penalty threshold, and $\mathbb I(g_{i,t}=1)$ specifies the instantaneous penalty when an agent occupies the bandwidth.

With the original sum of rewards as the optimization objective, the problem now becomes solving the constrained MDP
\begin{equation} \label{eq:optimization}
\begin{split}
\max\,\, \mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t r_{i,t} \right], \quad
s.t.\quad
\mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t  g_{i,t} \right] \leq C_{sup}.
\end{split}
\end{equation}

Note that IMAC \cite{imac} also gives an explicit mathematical transformation of bandwidth limitition.
It reduces bandwidth occupation through message compression, but the transmission frequency is unchanged.
Our ECNet keeps the completeness of information and reduces the frequency of sending messages in an event-triggered way.

\subsection{Optimization Algorithm}\label{sec:algorithm}

In the implementation of ETCNet, we define three networks for each agent: EncoderNet, GatingNet, and ActorNet. They correspond to the encoding function, the gating policy function, and the agent policy function, respectively. All the homogeneous agents share the same models rather than defining different network parameters.
To lower down the learning difficulty of three networks, we separate the training into two processes. First, we train the EncoderNet and ActorNet at the full communication, that is the event-triggered module always sends message. After obtaining the well-trained  EncoderNet and ActorNet, we keep them fixed and train the GatingNet.

\textbf{Training EncoderNet and ActorNet.} The Centralized Training and Decentralized Execution (CTDE) paradigm is adopted to train EncoderNet and ActorNet to overcome the non-stationary problem.
Typically we use a centralized CriticNet parameterized by $\theta_c$ to estimate the state value function $V_{\theta_c}(\upsilon_{i}) \approx \mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t r_{i,t} \right]$ where $\upsilon_{i}=[o_{i},\boldsymbol{o}_{-i}]$ is taken as the approximation of the global state $s$. The EncoderNet and ActorNet are parameterized by $\theta_e$ and $\theta_a$, respectively.
The critic value is updated based on temporal-difference as
\begin{align}
\delta_{c} & = r_{i,t} + \gamma V_{\theta_c}(\upsilon_{i,t+1}) - V_{\theta_c}(\upsilon_{i,t}), \\
\mathcal{L}_{i,t}^{critic} & =  \delta_{c}^{2}.
\end{align}
The EncoderNet and ActorNet train the parameters by back-propagation of the policy loss
\begin{multline}
\mathcal{L}_{i,t}^{act\_enc} = -\log{\pi_{i}(a_{i,t}|o_{i,t},\boldsymbol{m}_{-i,t},\theta_a)}\delta_{c}  \\
- \alpha H(\pi_{i}(a_{i,t}|o_{i,t},\boldsymbol{m}_{-i,t},\theta_a))
\end{multline}
where
$\boldsymbol{m}_{-i,t} = \left[ m_{i,t}, ..., m_{i,t}, m_{i+1,t}, ..., m_{N,t} \right]$ for full communication and $m_{i,t} = e(o_{i,t} | \theta_e)$ . An entropy term is used to discourage premature convergence.

\textbf{Training GatingNet.}
After learning the EncoderNet and ActorNet at the pretraining stage, we apply them to ETCNet framework and keep their parameters fixed. Now we learn the GatingNet parameterized by $\theta_g$, for the gating policy $\mu_{i} ( m_{i,t}, m_{i,\hat{t}^i} |\theta_g )$ to satisfy the constrained optimization as (\ref{eq:optimization}).
We use a Lagrangian multiplier $\lambda \ge 0$  to deal with the constraint and define the Lagrangian function
\begin{equation} \label{eq:lagrangian}
\begin{split}
L(\mu_{i},\lambda)=\mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t (r_{i,t} -\lambda g_{i,t}) \right]+\lambda C_{sup}.
\end{split}
\end{equation}
The dual objective $d(\lambda)$ of the primal problem is defined as
\begin{equation} \label{eq:d}
\begin{split}
d(\lambda) = \sup_{\mu_{i}} L(\mu_{i},\lambda).
\end{split}
\end{equation}

Suppose at step $t$, we have had an estimate of $\lambda$, denoted as $\lambda_t$. The optimal solution for (\ref{eq:d}) is to find $\mu_i^*=\arg\underset{\mu_{i}}{\max}\, L(\mu_{i},\lambda_t)$,
which is in fact reduced to solve a new MDP optimization with the new reward $r'_{i,t}=r_{i,t}-\lambda_t g_{i,t}$.
Reinforcement learning is able to optimize the multi-agent performance considering the new reward signal by updating the gating policy with the learned EncoderNet and ActorNet.
A centralized LagrangianNet parameterized by $\theta_L$ is used to estimate the state value function $V_{\theta_L}(\upsilon_{i}) \approx \mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t r'_{i,t} \right]$ for the GatingNet, and the value and policy losses are defined by
\begin{align}
\delta_L & = r'_{i,t} + \gamma V_{\theta_L}(\upsilon_{i,t+1}) - V_{\theta_L}(\upsilon_{i,t}), \\
\mathcal{L}_{i,t}^{Lagr} & = \delta_L ^2, \\
\mathcal{L}_{i,t}^{gate} & = -\log{\mu_{i}(g_{i,t}|m_{i,t}, m_{i,\hat{t}^i},\theta_g)}\delta_L \nonumber \\
& \qquad \qquad \quad - \alpha H(\mu_{i}(g_{i,t}|m_{i,t}, m_{i,\hat{t}^i},\theta_L)).
\end{align}
Note that $\lambda_t$ is an estimate of the true $\lambda$ and the optimal multiplier $\lambda^*$ satisfies
\begin{equation} \label{eq:lambda_star}
\begin{split}
\lambda^* = \arg\underset{\lambda \ge 0}{\min}\, g(\lambda).
\end{split}
\end{equation}
Then the $\lambda_t$ is updated following
\begin{align}
\lambda_{t+1} & = (\lambda_{t}-\eta_{\lambda} \nabla d(\lambda))^+ , \\
\nabla d(\lambda) & =\frac{\partial d(\lambda)}{\partial \lambda}=\left.\frac{\partial L\left(\mu_{i}, \lambda\right)}{\partial \lambda}\right|_{\mu_{i}=\mu_i^*} \nonumber \\
&=-\mathbb{E}\left[\sum_{t=0}^{\infty} \gamma^{t} g_{i, t}\right]+C_{sup}.
\end{align}
A PenaltyNet parameterized by $\theta_p$ is used to estimate the penalty value function to approximate $V_{\theta_p}(\upsilon_{i}) \approx \mathbb{E}\left[ \sum_{t=0}^\infty \gamma^t g_{i,t} \right]$. Its parameters are updated based on temporal-difference
\begin{equation} \label{eq:penalty_g}
\begin{gathered}
\mathcal{L}_{i,t}^{penalty} = \left[ g_{i,t} + \gamma V_{\theta_p}(\upsilon_{i,t+1}) - V_{\theta_p}(\upsilon_{i,t}) \right]^2.
\end{gathered}
\end{equation}
Then the update of $\lambda$ becomes
\begin{equation}
\begin{split}
\lambda_{t+1} = (\lambda_{t}-\eta_{\lambda} (-V_{\theta_p}+C_{sup}))^+.
\end{split}
\end{equation}

Considering the variance of messages varies with the change of the gating policy, we calculate the variance and update the penalty threshold $C_{sup}$ periodically throughout the training. At the end of gradient iterations, the optimal policy $\mu_{i}^*$ for the unconstrained problem is obtained. The pseudocode of training ETCNet is presented in the supplementary material.
Note that in the field of control theory, event-triggered control mainly forcus on reducing the update of control signals, and its triggering condition relies on a predefined energy function \cite{et2012, zhu2016event, zhang2017event}. Here the event-triggered module is applied to reduce the transmission of messages, and the communication protocols are learned from scratch.

\section {Experiments}\label{sec:Experiment}

Two variants of multi-agent partical environments are introduced to test the performance of ETCNet on saving communication resurces, that is Cooperative Navigation and Predator and Prey \cite{maddpg} shown in Figure \ref{fig:tasks}.
The detailed experimental configurations are expatiated in the following subsections as well as in the supplementary material. \\
\textbf{Cooperative Navigation.} There are two agents and each agent aims to arrive at a specified and dynamic destination by moving along discrete directions. Each agent only observes the destination and position of the other agent. Once an agent reaches its destination, it will get a positive reward. Otherwise there is always a negative reward until the end of episode. \\
\textbf{Predator and Prey.} In this scenario, $n$ predators chase $m$ preys within a certain area. The top and the bottom, and the left and the right of the area are interconnected.
Predators and preys have the same velocity, and preys are equipped with a fixed escape policy (running from the closest predator) with a complete map vision.
A predator only has a local view around itself, so they have to cooperate to capture preys and are required to avoid collision with other predators.

\begin{figure}[h]
\begin{center}
\subfigure[Cooperative Navigation]{
\includegraphics[width=0.4\linewidth]{ng.pdf}
}
\subfigure[Predator and Prey]{
\includegraphics[width=0.4\linewidth]{pp.pdf}
}
\caption{Illustrations of two multi-agent tasks. (a) A circle represents an agent and the pentagram with the same color represents its destination. (b) A green circle represents a prey and a small square represents a predator with the local view surrounded by a large square. }
\label{fig:tasks}
\end{center}
\end{figure}

Our baselines for limited bandwidth are (1) Gated-ACML \cite{pruning}, (2) SchedNet \cite{schednet}, and (3) Message-Dropout \cite{noisy}. In addition, A3C2 \cite{a3c2} is introduced as a full-communication version. In fairness, the communication protocols of all methods are set to the same as A3C2, that is, each agent sends the same encoded message to the others.
Ideally, A3C2 allows agents to fully communicate and should have the best performance. In Message-dropout, agents are allowed to send messages only with a certain probability, so it can be seen as the randomly-failed-communication version of A3C2. The gate module in Gated-ACML works much the same as the attention module in ATOC \cite{atoc}, so we choose Gated-ACML as the baseline to represent the class of $\Delta Q$-based algorithms.
Note that SchedNet selects top-$k$ agents to send messages at each step, leading to the discrete property of sending probability in the $N$-agent system, such as $1/N$ and $k/N$.
For the sake of fairness, we set identically desired sending probabilities for every baseline by adjusting some parameters, such as bandwidth of ETCNet, Q-difference threshold of Gated-ACML, and the probability of dropout in Message-Dropout.
Other configurations are given in detail in the supplementary material.

\subsection{Cooperative Navigation}\label{sec:ng}

In this task, we consider communication network with $K=2$ and $L=6$. The sampling frequency of the system has $F=45$Hz. The bandwidth is first limited to 170 bit/s. At the full communication of ETCNet, the variance of massages is $\sigma^2=0.69$, so the maximally allowed communication probability is about $50\%$. After calculating the penalty threshold and continuing the GatingNet training in the event-triggered architechture, we observe that the messsage variance varies slightly smaller to $\sigma^2=0.57$, which further relaxes the upper bound of communication probability. In fact, the final ETCNet agents send messages at each step only with 46\% probability, lower than the desired 50\%.

We take the number of steps accomplishing the task as the evaluation. The fewer steps, the better performance. Table \ref{tab:ng1} shows the results of ETCNet and baselines. It is observed that ETCNet is far superior to other methods under the same communication constraint, and is closest to the performance of full communication.
\begin{table}[h]
  \centering
  \caption{The performance of ETCNet and baselines for Cooperative Navigation with $50\%$ communication.}
    \begin{tabular}{|c|c|}
    \hline \hline
    {\bf Methods} & {\bf Steps}  \\ \hline
    \bf{ETCNet} & $\bf{16.50 \pm 6.03}$ \\
    Gated-ACML & $33.83 \pm 10.08$  \\
    SchedNet & $39.49 \pm 3.27$ \\
    Message-Dropout & $35.09 \pm 7.96$ \\ \hline
    A3C2(Full communication) & $16.20 \pm 5.85$  \\
    \hline \hline
    \end{tabular}%
  \label{tab:ng1}%
\end{table}%

To demonstrate that ETCNet can greatly reduce bandwidth consumption and preserve the multi-agent cooperation, we repeat the experiment with $B = 100$ bit/s and $B =60$ bit/s.
Figure \ref{fig:nglearningcurves} shows the learning curves of ETCNet with different bandwidths.
The top subgraph shows that the more limitation on bandwidth, the more degradation of performance in the early stage. But the curves are still stabilized back to near optimality through the later training. The bottom subgraph shows that all experiments optimize the sending penalties to satisfy corresponding thresholds. The final communication percentages are $14\%$, $24\%$ and $46\%$ for $B$ equal to 60, 100, and 170, respectively. It is concluded that ETCNet is capable of adjusting to different bandwidth constraints and preserving the best performance.
\begin{figure}[h]
\begin{center}
\subfigure[Step evaluation]{
\includegraphics[width=3in, height=1.2in]{ng_step.png}
}
\subfigure[Mean penalty per step]{
\includegraphics[width=3in, height=1.2in]{ng_penalty.png}
}
\caption{Learning curves for Cooperative Navigation with respect to different bandwidths.}
\label{fig:nglearningcurves}
\end{center}
\end{figure}

\subsection{Predator and Prey}\label{sec:pp}

In this task, we first consider 3-agent Predator and Prey. The communication networks has $K=2$ and $L=15$. The sampling frequency of the system has $F=45$Hz. The bandwidth is first limited to 580 bit/s. At the full communication of ETCNet, the variance of massages is $\sigma^2=0.330$, so the maximally allowed communication probability is about $33.3\%$. After calculating the penalty threshold and continuing the GatingNet training in the event-triggered architecture, we observe that the message variance barely changes, so the bandwidth limitation is not violated.

We take the number of steps accomplishing the task as the evaluation. The fewer steps, the better performance. Table \ref{tab:pp1} gives the performance of ETCNet and baselines with 3 agents and $33.3\%$ desired communication probability. It shows that SchedNet is competitive to our ETCNet and they both achieves the similar performance to the full-communication results.
We further compare their performance under different communication probabilities and agent numbers. The results are also listed in Table \ref{tab:pp1}. In some experiments, ETCNet outperforms all baselines, while in the others, it is competitive to the best SchedNet. The performance gap between ETCNet and the full communication is quite small. It is worth noting that ETCNet works in a variety of limited-bandwidth settings and optimally exploits the bandwidth. The communication probability in SchedNet is proportional to
$1/N$ since its mechanism is to select top-$k$ agents to send messages at each step.
\begin{table*}[t]
  \centering
  \caption{The performance of ETCNet and baselines for Predator and Prey with different agent numbers and communication upper bound.}
    \begin{tabular}{|c|c|c|c|c|c|}
    \hline \hline
      \multirow{2}{*} & {2 agents} & \multicolumn{2}{|c|}{3 agents} & \multicolumn{2}{|c|}{4 agents}  \\  \cline{2-6}
     ~ & {$\le $ 50\%} & {$\le $ 33.3\%} & {$\le $ 66.6\%} & {$\le $ 25\%} & {$\le $ 50\%}  \\ \hline
    \bf{ETCNet} & $\bf{52.19 \pm 14.63}$   & $\bf{54.59 \pm 18.63}$ & $53.52 \pm 17.81$ & $47.79 \pm 19.02$ & $\bf{46.04 \pm 17.39}$ \\
    {Gated-ACML} & $93.80 \pm 13.33$ & $65.87 \pm 22.05$ & $60.35 \pm 20.40$ & $70.67 \pm 28.86$ & $46.28 \pm 18.43$  \\
    {SchedNet} & $54.06 \pm 13.14$ & $54.95 \pm 17.41$ & $\bf{52.78 \pm 17.71}$ & $\bf{45.52 \pm 16.69}$ & $46.66 \pm 18.05$  \\
    {Message-Dropout} & $91.385 \pm 13.94$ & $85.33 \pm 19.32$ & $61.01 \pm 19.53$ & $73.933 \pm 26.404$ & $59.63 \pm 23.98$ \\ \hline
    {A3C2 (Full communication)} & {$51.08 \pm 12.68$} & \multicolumn{2}{|c|}{$50.68 \pm 16.75$} & \multicolumn{2}{|c|}{$45.256 \pm 16.755$}  \\
    \hline \hline
    \end{tabular}%
  \label{tab:pp1}%
\end{table*}%

Figure \ref{fig:pplearningcurves} shows the learning curves of ETCNet for 3-agent Predator and Prey with different bandwidths ($B = 580$ bit/s and $B = 1200$ bit/s). The learned gating policies send messages at probabilities $p=31.7\% $ and $p = 55.7\%$, respectively.
The plot shows that the learning process with the lower bandwidth has the lower frequency of sending messages, but the evaluation is worse than the learner with higher bandwidth. It is consistent with the fact that more communication is benificial to multi-agent cooperation.
\begin{figure}[h]
\begin{center}
\subfigure[Step evaluation]{
\includegraphics[width=2.9in, height=1.2in]{pp_step.png}
}
\subfigure[Mean penalty per step]{
\includegraphics[width=3in, height=1.2in]{pp_penalty.png}
}
\caption{Learning curves for Predator and Prey with respect to different bandwidths.}
\label{fig:pplearningcurves}
\end{center}
\end{figure}

\subsection{Event-triggered Gate}\label{sec:ET}
We argue that agents in ETCNet send messages only when necessary. We demonstrate this argument by analyzing system trajectories of Cooperative Navigation obtained by ETCNet.
Figure \ref{fig:time} shows gating actions and representative sceneries in an epoch.
We first focus on the sending behaviours of blue agent. It sends message at the starting point (a) and does not send at (b) because of no changes in its observations. It even refuses to send message at (c) when the red destination moves. It is because the lastly  received message of the red agent can still help in choosing the correct action (towards the left), especially considering that the red destination is likely to change later. It sends message at (d) because the red agent will be mislead to the wrong direction by the old message.
Now let us see the gating of red agent. It does not send message at (e) although its observation changes with the blue destination. The blue agent continues to utilize the old message and moves forwards the correct direction. When the blue agent reaches its destination at (f), the epoch terminates with both agents accomplishing their tasks. It suggests that ETCNet agents trigger
the gating policy only when the communication is important for cooperation, not simply determined by the change of observation.

\begin{figure*}[t]
   \centering
   \begin{center}
   \scriptsize
     \includegraphics*[width=6.5in]{time.pdf}\\
   \caption{An event-triggered gating display in an epoch. The circle represents an agent and the pentagram with the same color represents its destination. The yellow ring surrounding an agent indicates it is currently sending message. }
   \label{fig:time}
   \end{center}
\end{figure*}

\subsection{Ablation}\label{sec:Ablation}

To investigate the effect of the memorized messages in agent policy $\pi_i$ and gating policy $\mu_i$ , we conduct some ablation studies.
Figure \ref{fig:Ablation} shows the learning curves of ETCNet on Cooperative Navigation with different ablation.
\begin{figure}[htbp]
\begin{center}
\subfigure[Step evaluation]{
\includegraphics[width=2.94in, height=1.3in]{ab_step.pdf}
}
\subfigure[Mean penalty per step]{
\includegraphics[width=3in, height=1.3in]{ab_penalty.pdf}
}
\caption{Learning curves with respect to the ablation of memorized messages in $\pi_i$ and $\mu_i$.}
\label{fig:Ablation}
\end{center}
\end{figure}

First, we analyze the effect of the memorized messages in agent policy $\pi_i$, which is used to assist decision-making. In ablation, we pad the zero vector to ActorNet when an agent receives nothing.
We observe that the system tends to learn full communication because the failure of storing the lastly received messages greatly degrades the cooperation.

Next, we disentangle the influence of the memorized messages in gating policy $\mu_i$. The GatingNet only takes the current message as input, regardless of the lastly triggering message.
The two blue lines reveal that the ablation of the memorized messages in $\mu_i$ leads to performance degradation and bandwidth consumption.
Without knowing what has been sent in the past, the event-triggered learner has to increase sending frequency to ensure valuable messages are successfully received by others. It disturbs the learning of the multi-agent policy and deteriorates cooperation effects.

\section {Conclusion}\label{sec:Conclusion}

In this work, we propose a novel method, Event-Triggered Communication Network, for multi-agent reinforcement learning with limited-bandwidth communication.
Bandwidth limitation is mathematically transformed into a penalty threshold to restrict sending behaviours. Combined with the multi-agent optimization objective,
we establish a constrained MDP model and learns the event-triggered communication protocols. Through our experiments, we show that
ETCNet learns to send messages only when necessary while performing well in cooperation under different bandwidth constraints.

\section*{Acknowledgments}

This work was supported in part by the National Key Research and Development Program of China under Grant 2018AAA0101005 and Grant 2018AAA0102404 and in part by Huawei Noah's Ark Lab under Grant YBN2020075035.


\section{Introduction}\label{intro}
Working with time-series data is a perennial challenge in science.
In medicine, classification of electrocardiogram allows detecting anomalies associated with different illnesses, \cite{Jambukia2015}.
In ornithology, scientists can classify audio signals of recorded bird sounds for species recognition, \cite{Qian2015}.
These are, however, cases of univariate time series.
The multivariate case remains specially difficult, as it is seen on the established benchmark in the time-series community by \cite{Bagnall2017}.
Compared to univariate time-series, methods for multivariate data remain uncommon.
Similarly, their performance is often underwhelming.

Entries, such as \cite{M2016}, \cite{Nielson2015}, \cite{Anderson2018}, \cite{Rivera-Castro2019-nv}, validate the significance of topological data analysis (TDA) for time series data analysis.
TDA is a set of statistical methods aimed to analyse complex datasets.
It relies on ideas from topology and computational geometry, \cite{Edelsbrunner2002}.
Especially useful in the case of high-dimensional, complex datasets, as it allows to identify structure and shape.
This makes TDA robust to small perturbations in data, \cite{Cohen-Steiner2007}.
Similarly, TDA can extract seasonality and periodicity, \cite{Perea2015}, \cite{Dotko2019}.
Authors have found success with TDA for a wide array of applications ranging from marketing data, \cite{Rivera-Castro2019}, to cryptocurrencies, \cite{Rivera-Castro2019-lv}.
This highlights one of its benefits.
It does not require particular assumptions on the data.
These characteristics make TDA-based models interesting to address the classification of time series.

\paragraph{Innovation}
This work features several innovations compared to other works, such as \cite{Gidea2018}, \cite{Seversky2016}, \cite{Wu2019}, \cite{Pereira2015}, \cite{Umeda}.
Firstly, we compare our approach not only with baseline models, such as 1-Nearest Neighbour Euclidean or Dynamic Time Wrapping models.
Rather than that, we also compare TOTOPO with other TDA-based approaches.
These techniques are either widely discussed in the literature or are current state-of-the-art models for time series classification.
Secondly, common approaches for application of TDA to time series usually rely on creating point clouds from time series data using sliding windows.
Our approach, on the other hand, additionally includes direct persistence diagrams.
They enable us to enrich the model with additional relevant features.
Finally, we combine models based on different TDA descriptors as an ensemble.
This allows us to enrich TOTOPO even further.
With this, we make it more flexible for different types of signals.

\paragraph{Approach}
In general, our proposed approach consists of four steps.
First, we create Persistence Diagrams from time series.
We do this both directly and using sliding windows.
Second, we calculate TDA inputs from the PDs.
Third, we train base learners on these TDA inputs as well as the original time series.
Finally, we generate an ensemble of base learners.

\section{Topological ideas behind}
To create the topological inputs, first, we extract persistence diagrams from the original signals in two ways.
We do a direct extraction.
Also, we do a transformation of the time series with a sliding window.
\cite{Cohen-Steiner2007} described how to do a direct extraction of persistence diagrams from signals.
Such persistence diagrams contain information about extreme points and thus shape of time series.
Our work is the first in the time-series classification literature to use summaries from these diagrams.
Further, we transform time series into point clouds using sliding windows embedding and extract PDs from them.
This is akin to the works of \cite{Seversky2016} and \cite{Pereira2015}.

\section{TOTOPO}\label{feats_pds}
In order to include TDA into machine learning pipeline, we create TDA Inputs, which summarize relevant information from PDs.
In this work, three main sets of TDA Inputs are calculated: Betti series, $L^2$-norms series and TDA Summaries.

\paragraph{Feature extraction}
For each Persistence Diagram ($PD$), we calculate features separately for the set of 0-dim and 1-dim persistent holes.
We describe the following five features.
First, the number of holes of each dimension $d$.
It is the number of points in the $PD_d$, we define them as $N_d$.
Second, the maximal lifetime of holes of each dimension.
Lifetime is the difference between death and birth, $\max_{d} = \max_{z_i \in PD_d}(d_i - b_i)$.
Third, the number of "relevant" holes of each dimension. This is the number of points with a lifetime higher or equal than $ratio * \max_d$.
$ratio \in (0, 1)$ is the parameter describing the strength of relevancy.
Fourth, the average lifetime of holes of each dimension $mean_d = \frac{1}{N_d} \sum_{i=0}^{N_d}(d_i-b_i)$.
Fifth, the sum of lifetimes of all holes of each dimension $sum_d =  \sum_{i=0}^{N_d}(d_i-b_i)$.
Since direct PDs contain only $0$-dimensional holes, we extract only one set of five features from them.
Overall, TDA Summaries is a set of 15 numerical values.
They summarize information from persistence diagrams.

\paragraph{Betti series}
Another topological input that we use is Betti series, \cite{Umeda}.
For each dimension, Betti series is a sequence of Betti numbers $B_1, B_2, \ldots, B_k$.
$k$ is an arbitrary number.
In our case, it is set to 100.
Each Betti number shows how many persistent holes are "alive" at a particular radius.
Thus, Betti series is a sequence of Betti numbers for each radius.

\paragraph{L2 norms}
\cite{Gidea2018} and \cite{Gidea2018a} discuss the advantages of $L^2$-norms for time series analysis.
They serve as an inspiration for this work.
Firstly, we extract point clouds $P$ from each time series, using sliding windows of size $d$.
Further, we apply a larger sliding window of size $W$ to vectors from $P$.
It creates a series of small point clouds of size $d \times W$.
Then for each of these point clouds, we calculate the persistence diagrams and persistence landscapes.
Finally, we extract a series of $L^2$ norms from consequent persistence landscapes.
We can define it as $\|\lambda\|_2 = \bigg( \sum_{k=1}^{\infty} \|\lambda_k\|_2^2 \bigg)^{1/2}$.
Here, $\lambda_k$ are lines in a persistence landscape.

\paragraph{Base learner}
We summarize the architecture of the base learner in \autoref{base_learner}.
The architecture can process multi-channel input signals or features.
It contains three convolutional blocks, global averaging, and fully connected layer with SoftMax activation in the end.
As an activation function, we use LeakyReLu.
Also, we consider dropout for training at the end of two blocks out of three, to avoid overfitting.
Further, we train the base learners on four sets of inputs: TS — original time series, Betti series, $L^2$-norms, and TDA Summaries.

\paragraph{Ensembling}
The approaches of \cite{Bagnall2017} and \cite{Lines2018} inspire our ensemble.
When we train the base learners, we save their training losses for each dataset.
Then we weight the predictions of each model by its rank among other models.
We sum all weighted predictions to obtain the final one.
The final prediction is thus the class with a maximal weighted probability.
We outline the ensemble in Algorithm 1.

\begin{algorithm}[ht]
\caption{Ensemble}
\begin{algorithmic}
\State $train\_losses = $ train losses of each model
\State $predictions = $ predictions of each model
\State $model\_names = \{ts,\:Betti,\: TDAFeats, \:L2norms\}$
\State $number\_models = k$
\State $inputs=\{\textbf{ts},\:\textbf{Betti},\:\textbf{TDAFeats}, \:\textbf{L2norms}\}$
	\State Sort ascending $train\_losses$ $\to$ $models\_sorted$
	\State Assign votes: $k, k-1, \ldots, 1$ to best, second best, $\ldots$, worst models.
    \State $votes = \{models\_sorted_1:k, \: models\_sorted_2:k-1,\: \ldots,  \: models\_sorted_i:1\}$
	        \State $combinated\_prediction = 0$
	        \For{$key$, $vote$ in $votes$}
	            \State $combinated\_prediction \: += predictions\{key\} \times vote$
	        \EndFor
    \State $predicted\_classes = \{argmax(combinated\_prediction_i) \textnormal{ for i in range}(N)\}$
    \State \textbf{return} $predicted\_classes$
\end{algorithmic}
\end{algorithm}

\section{Discussion of results}
We compare TOTOPO against four techniques widely used in the time series community and other TDA-based time series classifiers.
From \cite{Bagnall2017}, we consider 78 datasets for the comparison.
To evaluate the overall performance of algorithms, we use the critical difference (CD) diagram.
We present it in \autoref{cd_diagram}.
TOTOPO is the third best classifier in the 78 datasets.
It is statistically similar to the state-of-the-art method BOSS.
TOTOPO beats BOSS on 44\% of the datasets.

\begin{figure}[ht]
\centering
\includegraphics[width=0.9\linewidth]{images/cd.pdf}
\caption{CD diagram with black horizontal lines showing statistically similar classifiers.}
\label{cd_diagram}
\end{figure}

\paragraph{Multivariate time series}
Multivariate time series methods have less entries in the literature than their univariate counterparts.
Most methods work either for univariate or for multivariate time series.
This is true for the state of the art models for univariate data, such as HIVE-COTE.
TOTOPO has an outstanding characteristic.
For multivariate time series, TOTOPO belongs to the state-of-the-art in terms of performance.
For this evaluationn, we use the datasets available in \cite{Bagnall2018}.
We can see the results in a comparison against two baselines and a novel SOTA deep-learning architecture by Franceschi, \cite{Franceschi}.
\autoref{tab:multivariate_ranking} presents the average ranks for each model.
Lower values are better.
Our model is the second best performer in the comparison.

\begin{table}[H]
    \centering
    \caption{Multivariate models, ordered by average rank. Less is better.}
    \begin{tabular}{cccc}
    Method & Average Rank & Method & Average Rank \\ \hline
    Franceschi & 1.79 & 1-NN DTW & 2.58\\
    TOTOPO & 2.5 & 1-NN Euclidean & 3.11\\
    \end{tabular}

    \label{tab:multivariate_ranking}
\end{table}

\paragraph{Best results}
We want to highlight the dataset 4-channel ERing.
The data is related to Human Activity Recognition.
On this dataset, our approach significantly outperforms all other methods.
It has an accuracy of 94.4\% as opposed to an average accuracy of 13.3\% in the comparison.
Similarly, TOTOPO excels on the univariate datasets of type DEVICE.
We show in \autoref{univariate_res} the improvement over Dynamic Time Warping.
Our approach outperforms all methods.
This is the case for both baselines and state-of-the-art, HIVE-COTE and BOSS, on five out of six datasets.
Time series from this source type correspond to different devices' electricity consumption.
These time series have regions of periodic behaviour.
We assume that TDA makes use of this characteristic.

\paragraph{Limitations}\label{limitations}
TOTOPO underperforms with datasets containing similar classes.
This is the case, for example, for the OliveOil dataset.
The differences between classes among the time series are extremely small.
Our explanation is that TDA is invariant to constant shifts and small perturbations in data.
It requires larger differences in the classes to generate different homologies.

\setlength{\tabcolsep}{2.5pt}
\begin{longtable}{ll|cccccc|c}

\captionsetup{width=15cm}
\caption{Accuracy improvement over DTW across classifiers. Results in \%. More is better.}
\label{univariate_res}\\
\toprule
        & & Seversky & Wu 30 & Pereira & Umeda & HC & BOSS & TOTOPO \\
        Type & Dataset & & & & &  \\
\midrule
\endhead
\midrule
\multicolumn{7}{r}{{Continued on next page}} \\
\midrule
\endfoot
\bottomrule
\endlastfoot
\multirow{6}{*}{\rotatebox[origin=c]{90}{DEVICE}}
& LargeKitchenAppliances  & -6,93    & -46,14 & 2,94    & -9,33    & 6.93  & -2.94  & \textbf{9.33}\\
& SmallKitchenAppliances  & -21,06  & -50,94  & -8,26    & -21,33  & 21.06  & 8.26  & \textbf{21.33}\\
& ScreenType        & -19,2    & -6,4  & -6,67    & -23,2    & 19.2  & 6.67  & \textbf{23.2}\\
& RefrigerationDevices    & -9,33    & -10,67  & -3,47    & -14,13  & 9.33  & 3.47  & \textbf{14.13}\\
& ElectricDevices      & -16,82  & N/A  & -19,71  & -16,29  & 16.82  & \textbf{19.71}  & 16.29\\
& Computers          & -6    & -20  & -5,59    & -12,8    & 6    & 5.6  & \textbf{12.8}\\
\end{longtable}

\paragraph{Noise robustness exploration}\label{noise_robustness}
We want to understand the robustness of TOTOPO.
For this, we apply three noise levels to the signals.
For each noise level we use a targeted Signal to Noise Ration (SNR).
We express it in decibels: 10dB, 15dB and 20dB.
For each time series $X$, we conver SNR into a linear one.
We follow the suggestions in \cite{Sherman2007}.
For the comparison, we plot average accuracies for each noise level on \autoref{noise}.
On the $x$-axis, we show noise levels.
"0" corresponds to the original time series.
"1" is equivalent to an SNR of 20dB, "2" for 15dB and "3" for 10dB.
For small noise levels, TOTOPO outperforms all other models.
Its accuracy deteriorates less.
Overall, this suggests that TDA is indeed robust to noise, especially in case of small changes.

\begin{figure}
\centering
\begin{minipage}{.49\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/base_learner.pdf}
  \captionof{figure}{Overview of the structure of the base learner classifier}
  \label{base_learner}
\end{minipage}%
\hfill
\begin{minipage}{.49\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/noise.pdf}
  \captionof{figure}{Comparison of different algorithms' performances for different noise levels}
  \label{noise}
\end{minipage}
\end{figure}

\section{Outlook}
We aim to improve the base learner classifiers in future works.
This includes the hyperparameter tuning for each of the datasets.
There is also space to optimize the architecture.
Additionally, we seek to evaluate the addition of non-topological descriptors into the model.
With this, we want to overcome the methodological problems of TOTOPO described in \autoref{limitations}.
These can be features to account for constant shifts in data.
We can also use other common descriptors for time series data, i.e., autocorrelation, Fourier coefficients, etc.
Moreover, we aim to create and evaluate new topological descriptors.
They will allow us to improve TOTOPO even further.

\newpage
\bibliographystyle{abbrv}


\section{Introduction}
Deep neural networks (DNNs) have exhibited significant success in many tasks, and the interpretability of DNNs has received increasing attention in recent years.
Most previous studies of post-hoc explanation of DNNs either explain DNN semantically/visually~\cite{lundberg2017unified, ribeiro2016should},
or analyze the representation capacity of DNNs~\cite{higgins2017beta, achille2018emergence, fort2019stiffness, Liang2020Knowledge}.

In this paper, we propose a new perspective to explain a trained DNN, \emph{i.e.} quantifying interactions among input variables that are used by the DNN during the inference process. Each input variable of a DNN usually does not work individually. Instead, input variables may cooperate with other variables to make inferences. We can consider the strongly interacted input variables to form a prototype feature (or a coalition), which is memorized by the DNN.
For example, the face is a prototype feature for person detection, which is comprised of the eyes, nose, and mouth. Each compositional part inside a face does not contribute to the inference individually, but they collaborate with each other to make the combination of all parts a meaningful feature for person detection.

In fact, our research group led by Dr. Q. Zhang have used game-theoretic interactions as new metrics to explain signal-processing behaviors in trained DNNs. Specifically, we have proved that the mathematic connections between interactions and network output~\cite{zhang2020game}, and we have used interactions to explain adversarial samples~\cite{unified2020wang} and the generalization power of DNNs~\cite{zhang2020interpreting}.

Previous studies mainly focused on the interaction between two variables~\cite{lundberg2018consistent, singh2018hierarchical, murdoch2018beyond, janizek2020explaining}.
~\citet{Grabisch1999AnAA} measured $2^n$ different interaction values for all the $2^n$ combinations of $n$ variables, and the cost of computing each interaction value is NP-hard.
In comparison to the unaffordable computational cost, our research summarizes all $2^n$ interactions into a single metric, which provides a global view to explain potential prototype features encoded by DNNs.

In this study, we define the significance of interactions among multiple input variables based on game theory.
In the game theory, each input variable can be viewed as a player. All players (variables) are supposed to obtain a high reward in a game. Let us consider interactions \emph{w.r.t.} the scalar output $y=f(I)$ of the DNN (or interactions \emph{w.r.t.} one dimension of the vectorized network output), given the input $I$ with overall $n$ players. The absolute change of $y$ \emph{w.r.t.} an empty input $\bf 0$, \emph{i.e.} $|f(I)-f(\bf 0)|$, can be used as the total reward of all players (input variables). The reward can be allocated to each player, which indicates the contribution made by the player (input variable), denoted by $\phi_1, \phi_2,\cdots,\phi_n$, satisfies $|f(I)-f({\bf 0})|=\sum_{i=1}^n\phi_i$.

A simple definition of the interaction is given as follows. If a set of $m$ players $S$ always participates in the game together, these players can be regarded to form a coalition. The coalition can obtain a reward, denoted by $\phi_S$. $\phi_S$ is usually different from the sum of rewards when players in the set $S$ participate in the game individually. The additional reward $\phi_S-\sum_{i\in S}\phi_i$ obtained by the coalition can be quantified as the interaction.
If $\phi_S-\sum_{i\in S}\phi_i>0$, we consider players in the coalition have a positive effect. Whereas, $\phi_S-\sum_{i\in S}\phi_i<0$ indicates an negative/adversarial effect among the set of variables in $S$.

However, $\phi_S-\sum_{i\in S}\phi_i$ mainly measures the interaction in a single coalition, whose interaction is either purely positive or purely negative. In real applications, for the set $S$ with $m$ players, players may form at most $2^m-m-1$ different coalitions.
Some coalitions have positive interaction effects, while others have negative interaction effects. Thus, interactions of different coalitions can counteract each other.

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\linewidth, height=0.12\linewidth]{./fig/overview6.pdf}
\caption{Overview of the proposed method. In this paper, we aim to quantify the significance of interactions among a set of input variables. We divide the set of input variables into various coalitions, which mainly encode either positive or negative interaction effects. In this figure, each red region represents a coalition that mainly encodes positive interaction effects, and each blue region represents a coalition that mainly encodes negative interaction effects. Thickness of the edge indicates the strength of the interaction.}
\label{fig:overview}
\end{figure*}

In order to objectively measure the significance of the overall interactions within a specific set $S$ with $m$ players (variables), we propose a new metric, which reflects both positive and negative interaction effects.
We develop a method to divide the $m$ players into a few coalitions, and ensures that each coalition mainly contains positive interaction effects. Similarly, we can also divide the $m$ players into another few coalitions, each of which mainly encodes negative interaction effects.
In this way, we can quantify the significance of both positive and negative interaction effects, as shown in Figure~\ref{fig:overview}.

In experiments, we have applied our method to five DNNs with different architectures for various tasks. Experimental results have provided new insights to understand these DNNs. Besides, our method can be used to mine prototype features towards correct and incorrect predictions of DNNs.

Contributions of this study can be summarized as follows. (1) In this paper, we define and quantify the significance of interactions among multiple input variables in the DNN, which can reflect both positive and negative interaction effects. (2) Our method can extract prototype features, which provides new perspectives to understand the DNN, and mines prototype features towards correct and incorrect predictions of the DNN. (3) We develop a method to efficiently approximate the significance of interactions.

\section{Related work}
\subsection{Interpretability of DNNs in NLP}
The interpretability of DNNs is an emerging direction in machine learning, and many methods have been proposed to explain the DNN in natural language processing (NLP). Some methods computed importance/attribution/saliency values for input variables \emph{w.r.t.} the output of the DNN~\cite{lundberg2017unified, ribeiro2016should, binder2016layerwise, simonyan2013deep, shrikumar2016not, springenberg2014striving, li2015visualizing, li2016understanding}.
Another kind of methods to understand DNNs are to measure the representation capacity of DNNs~\cite{higgins2017beta, achille2018emergence, fort2019stiffness, Liang2020Knowledge, cheng2020explaining}. \citet{guan2019towards} proposed a metric to measure how much information of an input variable was encoded in an intermediate-layer of the DNN.
Some studies designed a network to learn interpretable feature representations.~\citet{chung2017hierarchical} modified the architecture of recurrent neural network (RNN) to capture the latent hierarchical structure in the sentence.~\citet{shen2019ordered} revised the LSTM to learn interpretable representations.
Other studies explained the feature processing encoded in a DNN by embedding hierarchical structures into the DNN~\cite{tai2015improved-0, wang2019tree}.~\citet{wang2019self-attention} used structural position representations to model the latent structure of the input sentence by embedding a dependency tree into the DNN.~\citet{dyer2016recurrent-0} proposed a generative model to model hierarchical relationships among words and phrases.

In comparison, this paper focuses on the significance of interactions among multiple input variables, which provides a new perspective to understand the behavior of DNNs.

\subsection{Interactions}
Several studies explored the interaction between input variables.~\citet{bien2013a} developed an algorithm to learn hierarchical pairwise interactions inside an additive model.~\citet{sorokina2008detecting} proposed a method to detect the statistical interaction via an additive model-based ensemble of regression trees.
Some studies defined different types of interactions between variables to explain DNNs from different perspectives.~\citet{murdoch2018beyond} proposed to use the contextual decomposition to extract the interaction between gates in the LSTM. ~\citet{singh2018hierarchical, jin2019towards} extended the contextual decomposition to produce a hierarchical clustering of input features and contextual independence of input words, respectively.~\citet{janizek2020explaining} quantified the pairwise feature interaction in DNNs by extending the explanation method of Integrated Gradients~\cite{sundararajan2017axiomatic}.
~\citet{tsang2018detecting} measured the pairwise interaction by applying the learned weights of the DNN.

Unlike above studies,~\citet{lundberg2018consistent} defined the interaction between two variables based on the Shapley value for tree ensembles.
Because Shapley value was considered as the unique standard method to estimate contributions of input words to the prediction score with solid theoretical foundations~\cite{weber1988probabilistic}, this definition of the interaction can be regarded to objectively reflect the collaborative/adversarial effects between variables \emph{w.r.t.} the prediction score.
Furthermore,~\citet{Grabisch1999AnAA, dhamdhere2019the} extended this definition to the elementary interaction component and the Shapley-Taylor index among various variables, respectively. These studies quantified the significance of interaction for each possible specific combination of variables, instead of providing a single metric to quantify interactions among all combinations of variables. In comparison, this paper aims to use a single metric to quantify the significance of interactions among all combinations of input variables, and reveal prototype features modeled by the DNN.
\subsection{Shapley values}
The Shapley value was initially proposed in the game theory~\cite{shapley1953value}. Let us consider a game with multiple players. Each player can participate in the game and receive a reward individually. Besides, some players can form a coalition and play together to pursue a higher reward. Different players in a coalition usually contribute differently to the game, thereby being assigned with different proportions of the coalition's reward.
The Shapley value is considered as a unique method that fairly allocates the reward to players with certain desirable properties~\cite{weber1988probabilistic, ancona2019explaining}. Let $N=\{1,2,\cdots,n\}$ denote the set of all players, and $2^{N}$ represents all potential subsets of $N$. A game $v:2^{N}\rightarrow \mathbb{R}$ is implemented as a function that maps from a subset to a real number. When a subset of players $S\subseteq N$ plays the game, the subset can obtain a reward $v(S)$.
Specifically, $v(\emptyset)=0$. The Shapley value of the $i$-th player $\phi_v(i|N)$ can be considered as an unbiased contribution of the $i$-th player.
\begin{small}
\begin{equation}
\label{eqn:shapleyvalue}
\phi_v(i|N)=\sum_{S\subseteq N\setminus\{i\}}\frac{(n-|S|-1)!|S|!}{n!}\Big[v(S\cup\{i\})-v(S)\Big]
\end{equation}
\end{small}~\citet{weber1988probabilistic} have proved that the Shapley value is the only reward with the following axioms.\\
\textbf{Linearity axiom:} If the reward of a game $u$ satisfies $u(S)=v(S)+w(S)$, where $v$ and $w$ are another two games.
Then the Shapley value of each player $i\in N$ in the game $u$ is the sum of Shapley values of the player $i$ in the game $v$ and $w$, \emph{i.e.} $\phi_u(i|N) =\phi_v(i|N)+\phi_w(i|N)$.\\
\textbf{Dummy axiom:} The dummy player is defined as the player that satisfies $\forall S\subseteq N\setminus\{i\}, v(S\cup\{i\})=v(S)+v(\{i\})$. In this way, the dummy player $i$ satisfies $v(\{i\})=\phi_v(i|N)$, \emph{i.e.} the dummy player has no interaction with other players in $N$.\\
\textbf{Symmetry axiom:} If $\forall S\subseteq N\setminus\{i, j\}, v(S\cup\{i\})=v(S\cup\{j\})$, then $\phi_v(i|N)=\phi_v(j|N)$.\\
\textbf{Efficiency axiom:} $\sum_{i\in N}\phi_v(i|N)=v(N)$. The efficiency axiom can ensure the overall reward can be distributed to each player in the game.

\section{Algorithm}
\subsection{Multivariate interactions in game theory}

Given a game with $n$ players $N=\{1,2,\cdots,n\}$, in this paper, we define the interaction among multiple players based on game theory. Specifically, we focus on the significance of interactions of a set of $m$ players selected from all players. Let $2^N$ denote all potential subsets of $N$, \emph{e.g.} if $N=\{a,b\}$, then $2^N=\{\emptyset, \{a\}, \{b\}, \{a,b\}\}$. The game $v:2^N\rightarrow \mathbb{R}$ is a set function that maps a set of players to a scalar value. We consider this scalar value as the reward for the set of players.
For example, we consider the DNN with a scalar output as a game, and consider the input variables as a set of players. In this way, $v(S)$, $S\subseteq N$, represents the output of the DNN when all input variables in $N\setminus S$ are replaced by the baseline value. In this paper, the baseline value is set as zero, which is similar to~\cite{ancona2019explaining}.
If the DNN is trained for multi-category classification, $v(S)$ is implemented as the score of the true category before the softmax layer.
The overall reward can be represented as $v(N)-v(\emptyset)$, \emph{i.e.} the score obtained by all input variables subtract the score obtained by the zero input.
The overall reward can be allocated to each player as $\sum_{i=1}^n\phi_v(i|N)=v(N)-v(\emptyset)$, where $\phi_v(i|N)$ denotes the allocated reward of the player $i$ over the set $N$.
Each player is supposed to obtain a high reward in the game.
In this paper, we use $\phi(i|N)$ to represent $\phi_v(i|N)$ in the following manuscript for simplicity without causing ambiguity.
Specifically, $\phi(i|N)$ can be computed as the Shapley value of $i$ based on Equation~\eqref{eqn:shapleyvalue}. The Shapley value is an unbiased method to allocate the reward, so we use the Shapley value to define the interaction.

\paragraph{Interaction between two players:}
Let us first consider the interaction between two players in the game $v$. Given two players $i$ and $j$, $\phi(i|N)$ and $\phi(j|N)$ represent their Shapley values, respectively. If players $i$ and $j$ always participate in the game together and always be absent together, then they can be regarded to form a coalition. The reward obtained by the coalition is usually different from the sum of rewards when players $i$ and $j$ participate in the game individually.
This coalition can be considered as one singleton player, denoted by $S_{ij}$. Note that the Shapley value of the coalition can also be computed using Equation~\eqref{eqn:shapleyvalue} by replacing the player $i$ with $S_{ij}$.
In this way, we can consider there are only $n-1$ players in the game, including $S_{ij}$, and excluding players $i$ and $j$, $N'=N\setminus\{i,j\}\cup S_{ij}$.
The interaction between players $i$ and $j$ is defined as the additional reward brought by the coalition $S_{ij}$ \emph{w.r.t.} the sum of rewards of each player playing individually, \emph{i.e.}
\begin{small}
\begin{equation}
\label{eqn:TwoInteraction}
B(S_{ij})=\phi(S_{ij}|N')-
\left[\phi(i|N_i)+\phi(j|N_j)\right],
\end{equation}
\end{small}
where {\small$\phi(i|N_i)$}, {\small$\phi(j|N_j)$}, and {\small$\phi(S_{ij}|N')$} are computed over the set of players {\small$N_i=N\setminus\{j\}$}, {\small$N_j=N\setminus\{i\}$}, and {\small$N'$}, respectively.
$\textbf{(1)}$ If {\small$B(S_{ij})>0$}, then players $i$ and $j$ cooperate with each other for a higher reward, \emph{i.e.} the interaction is positive.
$\textbf{(2)}$ If {\small$B(S_{ij})<0$}, then the interaction between players $i$ and $j$ leads to a lower reward, \emph{i.e.} the interaction is negative.
\begin{figure}[t]
\centering
\includegraphics[width=0.99\linewidth]{./fig/contexts.pdf}
\caption{Visualization of salient contexts $S$ given different pairs of players.}
\label{fig:contexts}
\end{figure}
{Specifically, we can measure and visualize the interaction between two players in the computer vision (CV) task. We divide the input image into $16\times 16$ grids, and consider each grid as a player. As Figure~\ref{fig:contexts} shows, for a given pair of players $i$ and $j$, we visualized contexts $S$ that boost the strength of the interaction $B(S_{ij})$, \emph{i.e.} $\Delta v(i,j,S)\cdot B(S_{ij})>0$, where $\Delta v(i,j,S)=v(S\cup\{i,j\})-v(S\cup\{i\})-v(S\cup\{j\})+v(S)$.
Let $\text{map}(S)\in \{0,1\}^m$ denote the map corresponding to the context $S$. $\text{map}_k(S)=1$ if the $k$-th grid is in the context $S$; otherwise, $\text{map}_k(S)=0$. We visualize the weighted average contexts as $\sum_S|\Delta v(i,j,S)|\cdot \text{map}(S)$.}

\paragraph{Interaction among multiple players:}
We can extend the definition of interactions to multiple players. Let us focus on a subset of players $A\subsetneqq N$.
If players in $A$ always participate in the game together, then players in $A$ form a coalition, which can be considered as a singleton player, denoted by $[A]$.
The interaction in $A$ can be measured as
\begin{small}
\begin{equation}
\label{eqn:Interaction}
{B([A])=\phi([A]|N_{[A]})-\sum\nolimits_{i\in A}\phi(i|N_i)},
\end{equation}
\end{small}
where $\phi([A]|N_{[A]})$ is the Shapley value of the coalition $[A]$. We compute the Shapley value over the player set $N_{[A]}=N\setminus A\cup\{[A]\}$. Similarly, $\phi(i|N_i)$ represents the Shapley value computed over the set $N_i=N\setminus A\cup\{i\}$, where players in $A\setminus\{i\}$ are supposed not to attend the game.
In this way, $B([A])$ reflects all interactions inside $A$, including both positive and negative interactions, which counteract each other.
For example, let the set $A=\{i,j,k,l\}$, we assume that only the interaction between players $i$ and $j$ and the interaction between players $k$ and $l$ are positive. Interactions of other coalitions are negative. Then, the positive interaction effects are counteracted by negative interaction effects in $B([A])$.
We will introduce more details in Equation~\eqref{eqn:Relationship}.

However, the significance of interactions is supposed to reflect both positive and negative interaction effects.
Therefore, we hope to propose a new metric to measure the significance of interactions, which contains both positive and negative effects among a set of players $A\subsetneqq N$.
We aim to divide all players in $A$ into a number of coalitions $\Omega=\{C_1,C_2,\cdots,C_k\}$, which are ensured to mainly encode positive interaction effects. $\Omega$ is a partition of the set $A$. In the above example, the four players $\{i,j,k,l\}$ is supposed to be divided into $\Omega=\{C_1=\{i,j\},C_2=\{k,l\}\}$.
In this case, the following equation can better reflect positive interaction effects than Equation~\eqref{eqn:Interaction}, when we measure the reward by taking each coalition $C_i$ as a singleton player.
\begin{small}
\begin{equation}
\label{eqn:MaxEffect}
B_{\text{max}}([A])=\max\nolimits_{\Omega}\sum\nolimits_{C\in \Omega}\phi(C|N_C)-\sum\nolimits_{i\in A}\phi(i|N_i)
\end{equation}
\end{small}
where $\phi(C|N_C)$ is computed over $N_{C}=N\setminus A\cup\{[C]\}$, and $\phi(i|N_i)$ is computed over $N_i=N\setminus A\cup\{i\}$.
Similarly, we can use $B_{\text{min}}([A])=\min_{\Omega}\sum_{C\in \Omega}\phi(C|N_C)-\sum_{i\in A}\phi(i|N_i)$ to roughly quantify the negative interaction effects inside the set $A$.
Thus, we define the metric $T([A])$ to measure the significance of both positive and negative interaction effects, as follows.
\begin{small}
\begin{equation}
\begin{aligned}
\label{eqn:Approximate}
T([A])&=B_{\text{max}}([A])-B_{\text{min}}([A])\\
&=\max_{\Omega}\sum_{C\in \Omega}\phi(C|N_C)-\min_{\Omega}\sum_{C\in \Omega}\phi(C|N_C)
\end{aligned}
\end{equation}
\end{small}
\textbf{Relationship between $T([A])$ and $B([A])$:}
Theoretically, the interaction $B([A])$ in Equation~\eqref{eqn:Interaction} can be decomposed into elementary interaction components $I(A)$, which was firstly proposed by~\cite{Grabisch1999AnAA}.
$I(A)$ quantifies the marginal reward of $A$, which removes all marginal rewards from all combinations of $A$.
We derive the following equation to encode the relationship between $B([A])$ and $I(A)$.
\begin{small}
\begin{equation}
\label{eqn:Relationship}
B([A])=\sum\nolimits_{A'\subseteq A, |A'|>1}I(A').
\end{equation}
\end{small}
According to Equation~\eqref{eqn:Relationship}, there are a total of $2^{m}-m-1$ different elementary interaction components inside the interaction $B([A])$, where $m$ represents the number of players inside $A$. Some of these interaction components are positive, and others are negative, which leads to a counteraction among all possible interaction components. In order to quantify the significance of all interaction components, a simple idea is to sum up the absolute values of all elementary interaction components, as follows. $B'([A])=\sum_{A'\subseteq A, |A'|>1}|I(A')|=B^+-B^-$, where $B^+=\sum_{A'\subseteq A,I(A')>0}I(A')$ and $B^-=\sum_{A'\subseteq A,I(A')<0}I(A')$, subject to $|A'|>1$.
In Equation~\eqref{eqn:MaxEffect}, $B_\text{max}([A])$ is supposed to mainly encode most positive interaction effects within $A$, which is highly related to $B^+$. Similarly, $B_\text{min}([A])$ is highly related to $B^-$ for mainly encoding negative interaction effects within $A$. Therefore, $T([A])=B_\text{max}([A])-B_\text{min}([A])$ is highly related to $B'([A])=B^+-B^-$.

\textbf{Why we need a single metric?}~\citet{Grabisch1999AnAA} designed $I(A)$, which computed $2^m$ different values of $I(A')$ for all possible $A'\subseteq A$.
In comparison, the proposed $B([A])$ is a single metric to represent the overall interaction significance among all $m$ input variables, which provides a global view to explain DNNs.

\subsection{Explanation of DNNs using multivariate interactions in NLP}
In this paper, we use the interaction based on game theory to explain DNNs for NLP tasks. Given an input sentence with $n$ words $N=\{1,2,\cdots,n\}$, each word in the sentence can be considered as a player, and the output score of DNN can be taken as the game score $v(S)$. $v(S)$ represents the output of DNN when we mask all words in the set $N\setminus S$, \emph{i.e.} setting word vectors of masked words to zero vectors. For the DNN with a scalar output, $v(S)$ denotes the output of the DNN. If the DNN was trained for the multi-category classification task, we take the output of the true category before the softmax layer as $v(S)$. Strongly interacted words usually cooperate with each other and form a prototype feature, which is memorized by the DNN for inference. Thus, the multivariate interaction can be used to analyze prototype features in DNNs.

Among all input words in $N$, we aim to quantify the significance of interactions among a subset of $m$ sequential words $A \subsetneqq N$.~\citet{chen2018lshapley} showed that non-successive words usually contain much fewer interactions than successive words. Thus, we require each coalition consists of several sequential words to simplify the implementation.
Although $T([A])$ in Equation~\eqref{eqn:Approximate} can be computed by enumerating all possible partitions in $A$, $\Omega$, the computational cost of such a method is too unaffordable.

Therefore, we develop a sampling-based method to efficiently approximate $T([A])$ in Equation~\eqref{eqn:Approximate}.
The computation of $\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)$ requires us to enumerate all potential partitions to determine the maximal value. In order to avoid such computational expensive enumeration, we propose to use $\textbf{p}=\{p_1, p_2, \cdots, p_{m-1}\}$ to represent all possible partitions. $p_i\,(0\le p_i\le 1)$ denotes the probability of the $i$-th word and the $(i+1)$-th word belonging to the same coalition.
We can sample $\textbf{g}=\{g_1, g_2, \cdots g_{m-1}\}$ based on the $\textbf{p}$, \emph{i.e.} $g_i\in\{0, 1\}, \,g_i\sim\text{Bernoulli}(p_i)$, to represent a specific partition $\Omega$. $g_i=1$ indicates the $i$-th word and the $(i+1)$-th word belong to the same coalition.
$g_i=0$ represents these words belong to different coalitions.
In this way, we can use $\text{max}_{\bf p}\mathbb{E}_{\bf g\sim \text{Bernoulli}(\bf p)}\sum_{C\in \Omega_{\bf g}}\phi(C|N_C)$ to approximate $\text{max}_\Omega\sum_{C\in\Omega}\phi(C|N_C)$, where $\Omega_{\bf g}$ denotes the partition determined by $\bf g$.

In addition, $\phi(C|N_C)$ can be approximated using a sampling-based method~\cite{castro2009polynomial}.
\begin{small}
\begin{equation}
\label{eqn:NLPphic}
\begin{aligned}
\phi(C|N_C)&=\underset{r}{\mathbb{E}}\Big[\underset{\substack{S\subseteq N_C,\\|S|=r}}{\mathbb{E}}[v(S\cup C)-v(S)]\Big]\\
&=\underset{r}{\mathbb{E}}\Big\{\underset{\substack{|S|=r+1,\\S\subseteq N_C,S\ni [C]}}{\mathbb{E}}[v(S)]-\underset{\substack{|S|=r,\\S\subseteq N_C,S\not\ni [C]}}{\mathbb{E}}[v(S)]\Big\}
\end{aligned}
\end{equation}
\end{small}
where $r$ represents the number of players fed into the DNN, $N_C=N\setminus A\cup\{[C]\}$.
In this way, we can approximate $\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)$ as $\max_{\bf p}\mathbb{E}_{{\bf g}\sim \text{Bernoulli}({\bf p})}\sum_{C\in\Omega_{\bf g}}\phi(C|N_C)$, as follows.
\begin{small}
\begin{equation}
\begin{aligned}
\label{eqn:NLPphic1}
\mathcal{L}\!\!&=\!\!\!\!\underset{{\bf g}\sim \text{Bernoulli}({\bf p})}{\mathbb{E}}\!\!\sum_{C\in\Omega_{\bf g}}\!\!\!\!\phi(C|N_C),
\max_{\bf p}\mathcal{L}\!=\!
\max_{\bf p}\sum_{i\in A}\\
&\underset{r}{\mathbb{E}}
\Bigg\{\underset{{\bf g}\sim \text{Bernoulli}({\bf p})}{\mathbb{E}}\!\!\!\!\!\!\!\!\lambda_i({\bf g})\Bigg[
\underset{\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\ni i, |S|=r+1}}{\mathbb{E}}\!\!\![v(S)]-\!\!\!\!\!\!\underset{\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\not\ni i, |S|=r}}{\mathbb{E}}\!\!\![v(S)]
\Bigg]\Bigg\}
\end{aligned}
\end{equation}
\begin{equation}
\label{eqn:gradient}
\begin{aligned}
\frac{\partial \mathcal{L}}{\partial p_i}\!\!=\!\!&\sum_{j\in A}\!\mathbb{E}_r\Big\{\underset{\substack{{\bf g}\sim\text{Bernoulli}({\bf p})\\g_i=1}}{\mathbb{E}}\!\!\!\!\!\!\!\lambda_i({\bf g})\!\Big[\!\!\underset{
\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\ni j, |S|=r+1}}{\mathbb{E}}\!\!\!\!\!\!\![v(S)]\!\!-\!\!\!\!\!\underset{\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\not\ni j, |S|=r}}{\mathbb{E}}\!\!\!\!\![v(S)]
\Big]\\
-&\!\!\!\!\!\!
\underset{\substack{{\bf g}\sim\text{Bernoulli}({\bf p})\\g_i=0}}{\mathbb{E}}\!\!\!\!\!\!\!\!\lambda_i({\bf g})\Big[\!\!\!\!\underset{
\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\ni j, |S|=r+1}}{\mathbb{E}}\!\!\!\!\![v(S)]-\!\!\!\!\!
\underset{\substack{S\in \text{Sam} (\Omega_\textbf{g}),\\S\not\ni j, |S|=r}}{\mathbb{E}}\!\!\!\!\![v(S)]
\Big]\!\!
\Big\}
\end{aligned}
\end{equation}
\end{small}where $\lambda_i({\bf g})=\frac{1}{|C_i|}$, and $S\in \text{Sam}(\Omega_{\bf g})$ represents the sampled set of words, which contains all words in $N\setminus A$ and the coalition $C_i$. $C_i\in\Omega_{\bf g}$ denotes the coalition determined by ${\bf g}$ that contains the word $i$. We learn $\bf p$ to maximize the above equation.
For example, let $A=\{x_1,x_2,x_3,x_4,x_5,x_6\}$, and the sampled ${\bf g}=\{g_1=1, g_2=1, g_3=0, g_4=0, g_5=1\}$. We consider the first three words in $A$ as a coalition, and the last two words in $A$ as another coalition, \emph{i.e.} $\Omega_{\bf g}=\{C_1=\{x_1,x_2,x_3\}, C_2=\{x_4\}, C_3=\{x_5,x_6\}\}$.
The set $S$ is supposed to be sampled over these coalitions in $\Omega_{\bf g}$, instead of over individual words. In this case, we have $\lambda_1({\bf g})=\lambda_2({\bf g})=\lambda_3({\bf g})=1/3$, $\lambda_4({\bf g})=1$, $\lambda_5({\bf g})=\lambda_6({\bf g})=1/2$.
In order to learn $p_i$, we compute $\partial \mathcal{L}/\partial p_i$, which is given in Equation~\eqref{eqn:gradient}.
In this way, we can use $\max_{\bf p}\mathcal{L}$ in Equation~\eqref{eqn:NLPphic1} to approximate the solution to $\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)$. Similarly, we can use $\min_{\bf p}\mathcal{L}$ to approximate the solution to $\min_\Omega\sum_{C\in\Omega}\phi(C|N_C)$, thereby obtaining $T([A])=\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)-\min_\Omega\sum_{C\in\Omega}\phi(C|N_C)$.

\begin{table}[t]
\begin{center}
\caption{Accuracy of the estimated partition.}
\label{tab:toy_example}
\resizebox{0.85\linewidth}{!}{
\begin{tabular}{lccc}
\toprule
& \tabincell{c}{Add-Multiple\\ Dataset} & \tabincell{c}{AND-OR\\ Dataset} & \tabincell{c}{Exponential\\ Dataset} \\
\midrule
Baseline 1 & 0.500 & 0.503 & 0.506 \\
Baseline 2 & 1.000 & 0.996 & 1.000\\
Baseline 3 & 1.000 & 0.523 & 0.846\\
\midrule
Our method & \textbf{1.000} & \textbf{0.999} & \textbf{1.000}\\
\bottomrule
\end{tabular}
}
\end{center}
\end{table}

\textbf{Comparison of computational cost:} Given a set of $m$ words in a sentence with $n$ words, the computation of $T([A])$ based on its definition in Equation~\eqref{eqn:Approximate} and Equation~\eqref{eqn:shapleyvalue} is NP-hard, \emph{i.e.} $O(2^nm)$\footnote[1]{Please see supplementary materials for the proof.}.
Instead, we propose a polynomial method to approximate $T([A])$ with computational cost $O(K_1K_2K_3nm)$\footnotemark[1], where $K_1$ denotes the number of updating the probability $\bf p$, $K_2$ and $K_3$ represent the sampling number of $\bf g\sim\text{Bernoulli}({\bf p})$ and $S\in\text{Sam}(\Omega_{\bf g})$, respectively.
In addition, we conducted experiments to show the accuracy of the estimation of $T([A])$ increases along with the increase of the sampling number.

\section{Experiments}
\textbf{Evaluation of the correctness of the estimated partition of coalitions:}
The core challenge of evaluating the correctness of the partition of coalitions was that we had no ground-truth annotations of inter-word interactions, which were encoded in DNNs. To this end, we constructed three new datasets with ground-truth partitions between input variables, \emph{i.e.} \emph{the Addition-Multiplication Dataset}, \emph{the AND-OR Dataset}, and \emph{the Exponential Dataset}.

\begin{figure}
\begin{center}
\includegraphics[width=0.75\linewidth]{./fig/toy_example6.pdf}
\captionof{figure}{Convergence of $p_i$ during the training process.}
\label{fig:toy_example}
\end{center}
\end{figure}

$\bullet$ \emph{The Addition-Multiplication Dataset} contained $10000$ addition-multiplication models, each of which only consisted of addition operations and multiplication operations. For example, $y=f({\bf x})=x_1+x_2\times x_3+x_4\times x_5+x_6+x_7$.
Each variable $x_i$ was a binary variable, \emph{i.e.} $x_i\in\{0,1\}$. For each model, we selected a number of sequential input variables to construct $A$, \emph{e.g.} $A=\{x_2,x_3,x_4,x_5\}$. We applied our method to extract the partition of coalitions \emph{w.r.t.} $y$, which maximized $\sum_{C\in\Omega}\phi(C|N_C)$.

\begin{figure*}[!ht]
\centering
\includegraphics[width=0.99\linewidth, height=0.16\linewidth]{./fig/InstabilityAndDifferenceAndInter4.pdf}
\caption{(a) The decrease of error of $T([A])$ during the learning process. We obtained accurate $T([A])$ after 100 epochs. (b) The instability of interactions with different numbers of sampling. When the number of sampling was larger than $1000$, the instability of $T([A])$ was less than $0.1$ on all models and datasets.
(c) The significance of interactions in the computation of intermediate-layer features gradually increased during the forward propagation in the BERT.}
\label{fig:instabilityAndCorrectness}
\end{figure*}

$\bullet$ \emph{The AND-OR Dataset} contained $10000$ AND-OR models, each of which only contained AND operations and OR operations. For simplicity, we used $\&$ to denote the AND operation, and used $|$ to denote the OR operation. For example, $y=f({\bf x})=x_1\; |\;( x_2\; \&\; x_3)\;|\;( x_4\; \&\; x_5)\;|\; x_6\; |\; x_7$.
Each variable $x_i$ was a binary variables, \emph{i.e.} $x_i\in\{0,1\}$.
For each model, we selected a number of sequential input variables to construct $A$, \emph{e.g.} $A=\{x_2,x_3,x_4,x_5\}$. We used our method to extract the partition of coalitions, which maximized $\sum_{C\in\Omega}\phi(C|N_C)$.

$\bullet$ \emph{The Exponential Dataset} contained $10000$ models, each of which contained exponential operations and addition operations.  For example, $y=f({\bf x})=x_1^{x_2}+x_3^{x_4}+x_5+x_6$.
Each variable $x_i$ was a binary variables, \emph{i.e.} $x_i\in\{0,1\}$.
For each model, we selected a number of sequential input variables to construct $A$, which always contained the exponential operations, \emph{e.g.} $A=\{x_1,x_2,x_3,x_4\}$. We used our method to extract the partition of coalitions, which maximized $\sum_{C\in\Omega}\phi(C|N_C)$.

We compared the extracted partition of coalitions with the ground-truth partition of coalitions. Variables in multiplication operations and AND operations usually had positive interactions, which were supposed to be allocated into one coalition during $\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)$. Variables in OR operations had negative interactions, which were supposed to be allocated to different coalitions during $\max_\Omega\sum_{C\in\Omega}\phi(C|N_C)$. Note that the addition operation did not cause any interactions. Therefore, we did not consider interactions of addition operations when we evaluated the correctness of the partition.
In the above example of addition-multiplication model, if we considered the set $A=\{x_2,x_3,x_4,x_5\}$, then the ground-truth partition was supposed to be either $\{\{x_2,x_3\},\{x_4,x_5\}\}$ or $\{\{x_2,x_3,x_4,x_5\}\}$.
For the above example of AND-OR model and $A=\{x_2,x_3,x_4,x_5\}$, the ground-truth partition should be $\{\{x_2,x_3\},\{x_4,x_5\}\}$.
For each operation, if the proposed method allocated variables of an operation in the same way as the ground-truth partition, then we considered the operation was correctly allocated by the proposed method. The average rate of correctly allocated operations over all operations was reported in Table~\ref{tab:toy_example}.

\begin{figure*}[!ht]
\centering
\includegraphics[width=0.87\linewidth]{./fig/example8.pdf}
\caption{Prototype features modeled by the BERT trained using the SST-2 dataset. Grey box indicated words in $A$. We used different colors to indicate the extracted coalitions. Words in the same coalition formed a prototype feature modeled by the DNN. The set of words $A$ was randomly selected.}
\label{fig:interaction_example}
\end{figure*}

We also designed three baselines.
\textbf{Baseline 1} randomly combined input variables to generate coalitions.~\citet{lundberg2018consistent} defined the interaction between two input variables, which was used as the \textbf{Baseline 2}.~\citet{li2016understanding} proposed a method to estimate the importance of input variables, which was taken as \textbf{Baseline 3}. \textbf{Baseline 2} and \textbf{Baseline 3} merged variables whose interactions were greater than zero to generate coalitions.

Table~\ref{tab:toy_example} compares the accuracy of coalitions generated by our method and that of other baselines. Our method achieved a better accuracy than other baselines.
Figure~\ref{fig:toy_example} shows the change of the probability $p_i$ during the training process.
Our method successfully merged variables into a single coalition in multiplication operations and AND operations. Besides, our method did not merge variables in OR operations. Note that the probability $p_3$ of the addition operation in Figure~\ref{fig:toy_example} did not converge. It was because the addition operation did not cause any interactions, \emph{i.e.} arbitrary value of $p_3$ satisfied the ground-truth.

\textbf{Evaluation of the accuracy of $T([A])$:}
We applied our method to DNNs in NLP to quantify the interaction among a set of input words. We trained DNNs for two tasks, \emph{i.e.} binary sentiment classification based on the SST-2 dataset~\cite{socher2013recursive} and prediction of linguistic acceptability based on the CoLA dataset~\cite{warstadt2018neural}.
For each task, we trained five DNNs, including the BERT~\cite{devlin2018bert}, the ELMo~\cite{peters2018deep}, the CNN proposed by~\cite{kim2014convolutional}, the two-layer unidirectional LSTM~\cite{hochreiter1997long}, and the Transformer~\cite{vaswani2017attention}.
For each sentence, the set of successive words $A$ was randomly selected.

We compared the extracted significance of interactions $T([A])$ with the accurate significance of interactions. The accurate significance of interactions was quantified based on the definition in Equation~\eqref{eqn:Approximate} and Equation~\eqref{eqn:shapleyvalue}, which was computed by enumerating all partitions $\Omega$ and all subsets $S$ with an extremely high computational cost. Considering the unaffordable computational cost, such evaluation could only be applied to sentences with less than $12$ words.
Figure~\ref{fig:instabilityAndCorrectness} (a) reports the error of $T([A])$, \emph{i.e.} $|T_\text{truth}([A])-T([A])|$, where $T_\text{truth}([A])$ was the true interaction significance accurately computed via massive enumerations. We found that the estimated significance of interactions was accurate enough after the training of $100$ epochs.

\textbf{Stability of $T([A])$:}
We also measured the stability of $T([A])$, when we computed $T([A])$ multiple times with different sampled sets of $\bf g$ and $S$. The instability of $T([A])$ was computed as
$\text{instability}=\mathbb{E}_{I}[\frac{\mathbb{E}_{u,v:u\neq v}|T_{(u)}([A])-T_{(v)}([A])|}{\mathbb{E}_{w}|T_{(w)}([A])|}]$, where $T_{(u)}([A])$ denotes the $u$-th computation result of $T([A])$.
Figure~\ref{fig:instabilityAndCorrectness} (b) shows the instability of interactions in different DNNs and datasets.
Experimental results showed that interactions in CNN and ELMo converged more quickly. Moreover, instability decreased quickly along with the increase of the sampling numbers. Our metric was stable enough when the number of sampling was larger than $2000$.

\begin{figure*}
  \centering
  \includegraphics[width=0.99\linewidth]{./fig/fail_example5.pdf}
  \caption{Positive and negative interaction effects of sentences, which were mistakenly classified by the DNN. These results reflect prototype features towards correct and incorrect predictions, and representation flaws of the DNN.}
  \label{fig:fail_example}
\end{figure*}

\textbf{Prototype features modeled by DNNs:}
We analyzed results obtained by the proposed method to explore prototype features modeled by DNNs. Figure~\ref{fig:interaction_example} shows several results on the BERT. We found that: (1) For an entire constituent in the sentence, such as a short clause or a noun phrase, the BERT usually took the whole constituent as a single coalition, \emph{i.e.} a prototype feature. (2) For the set of words that contained conjunctions or punctuations, such as ``and'', ``or'', ``,'', \emph{etc.}, the BERT divided the set at conjunctions or punctuations to generate different coalitions (prototype features). (3) For the constituent modified by multiple adjectives and adverbs, the BERT usually merged all adjectives and adverbs into one coalition, and took the modified constituent as another coalition.
Such phenomena fit the syntactic structure parsing according to human knowledge.

\textbf{Interactions \emph{w.r.t.} the intermediate-layer feature:}
We could also compute the significance of interactions among a set of input words \emph{w.r.t.} the computation of an intermediate-layer feature $\bm f$.
The reward was computed using the intermediate-layer feature. Let ${\bm f}_N$ and ${\bm f}_S$ represent intermediate-layer features obtained by the set of input variables $N$ and $S$, respectively.
Since the intermediate-layer feature usually can be represented as high dimensional vertors, instead of a scalar value. We computed the output $v(S)$ as $v_{\bm f}(S)=\langle{\bm f}_N, {\bm f}_S\rangle/\Vert{\bm f}_N\Vert_2$, where $\Vert{\bm f}_N\Vert_2$ was the L2-norm and was used for normalization.
Figure~\ref{fig:instabilityAndCorrectness} (c) shows the significance of interactions computed with output features of different layers. The significance of interactions among a set of input variables increased along with the layer number. This phenomenon showed that the prototype feature was gradually modeled in deep layers of the DNN.

\textbf{Mining prototype features towards incorrect predictions of DNNs:} The proposed method could be used to analyze sentences, which were mistakenly classified by the DNN. We used the BERT learned using the SST-2 dataset.
Figure~\ref{fig:fail_example} shows two partitions that maximized and minimized the Shapley value of coalitions for each sentence.
For the partition that maximized the Shapley value, the generated coalitions were usually used to explain prototype features toward incorrect predictions of the DNN. In comparison, coalitions generated by minimizing the Shapley value represented prototype features towards correct predictions. For example, we tested two different sets of input variables $A$ for the sentence ``on the heels of ...... to creepiness.''. This sentence was incorrectly predicted to be negative by the DNN. As Figure~\ref{fig:fail_example} shows, the DNN used ``humorless horror movie'' as a prototype features, which led to the incorrect negative prediction. In comparison, if we minimized the Shapley value, we got the coalition ``is to be recommanded for its straight-ahead approach to creep'', which towards the correct positive prediction.

\section{Conclusion}
In this paper, we have defined the multivariate interaction in the DNN based on game theory. We quantify the interaction among multiple input variables, which reflects both positive and negative interaction effects inside these input variables. Furthermore, we have proposed a method to approximate the interaction efficiently. Our method can be applied to various DNNs for different tasks in NLP.
Note that the quantified interaction is just an approximation of the accurate interaction in Equation~\eqref{eqn:Approximate}. Nevertheless, experimental results have verified high accuracy of the approximation.
The proposed method can extract prototype features modeled by the DNN, which provides a new perspective to analyze the DNN.

\section*{Acknowledgements}
This work is partially supported by National Natural Science Foundation of China (61906120 and U19B2043).

\section*{Ethical Impact}
This study has broad impacts on the understanding of signal processing in DNNs. Our work provides researchers in the field of explainable AI with new mathematical tools to analyze DNNs. Currently, existing methods mainly focus on interactions between two input variables. Our research proposes a new metric to quantify interactions among multiple input variables, which sheds new light on the understanding of prototype features in a DNN. We also develop a method efficiently approximate such interactions. As a generic tool to analyze DNNs, we have applied our method to classic DNNs and have obtained several new insights on signal processing encoded in DNNs for NLP tasks.