\section{Introduction}

Magnetic resonance imaging (MRI) is the gold standard technique to diagnose brain tumors, such as glioblastoma (GBM). While being able to reliably highlight areas of sufficiently high tumor cell concentration in GBMs, it can lack the capacity to visualize areas of lower tumor cell density at the tumor border and most importantly areas of diffuse tumor infiltration, a key biological property of GBM. Current radiotherapy plans try to account for the unknown infiltration by targeting a uniformly extended volume around the tumor outlines visible in MRI. While decreasing the probability of tumor recurrence, such a treatment planning has an obvious drawback of unnecessarily damaging healthy tissue, which in turn has a negative impact on the patient’s life quality. Personalizing the target of radiotherapy by complementing the MRI scans with individual tumor simulation, that models a complete spatial distribution of tumor cell concentration, could preserve healthy tissue and reduce the risk for secondary malignancies and side effects \cite{Le_tmi,le_miccai,jana_tmi,hormuth2021image}.

Conventional approaches for simulation-based personalization attempt to model the tumor growth for each individual patient using differential equation solvers. The personlaization is achieved by solving the inverse problem - identifying tumor model's parameters best matching tumor signal from MRI. However, utilizing the numerical solvers for solving an inverse problem still results in extreme runtimes which obstruct transfer into clinical practice. To address this issue, highly efficient model- and data-driven approaches were developed over the recent years \cite{subramanian2020multiatlas,scheufele2019coupling,hormuth2018mechanically,hormuth2021image,jana_tmi,petersen2019deep,petersen2021continuous}. The time for solving an inverse problem using numerical model-based solvers hit one hour of compute \cite{subramanian2020multiatlas,scheufele2019coupling,hormuth2021image}. For data-driven approaches the computing time can be reduced to minutes \cite{petersen2019deep,petersen2021continuous,Ezhov_2019,Ezhov2020GeometryawareNS}, but they all rely on some type of a neural network to predict tumor progression. Neural networks are known to be non-robust when extrapolated to out-of-training distributions - the key obstacle in integrating such algorithms into safety-critical applications.

In this paper, we propose an image retrieval based approach that performs a query of
patient specific scans to a database of synthetic tumors, returning the closest resemblance of the patient’s tumor. As a baseline, this image retrieval process is implemented via a primitive iterative pair-wise comparison. Further, we investigate how a retrieval in low dimensional embeddings of the simulations -- that we obtain by using downsampling as well as autoencoders, variational autoencoders, unsupervised hashing, and radiomics feature representation -- can improve runtimes. As a result, our work shows that our query approach can yield accurate and, depending on the chosen optimization, also deterministic results in the order of seconds (correspondingly, minutes for the inverse problem).

\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figs/db_query.pdf}
\caption{\small{Schematic illustration of the proposed pipeline. First, we register a patient MRI scan to the brain atlas. Next, we use the obtained forward transformation matrix to morph the patient segmentation to the atlas space. Then, we perform a query to a database of 50k pre-simulated tumors in the atlas space. Essentially, in this step we find the closest to the query sample binary volume in the database. Finally, we transform a tumor simulation, corresponding to the best matching binary volume, back to the patient space.}}
\label{qualcomp}
\end{figure}
\section{Conclusion}

We present a method for inferring a tumor simulation from information available on medical scans relying on a simple database query strategy. A reader might say that the proposed method is a joke due to its offensively simple nature. While we would again agree with such a statement, we do believe that the method provides a generalizable, fast, and robust solution to the inverse problem. The runtimes for solving the inverse problem are in the order of a few minutes which is faster than any existing inverse modeling method. Even though the used database is composed of simulations from the reaction-diffusion model that might not be sufficient to provide a close match to real tumor growth in general, it is evident that the proposed method is generalizable to more sophisticated models (one just needs to resimulate tumor database in the atlas space with more complicated tumor model). As opposed to data-driven approaches, the method does not fail at extrapolation as it is based on a plain compute of the overlap measure between two binary volumes.\\

\bibliographystyle{splncs}

\section{Method}
\subsection{Tumor growth model.}

Our forward model is based on the most used within the medical imaging community mathematical formalism - the reaction-diffusion partial differential equation. The equation describes the evolution of tumor cell density through proliferation and diffusion of the cells:

\begin{equation}
\frac{\partial u}{\partial t}= \nabla (\textbf{D}\nabla u)+\rho u(1-u), \quad in \quad\Omega \\
\end{equation}
\begin{equation}
\textbf{n} \cdot \nabla u = 0, \quad in \quad\Gamma_{\Omega}.
\end{equation}

Here, the first term on the right-hand site of the equation denotes the diffusion process, where \textbf{D} is the diffusion tensor. The second term models the logistic proliferation of cells, and $\rho$ denotes the rate of cell proliferation. We impose the Neumann boundary condition on the boundary of the simulation domain $\Gamma_{\Omega}$, while $\Omega$ denotes the simulation domain itself, i.e. volume of the brain, and \textbf{n} is the unit vector orthogonal to $\Gamma_{\Omega}$.
As an initial condition for the cell density $u$ we use a seed point at a spatial location $\pmb{r^*}$: $u(\pmb{r},0)=u_0$ if $\pmb{r}=\pmb{r^*}$, $u(\pmb{r},0)=0$ if $\pmb{r}\neq\pmb{r^*}$. The solver outputs a simulated tumor density $u$ of size $128^3$.

We assume that visible outlines in an MRI image correspond to isolines in $u$, and to link the tumor simulations with an imaging signal, we threshold the simulations at $u_{T1Gd}, u_{FLAIR}$. This allows us to reproduce the binary segmentation obtained from an MRI scan (T1Gd, FLAIR) delineation \cite{le_miccai,Ezhov_2019}.

\subsection{Inverse model.}

The proposed inverse model pipeline consists of the following steps, Fig. 1:

1. A patient MRI scan is registered to the
atlas brain anatomy \cite{rohlfing2010sri24}, yielding a transformation matrix $M$.

2. $M$ is used to morph the corresponding segmentation $Y$ into the atlas space resulting in $Y_a$.

3. A query is made to a database of synthetic tumors which were simulated beforehand in the brain atlas. This query finds the closest match\footnote[1]{The closest match can be determined in different ways which we are comparing in this study} to $Y_a$, and returns the corresponding tumor simulation in the atlas space.

4. As a last step, utilizing the inverse of $M$, the atlas simulation is morphed back to the patient's brain anatomy to obtain the patient-specific simulation.

\subsection{Query to a database of synthetic tumors. }
The input MRI scans consist of two segmentations while the synthetic tumors are a spatial distribution containing continuous values in the range from 0 to 1. Therefore, we have to extract the desired segmentations from the raw synthetic tumor data. For this purpose, we threshold the
synthetic tumors at $u_{FLAIR}=0.2$ to receive the FLAIR and at $u_{T1Gd}=0.6$ for the T1Gd segmentations \cite{Le_tmi,martens2021deep}.

After preprocessing the dataset to make it comparable, we need to select a fitting metric. The Dice coefficient is a common measure for medical imaging segmentations, we therefore utilized it as the main metric.

Even though we believe that a for-loop is all you need, in order to please reviewers spoiled with deep learning sugar, we evaluated different learnable strategies to perform a query to a database of simulated tumors. For some of these query strategies, in order to evaluate similarities between the compressed vectors, we use the L2 = $\left\| \cdot \right\|_2 $ metric, as Dice is barely meaningful e.g. for comparing latent representations. Furthermore, since we have to compare two segmentations for each tumor, it is possible that the separate queries result in different best matches. Hence, the two query results have to be combined to determine a joint best match. For simplicity, we decided to simply add the two individual similarity values and choose the highest combined value for Dice ($Dice_{T1Gd}+Dice_{FLAIR}$) or the lowest for L2 ($L2_{T1Gd}+L2_{FLAIR}$) as the best match.

\subsubsection{Baseline for-loop search.} The most basic way to implement the query is by naively looping over the dataset and performing a pair-wise comparison between input and database tumors.
We implemented this basic brute-force loop to create a ground truth that serves
as a reference for later improvements and experiments. Since the individual comparisons
are not dependent on each other, we implemented the query by parallelizing the loop to
shorten its runtime.

Although parallelization can help to reduce the time complexity of the database query, the
speed of this operation is still limited by the significant memory usage. This restriction also
prohibits the use of efficient query frameworks. For instance, the FAISS library is designed
for vector sizes in the range of roughly 20 - 2000 \cite{johnson2019billion}. However, a plain flattening of our 3D volume of $128^3$ voxels exceeds this suggested limit by a factor of more than 1000. Therefore, we decided to experiment with a compression of the data.

As a first step, we try a naive downsampling of the volume to lower resolutions. In our case, we tried a down sampling to $64^3$ and $32^3$ voxels. We downsample the original volume by applying a spline interpolation of order 0 to the original volume \cite{virtanen2020scipy}. We then compare downsampled database with a downsampled query volume, and choose the one that has the highest Dice in the downsampled volume.

\subsubsection{Autoencoders.} To obtain even more compact data representation, we experimented with convolutional autoencoders (AE) and variational autoencoders (VAE). We examined if these techniques are capable of creating a problem-specific encoding that might be able to sufficiently preserve similarity relationships between tumors, while drastically reducing the runtime. The core idea for our use-case is to train the AE to extract meaningful properties of the tumors and then apply the encoder function to the dataset $\{Y\}$, yielding an encoded dataset $\{Y_{enc}\}$. For each query, we utilize $\{Y_{enc}\}$ and encode the input tumor using the same encoder, resulting in a data size drop from $128^3$ to the dimension of the latent space (1024). Since the segmentations of T1Gd and FLAIR scans differ in their properties, two separate AEs are needed, leading to two encoded datasets $\{Y_{enc}^{T1Gd}\}$ and $\{Y_{enc}^{FLAIR}\}$.

To further improve the performance of our AEs, we experimented with VAEs. Basic AEs have the potential disadvantage of limited control over the latent space. VAEs enforce a distribution over the latent space by penalizing deviations from a previously selected distribution during training, which can produce more continuous and organized latent spaces while still providing adequate reconstructions \cite{kingma2013auto}. Usually, this property is used to generate previously unseen output by randomly sampling from the latent space. However, we investigated if this characteristic is also beneficial for similarity preservation.

In the end, we compare latent representation of the database with a latent representation of the query volume, and choose the one that has the lowest L2.

\subsubsection{Unsupervised hashing.} An alternative compression strategy is to use a hashing - compression to a bit string representation \cite{hashing1,hashing2,hashing3,hashing4}. We resorted to the most recent deep unsupervised hashing technique \cite{liAAAI2021} for our experiments. The methodological idea is a replica of the AE with a difference in how the the latent space is formed. In \cite{liAAAI2021}, the continuous output of the AE's encoder is quantized via a so-called "Bi-half" hash layer into a binary representation. This representation serves then as an input to the AE's decoder. All other network elements are identical to the conventional AE setup. Similarly to the AE and VAE, we choose the best match that has the lowest L2 when compared to the binary representation of the query sample.

\subsubsection{Radiomics shape feature representation.} The ultimate compression strategy is to use single value characteristics of a 3D shape and morphology. We experimented with shape features from the pyradiomics library \cite{van2017computational}, namely voxel volume, elongation, major axis length, etc. Additionally, we computed the tumor's center of mass to account for its location. For the query, we used the following strategy. We selected 1000 samples from our database which have the closest center of mass to our query sample in the atlas reference space. Out of these 1000 samples, we selected the one that is closest to the query sample under L2 measure evaluated for the shape features: $\sum_i \left\| f_q^i, f_d^i \right\|_2$. Here, the summation is over the radiomics features, $f_q^i$ is the feature representation computed for the query sample, and $f_d^i$ is the feature computed for the database sample.

\subsection{Implementation.}

The base architecture for our AE is a convolutional neural network. In the encoder, we scale down
the volume layer-by-layer using strided 3D convolutions (stride=2, kernel size=3) followed by
a convolution without stride, while increasing the number of channels. After downscaling
the volume, we flatten the features and apply a linear layer to obtain a latent representation of size 1024. Afterwards, the decoder mirrors this process with strided 3D transposed convolutions. We executed multiple test runs searching for the best architecture that suggested that a downscaling to a resolution of $16^3$ before flattening the volume and 1024 as the size of the latent representation yield the most reasonable results. We use this resolution throughout all consecutive experiments. For comparability, our VAE network replicates the AE network (except the linear layer before the latent space that is replaced by 2 separate linear layers which output the mean and the logarithm of the variance). AE, VAE, and the hashing method were trained on a dataset of 1500 tumors (larger training set sizes showed no significant improvements) along with a validation set of size 150.

We used the Advanced Normalization Tools (ANTs) \cite{avants2009advanced} for the registration. with basic settings. For the looping-based queries, we performed the compute on AMD EPYC 7452 CPU and parallelized the search over 32 processes. For AE, VAE, and BF-HASH, the training and evaluation were performed on an   NVIDIA QUADRO P8000 GPU.

\section{Experiments}
\subsubsection{Data.}
The basis for our work is a dataset of 50,000 synthetic tumor volumes of $128^3$ resolution. In order to generate synthetic tumors, a common brain anatomy is required that is able to represent a broad range of realistic human brains. For this purpose, the atlas brain anatomy introduced in \cite{rohlfing2010sri24} was utilized. This represents a statistical average for the relevant distributions of cerebrospinal fluid (CSF), white matter (WM) and gray matter (GM), combined into a single atlas space. In this atlas space, a tumor can then be simulated using the reaction diffusion model by setting values for the initial location $ x, y, z $ , the proliferation rate $\rho$, the diffusion rate in WM $D_w$ as well as the simulation end time $T_{end}$, which corresponds to the tumor’s age. The synthetic dataset was generated by randomly sampling patient-specific parameters from the subsequent ranges and feeding the resulting parameter sets
$\theta = \{ x, y, z, D_w , \rho, T_{end} \}$ to the reaction diffusion model. Analogous to \cite{ezhov2021learn}, we discarded unrealistic in size tumors based on the range of real tumor sizes (BraTS dataset \cite{menzebrats}).

In all experiments, we evaluate the inverse model pipeline on 62 real MRI scans of patients diagnosed with gliomas for which T1Gd and FLAIR MRI scanning were performed. For these patients, the corresponding tumor annotations were obtained by expert radiologists' labeling. In addition, for evaluation on synthetic data, we used 1k tumors simulated in the atlas space but not present in the query database.

\begin{table}
  \centering
  \begin{tabular}{ |c||c|c|c|c|c|c|c| }
         \hline
          & direct query & DS 64 & DS 32 & VAE & AE & BF-HASH &  RF\\
         \hline
         Top-1              & 100\% & 88.7\% & 66.1\% & 41.9\% & 33.8\%  & 32.2\% & 4.92\% \\
         Top-5              & 100\% & 100\% & 96.7\% & 61.2\% & 64.5\%  & 67.7\% & 4.92\%\\
         Top-15             & 100\% & 100\% & 100\% & 80.6\% & 77.4\%  & 82.2\% & 9.83\% \\
         \hline
         Runtime     & 340s & 130s & 116s & 34s & 22s & 22s & 2.8s\\
         \hline
  \end{tabular}
  \centering
  \caption{Comparison between different best match query strategies: looping over the dataset with resolution $128^3$ (direct query), looping over downsampled to $64^3$ dataset (DS 64), looping over downsampled to $32^3$ dataset (DS 32), using embeddings obtained by the VAE and AE, using the radiomics features (RF), and using hashing (BF-HASH). Top N accuracy represents a measure whether the predicted best match falls in the top N closest predictions. For the query using direct, DS64, and DS32 methods the Dice score is used, while for the remaining methods we use the $L2$ measure. We expectedly observe that an improvement in runtime comes at a cost of query accuracy.}
\end{table}

\begin{table}
 \centering
 \begin{tabular}{c|c|c|c||c|c|c}

        \hline
        & \multicolumn{3}{c}{Real} &  \multicolumn{3}{c}{Synthetic} \\
        \cline{2-7}
         & Top-1 & Top-5 & Top-15  & Top-1 & Top-5 & Top-15 \\
        \hline
        direct query& 100\% & 100\% & 100\%  & 100\% & 100\% & 100\% \\
        \hline
        DS64 & 88.7\% & 100\% & 100\%   & 93.5\% & 100\% & 100\% \\
        \hline
        DS32 & 66.1\% & 96.7\% & 100\%  & 78.9\% & 99.6\% & 100\% \\

        \hline
 \end{tabular}
 \centering
 \caption{Comparison between different loop-based query strategies on real (62 cases) and synthetic (1000 cases): looping over the dataset with resolution $128^3$ (direct query), looping over downsampled to $64^3$ dataset (DS 64), looping over downsampled to $32^3$ dataset (DS 32). Top-N accuracy represents a measure whether the predicted best match falls in the top N closest predictions. The average DICE (sum of DICEs for T1Gd and FLAIR) for best matching samples is 1.07 for real, and 1.69 for synthetic samples.}
\end{table}

\subsubsection{Quantitative comparison between query methods.} Table 1 compares the presented database query approaches on the real dataset: 1) simple looping over the dataset with original resolution $128^3$ (direct query), 2) looping over downsampled to $64^3$ dataset (DS 64), 3) looping over downsampled to $32^3$ dataset (DS 32), 4) using embeddings obtained by the VAE and AE, 5) using the radiomics features (RF), and 6) and unsupervised hashing (BF-HASH). The top-N accuracy is computed with respect to the best match samples found by the baseline. As expected, we observe that an improvement in runtime comes at a cost of losing accuracy of similarity preservation.

\begin{figure}[ht]
\centering
\includegraphics[width=1.0\textwidth]{figs/real.pdf}
\caption{\small{DICE between real tumor segmentation (62 patients) and best matching synthetic tumor in the atlas space. DICE is computed for T1Gd (upper row), FLAIR (middle row), and combined T1Gd+FLAIR (down row).}}
\label{distr_real}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=1.0\textwidth]{figs/syn.pdf}
\caption{\small{DICE between synthetic tumor segmentation (1000 examples) and best matching synthetic tumor in the atlas space. DICE is computed for T1Gd (upper row), FLAIR (middle row), and combined T1Gd+FLAIR (down row).}}
\label{distr_syn}
\end{figure}

Tab.2 showcases how the difference between real and synthetic data distributions affects the query accuracy. Given there is an evident gap between real tumor progression and the tumor growth trajectory modeled by the reaction-diffusion model, we expectedly observe ca. 10 percent drop of the accuracy. Fig. \ref{distr_real} and \ref{distr_syn} demonstrate the same difference in real and synthetic data distributions for every example used in the study.
In a way, in future works probing more complicated tumor descriptions, such analysis can serve as a measure of plausibility of a tumor growth model.

\begin{figure}[hb]
\centering
\includegraphics[width=1.0\textwidth]{figs/tumors.pdf}
\caption{\small{Qualitative examples of the tumor simulations obtained from our proposed pipeline. The brain images correspond to five different patients.}}
\label{qualcomp}
\end{figure}

\section{Discussion}

Here we want to discuss the pros and cons of the proposed method. First, the method is deterministic - it relies on direct computation of the best match among tens of thousands of precomputed simulations under the measure of interest. Even though we provided a comparison with different query strategies analysing the trade-off between accuracy and time, by using the direct query strategy one obtains "global optimum" for the query at an acceptable runtime, so we suggest to use it in practice. Also despite using the DICE score as such a measure (since it is the first metric of choice in our community), the method is evidently generalizable to any arbitrary measure. The overall runtimes for the solving the inverse problem, including registration, query, and back transform vary from two (for the inaccurate RF query) to eight minutes (for the optimal direct query). This is an order of magnitude faster than any existing brain tumor inverse problem solving methods (starting from ca. 1 hour for the same resolution by \cite{subramanian2020multiatlas,scheufele2019coupling}). One may question whether such a speed-up comes at the cost of a significant error, which in turn comes from the registration to and finding the best matching in the atlas space. While we would agree with such a statement, it is worth to mention that dominant majority of existing literature resorts to such atlas registration in order to approximate the patient specific anatomy.


\section{Appendix}

\input{sections/VPUNN_07_Appendix.tex}

\end{document}

\section{Introduction}

\subsection{Overview}

This paper describes a new method of accurately modeling the decision costs within a neural network compiler (henceforth referred to as \lq a compiler\rq ) using a neural network trained on task-level profiling recordings. Specifically, this paper's research describes this new technique as applied to Intel's range of Versatile Processor Units (VPUs) and their corresponding compiler, however, the concepts presented herein can be applied to other compilers, processor types, or other applications that also make use of a cost model (e.g. performance prediction).

The first section of this paper describes existing methods for decision tree resolution and other relevant background information for understanding the context of the research reported in this paper. The following section contains a technical detailing of the VPUNN cost model (and variants thereof). Results gathered during the development are presented and discussed in the penultimate section. Finally, the paper is concluded with a brief discussion on the overall merits of the technique and the possible future direction of this work.

\subsection{Background Information}
\label{background_info}

At its most basic, compiling a neural network for execution on the VPU is a matter of scheduling smaller work items for the processor to complete in an efficient order. This is a common scenario for compilers, independent of research field.

\begin{figure}
  \centering
  \includegraphics[width=0.7\textwidth, bb=-40 0 285 248]{dog.png}

  \caption{A sample scheduling of neural network tasks.}
  \label{dog_image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.5\textwidth, bb=-20 0 204 63]{tensors.png}
  \caption{Splitting a tensor on three different axes.}
  \label{split_image}
\end{figure}

Figure \ref{dog_image} shows an example of a hypothetical neural network that classifies an image as either containing a dog or a cat. To execute this network on the VPU, the network is compiled into a queue of low-level tasks such as memory transfers and discrete hardware tensor operations. Scheduling some of these operations can be straightforward --- Convolution A can be completed in a single transfer --- but Convolution B has been split into multiple parts to comply with a limitation of the target platform's memory capacity. This limitation might lead to the compiler subdividing the work items into more compatible sizes. Figure \ref{split_image} shows three options for how this division might be done. In this simple depiction of a 3-dimensional tensor, deciding which axis of division (or \lq split\rq) to choose could be determined by a simple benchmarking of all three options. However, in a real-world scenario where there are more granular division options, optimization toggles, and other significant choices involved, the search space grows rapidly. This is especially true when decisions about multiple adjacent operations are highly correlated with each other. In that case, the inherently combinatorial nature of the problem makes it intractable by brute-forcing all possible decisions.

\begin{figure}
  \centering
  \includegraphics[width=\textwidth, bb=0 0 781 229]{sample_decision_diagram.png}
  \caption{Example of the increasing complexity of a task's decision tree scaling with the amount of options available }
  \label{decision_tree}
\end{figure}

Figure \ref{decision_tree} shows a set of these choices and the resultant decision tree. It can be seen how each additional configurable property markedly increases the number of possible path choices. Because the VPU is an embedded processor with significantly restricted resources (to achieve a low power consumption), the execution schedule can have a very high number of tasks compared to a processor with larger capabilities (e.g. a GPU). This leads to an increased number of compiler decisions to determine how to divide a tensor into transferrable chunks, how much work each distributed unit might consume, configuration of different optimizations (e.g. data compression/decompression, sparse acceleration), etc. all compound to create an exploration space which is unwieldy to resolve in a timely manner.

\subsection{Existing Methods}
\label{existing_methods}

This is not a problem faced solely by VPU compilers and indeed, there is a reasonable body of research dedicated to alleviating the problem. These solutions have primarily taken two forms:

\begin{itemize}
  \item Rule-Based decision: Determining the best scheduling decision by following a series of rules
  \item Cost-based decision: Selecting an outcome based on an algorithm or cost model of the underlying hardware that roughly models the behavior in terms of cycles/power of the target device
\end{itemize}

\subsubsection{Rule Based Methodologies}

Some projects with a decision tree may choose their path based on a series of rules (e.g \citet{Hajek1986KNOWLEDGEBASEDES} describes a rule-based decision system for the best method to maintain a given pavement). An example of a rule might be "Choose to split a tensor on the x-axis if the x-axis is larger than 75 elements". The compiler will test each rule as required until there is only one possible configuration left to choose from. This approach lends itself well to decisions where the choices of schedules are limited and those choices have clear positive and negative effects on the resultant schedule. However, this method becomes difficult to manage when the number of parameters is high, or the effect of each parameter is not significantly pronounced to allow for elimination-style decisions.

\subsubsection{Algorithmic Cost Modeling}
A more advanced method is for compilers to calculate the final configuration through an algorithmic, or mathematical, methodology. The CoRA Tensor Compiler (as described in \citet{CoRA}) for example operates by computing the schedule which minimizes the amount of memory reloading --- on the grounds that this will minimize wasted computation. A similar modeling function was used as the original cost model for the VPU compiler before VPUNN was integrated. However, this method is only as effective as the correlation of the algorithm to the actual performance of a system.

\subsubsection{Performance Prediction}
Performance prediction is a closely related topic to cost modeling. If the compiler knows the predicted performance of two different decisions, it can pick the preferrable one to proceed with.

\citet{Qi2017PaleoAP} provide algorithmic modeling of a hardware target's performance based upon the arithmetic complexity of the network. Their proposal is a generic approach across multiple hardware targets, though it requires to be provided some high-level information about the target's capabilities and processing paradigms.

More recent work leverages neural networks to model performance. \citet{PredictCost} discuss a prediction system operating on a range of Nvidia GPUs. Similarly in 2018, Intel Habana (\citet{HabanaPaper}) trained a neural network to predict the performance of convolutions and fully-connected layers on the Gaudi processor. However, to the authors' knowledge, the use of these performance predictors as a compiler cost model has not yet been published.

The VPUNN system is differentiated from prior work in several additional aspects

\begin{itemize}
  \item the neural network is trained at a hardware \lq task\rq \  level, rather than at a network \lq operation\rq\  level, allowing for more fine-grained modeling
  \item the VPUNN model includes hardware-specific parameters particular to Intel VPUs
  \item the model supports all available hardware-native operations of the VPU accelerator, rather than a subset.
  \item the VPUNN model is a product-grade cost model for the VPU compiler
\end{itemize}

VPUNN determines the most suitable selection of compiler decisions based on the projected cycle time as computed by its neural network. How this is achieved is described in Section \ref{tech}.

As of the time of writing, the VPUNN system is the cost modeling technology of choice for Intel VPU neural network compilers.

\section{Technical Overview}
\label{tech}

\begin{figure}
  \centering
  \includegraphics[width=0.8\textwidth, bb=-50 0 327 319]{phases.png}
  \caption{The two phases of VPUNN integration}
  \label{phases_fig}
\end{figure}

A large dataset is key to the successful training of a neural network. Data for training the neural network behind VPUNN was captured using two different methodologies:
\begin{enumerate}
  \item Compiling real-world neural networks into VPU tasks. Each network can generate several hundred tasks, allowing for quick initialization of data points into the database
  \item Creating a wide range of processor tasks through a generation suite. This has been particularly important to unbias the database, as it allows to balance all of the operations and datatype classes.
\end{enumerate}

The generated tasks are then profiled to collect their performance. It is possible to profile individual workloads with very hard granularity using both hardware emulators (Register-transfer level (RTL) simulation, Field-programmable gate array (FPGA) emulation, VPU behavioral model) and directly from the final silicon. This flexibility allows to have a very accurate cost model of the hardware well before silicon tape out, greatly increasing and simplifying the software and compiler development effort. Having such an accurate cost model of the hardware was very important to drive hardware/software co-design decisions at an early stage of product development.

Finally, task descriptions and cycle times are formatted to conform to a prescribed schema and submitted as records to a database. The top half of figure \ref{phases_fig} provides a high-level view of this workflow.

Once the database has enough submissions for an initial deployment of the cost model, the values are downloaded to a single location, and using the TensorFlow framework, the training phase of a neural network is initiated. This trained neural network can be versioned and released for consumption. VPUNN provides two different versions of the cost model that can be used by a VPU NN compiler. These two versions of VPUNN allow a compiler to choose between accuracy and speed depending on the preference of the host compiler and the application requirements. Both model versions expose the same Application Programming Interface (API) so the choice is completely transparent to the end-user of the library. Figure \ref{phases_fig} shows how the model can be deployed and used by the VPU compiler.

\subsection{Data preprocessing and feature engineering}

Initially, there are 26 raw input features which are then preprocessed and transformed into the final 71 features. The preprocessing step consists of categorically encoding the features that do not represent values, and logarithm the remaining numeric features for normalization. A more detailed analysis of the categorical encoding can be seen in Table \ref{categorical_table}. After these transformations have been done, the outlier data points are removed from the dataset before doing the train-validation split.

The output data is also transformed: The network is not trained to predict cycles values, which take values in the 100--10M range, but the utilization rate, which takes values in the 0--1 interval for dense workloads, a value that is much easier for the network to learn. The computed utilization rate can then be used to obtain the cycle value by understanding the cycles in a perfect utilization rate scenario.

\subsection{VPUNN model variants}
\subsubsection{Neural Network using Multi Layer Perceptron (MLP)}

For scenarios where a compiler may require a low-latency cost model, it can choose the VPUNN distribution based on an MLP network topology. This neural network has an input for each relevant aspect of a hardware task's description, has trained its parameters based on the database entries, and predicts cycle cost as an output. A simplified diagram can be seen in Figure \ref{Matrix Mul Diagram}(i). In addition to other performance-enabling features such as caching, this method can save a lot of time for a compiler - considering that the cost model may need to be queried tens of thousands of times for a single network compilation (due to a large number of tasks and configuration options as described in \ref{background_info}). Here, VPUNN takes advantage of the property of MLP network topologies to be converted to matrix-multiplication equivalents which can often take advantage of processor SIMD or vector instructions (\citet{conv_as_matmul}).

\begin{figure}
  \centering
  \includegraphics[width=0.9\columnwidth, bb=-60 0 479 223]{methods.png}
  \caption{High level concept: (i)MLP method (ii) Embedding method}
  \label{Matrix Mul Diagram}
\end{figure}

\subsubsection{Neural Network using Embeddings}

During training, the second style of the VPUNN model embeds costs from the database in an N-dimensional space, where each dimension corresponds with a property of the dataset (e.g., output channels). A simplified depiction of this insertion can be seen in Figure \ref{Matrix Mul Diagram}(ii).

When an inference of the neural network occurs, it will either return an exact value of a previous recording or interpolate between the nearest neighboring values to return an approximation. This technique is inspired by the similarity searches used by \citet{Guo2020AcceleratingLI} and \citet{johnson}. This style of deployment provides higher confidence results due to its ability to return an exact match (where one exists), as opposed to the MLP  method which can only return predictions.

\begin{figure}
\centering
\subfloat{{\includegraphics[width=0.15\columnwidth,height=6cm, bb=0 0 174 764]{matrix_multiplication_architecture.png} }}
\hspace{0.1\linewidth}
\subfloat {{\includegraphics[width=0.15\columnwidth, height=5cm, bb=0 0 174 572]{embeddings_architecture.png} }}
\caption{Network Architectures: (i)MLP method (ii) Embedding method}
\label{architectures}
\end{figure}

\section{Results \& Discussion}

\subsection{Training VPUNN to approximate the hardware model}
\label{hyperparams}

The VPUNN model analyzed in this paper was trained on 31366 recordings of simulated workloads, using 80\% of the data for training and 20\% for validation.
The network was trained on an NVIDIA A100-SXM4 40GB GPU. On this system the training process takes 36 minutes - The network quickly learns features in the first few iterations and slowly refines them until we no longer see a significant benefit. The MAPE curve for the network can be seen in Figure \ref{AC}.

The Early Stopping callback is used, with a patience of 100 epochs and a minimum delta of 0.1. The network is trained with a batch size of 64, for 1000 epochs, however, it only requires approximately 700 epochs to fully train, when the training stops because of the callback. These hyperparameters were obtained through repeated experiments and proved themselves to offer the best validation results.

At the end of the training, our model has a median validation Absolute Percentage Error (APE) of 0.72\%. This signals a high rate of confidence in the modeling approximation, and the resulted performance has high reproducibility throughout multiple experiments with different seeds as seen in Figure \ref{experiments_table}.

The Embedding version of VPUNN has been trained in the same fashion, however, there is a significant improvement in results. The Embedding model has a median validation APE of 0.31\%. As expected, this variant performs better.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth, bb=-350 0 2163 833]{epoch_loss_diagram.png}
  \caption{MAPE curve for matrix multiplication network}
\label{AC}
\end{figure}

\subsection{Effectiveness as a cost model}
\label{EF}

The fact that the VPUNN model is comparable to the hardware modeling does not necessarily imply a guarantee of improved cost modeling and resultant task schedules. For instance, a poor hardware modeling may be an inferior cost model when compared to a mathematical formula that better correlates to the costs in the system. The only way to verify whether it is a better cost model is by direct comparison of the VPU's performance when using VPUNN in place of the existing, rule-based, methodology. Figure \ref{vpux_table} shows the improvements in VPU execution time for a range of networks (in Frames per Second (FPS)) due to the integration of the VPUNN cost model inside the compiler.

The findings are significant in that every neural network that was profiled either improved or did not degrade. This means that the cost model was capable of guiding the compiler to make better decisions than the previous rule-based approach.

\pgfplotstableread[row sep=\\,col sep=&]{
  interval & amount & carD & carR \\
  0\%     & 5  & 0.1  & 0.2  \\
  3\%     & 6 & 3.8  & 4.9  \\
  5\%    & 6 & 10.4 & 13.4 \\
  8\%   & 4 & 17.3 & 22.2 \\
  10\%   & 1  & 21.1 & 27.0 \\
  13\%      & 2  & 22.3 & 28.6 \\
  15\%      & 2  & 22.3 & 28.6 \\
  18\%      & 0  & 22.3 & 28.6 \\
  20\%      & 2  & 22.3 & 28.6 \\
  23\%      & 1  & 22.3 & 28.6 \\
  25\%      & 0  & 22.3 & 28.6 \\
}\mydata

\begin{figure}
  \begin{tikzpicture}
    \begin{axis}[
        ybar,
        symbolic x coords={0\%,3\%,5\%,8\%,10\%,13\%,15\%,18\%,20\%,23\%,25\%},
        xtick=data,
        width=\textwidth,
        height=5cm,
      ]
      \addplot table[x=interval,y=amount]{\mydata};
    \end{axis}
  \end{tikzpicture}
  \caption{Observed improvements (x) in a quantity of neural networks(y)  }
  \label{vpux_table}

\end{figure}

\subsection{Execution Time}
\label{ET}

One might have been able to achieve similar results if they had used the complex hardware modeling directly as a cost model. However, these models are often very time-consuming and would be too slow for a standard user of the VPU compiler.

Table \ref{ms_table} shows the time to query both the complex hardware model and VPUNN for the different costs through a selection of networks. It can be seen that VPUNN takes far less time to return a result. Without this level of performance, using VPUNN as a cost model would have limited application, perhaps used as an optimization toggle to the compiler, but it is fast enough to be considered as a standard cost model for a compiler.

\begin{table}
  \caption{VPUNN execution times versus hardware model}
  \label{ms_table}
  \centering
  \begin{tabular}{lll}
    \toprule
    Network      & MLP VPUNN & Embeddings VPUNN \\
    \midrule
    ResNet-50    & 28.9x   & 21.1x              \\
    MobileNet v2 & 65.9x   & 31.2x              \\
    Yolo Tiny v2 & 47.8x   & 31.5x              \\
    \bottomrule
   \end{tabular}
\end{table}

Because of the widespread improvements were seen in Sections \ref{EF} and \ref{ET}, VPUNN has since been established as the primary cost model for all future VPU products.

\subsection{Limitations}
\label{limits}

There are two main limitations to the VPUNN approach to cost modeling. Firstly, VPUNN cannot be more accurate than the modeling of hardware it is approximating. The VPU target used in this paper was the third generation VPU in which the gathering of task-level profilings on hardware doesn't allow a zero observer effect. This made the profiling of individual workloads a difficult and error-prone task. An accurate modeling of the VPU runtime has been done to mitigate the effect of measurement on task profiling but the whole process is tricky at best to fine-tune.

To overcome these limitations, future VPU generations have a cycle-accurate hardware timer specifically tailored for individual task profiling. We include some preliminary findings in Appendix \ref{A_AB}, strongly indicating that using cycle costs from the actual hardware significantly reduces the impact of this limitation.

The second limitation is that gaps in training data could result in bad cost predictions by VPUNN. It is important, therefore, to have a training data set that covers as much of the network's design space in as much detail as possible. This large problem space needs a quite large training set to enable the network to properly generalize across unseen workloads.

\section{Conclusion and Future Plans / Use}
There is a clear improvement in the performance of the VPU when task schedules are compiled using VPUNN over the existing mathematical methods. VPUNN has been integrated into the OpenVINO VPU plugin because it is a very accurate model of complex hardware without the performance overhead included in a transaction level model.

A significant benefit of VPUNN is its enhanced ability to generalize to new and unseen networks (as discussed in Section \ref{EF}) and to efficiently scale to new workloads without much technical overhead. While an algorithmic model may require handcrafted analysis for every new type of workload, updating a neural network is much simpler and automatable once new workloads are collected in the database. It is possible to iteratively collect the data directly during pre-silicon validation and directly from end-users to allow the model to improve over time.

Though the topic of this paper has been focused on VPU, the concept is generic enough to apply VPUNN's successes to other targets provided the training is updated to appropriate task parameters.

\subsection{Cost Model Tuning}
As the VPU moves through the different stages in its lifecycle - from design to delivery, the data used to train VPUNN can change too. Early-stage simulations or software modelings can be replaced with actual recordings on physical devices. This means that over time, releases of VPUNN for a target will be incrementally more accurate and representative.

If data gathering is straightforward, software modeling may not even be needed. Provided with sufficient data gathering, VPUNN is very reusable between generations of a target, or even between different targets. To update the system, a new set of model parameters may need to be defined - whereas software modeling may need a lot of engineering time to understand the new properties of the system, assert the continued properties of continued features, debug mistakes, and so on. There is also a significant advantage in having a system less prone to the human element --- gathering extra training data is simpler than rewriting an entire software project!

\subsection{Dataset Exploration}
Outside of the improvements in cost model accuracy, having a singular repository of performance data is valuable for other uses for early-stage hardware developers. Exploration within the dataset can unveil unexpected performance characteristics, corner cases, or unmet expectations. Learnings from such investigations can be fed back into different parts of the organization to inform architectural and hardware design decisions, as well as firmware development choices. Collecting these datasets will enable future VPU generations to better understandings about performance and power at a much earlier stage of their lifecycles.

\subsection{Future Work: Neural Architecture Search}
The gathered dataset of tasks for VPUNN can be reused for other machine learning initiatives. One prominent example is using Neural Architecture Search (as described by \citet{ZophL16}) to vary the parameters of neural networks to contain operations and tensors that better suit the hardware whilst not compromising on accuracy. The work by \citet{NAS_TPU} has shown this to be a successful technique for optimizing the power utilization and performance throughput of a given network structure for Google's Tensor Processing Units (TPUs)--- a processor with similar trade-off considerations as the Intel VPUs. \citet{hamza_and_friends} have some established work in this area on a previous version of the VPU.

\subsection{Future Work: Optimizing for different costs}
For performance, the cycle count is an effective metric to train our neural networks for. However, particularly for Intel's range of VPU and other processors in the same market, power can be an important factor in competitive comparisons. Should this work be extended to gather the power cost for each task, a compiler could produce schedules that are optimized for power, rather than speed --- or a hybrid combination.

\paragraph{Summary}
The invention of VPUNN significantly improves the ability of the VPU compiler to produce the most performant task schedule for a target VPU. It provides the benefit of a complex hardware modeling system at a fraction of the execution overhead. Due to these merits, VPUNN is now the standard costing system for the VPU compiler, already integrated into the Intel OpenVINO toolkit.

\section{Acknowledgements}

\section*{Checklist}

\begin{enumerate}

  \item For all authors...
        \begin{enumerate}
          \item Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
                \answerYes{These sections claim novel improvements in VPU compiler scheduling due to the paper and the following sections describe how we achieved this and our results backing up this claim}
          \item Did you describe the limitations of your work?
                \answerYes{Limitations are discussed in Section \ref{limits}}
          \item Did you discuss any potential negative societal impacts of your work?
                \answerNA{This ideas in this paper are purely performance improvements to existing systems}
          \item Have you read the ethics review guidelines and ensured that your paper conforms to them?
                \answerYes{}
        \end{enumerate}

  \item If you are including theoretical results...
        \begin{enumerate}
          \item Did you state the full set of assumptions of all theoretical results?
                \answerNA{results are not theoretical, but actual gathered from hardware}
          \item Did you include complete proofs of all theoretical results?
                \answerNA{results are not theoretical, but actual gathered from hardware}
        \end{enumerate}

  \item If you ran experiments...
        \begin{enumerate}
          \item Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?
                \answerNo{Code is proprietary, however we have provided in the appendix a list of the tests used and users of the VPU compiler should be able to observe these performance improvements on these networks}
          \item Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?
                \answerYes{All the training details are specified in Section \ref{hyperparams}}
          \item Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?
                \answerYes{Multiple results of the same experiment can be observed in Figure \ref{experiments_table}}
          \item Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?
                \answerYes{The resources used are specified in Section \ref{hyperparams}}
        \end{enumerate}

  \item If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
        \begin{enumerate}
          \item If your work uses existing assets, did you cite the creators?
                \answerNA{This work used performance recordings, compiled by the researchers}
          \item Did you mention the license of the assets?
                \answerNA{The assets used are not available outside of Intel.}
          \item Did you include any new assets either in the supplemental material or as a URL?
                \answerNA{}
          \item Did you discuss whether and how consent was obtained from people whose data you're using/curating?
                \answerNA{}
          \item Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?
                \answerNA{}
        \end{enumerate}

  \item If you used crowdsourcing or conducted research with human subjects...
        \begin{enumerate}
          \item Did you include the full text of instructions given to participants and screenshots, if applicable?
                \answerNA{No human subject research was conducted}
          \item Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?
                \answerNA{}
          \item Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?
                \answerNA{}
        \end{enumerate}

\end{enumerate}

\subsection{Future Work: Improved Accuracy }
\label{A_AB}

Comparison of the Absolute Percentage Error between previous modelling tools and VPUNN when trained on hardware recordings --- versus the actual cycle count of hardware recordings. These were gathered on early versions of Generation 4 Intel VPUs, as described in Section \ref{limits}.
These results shows promise of improved accuracy of results if VPUNN were updated to use real hardware timings rather than simulated timings.

\begin{figure}[ht]
  \centering
  \includegraphics[width=0.7\textwidth, bb=-50 0 307 263]{MAPE_AB.png}
  \caption{ArchBench absolute percentage error}
\label{MAPE_AB}
\end{figure}

\begin{figure}[ht]
  \centering
  \includegraphics[width=0.7\textwidth, bb=-50 0 307 263]{MAPE_VPUNN.png}

  \caption{VPUNN absolute percentage Error}
\label{MAPE_VPUNN}
\end{figure}

\subsection{Preprocessing}

Along with the fields presented here, there are some excluded fields set to zero, that are reserved for future versions of VPUs.

\FloatBarrier
\begin{table}[h]
  \caption{Existing values for categorical features}
  \label{categorical_table}
  \centering
  \begin{tabular}{ll}
    \toprule
    Feature name & Values \\
    \midrule
    input datatype   &  [uint8; int8; float16; bfloat16]          \\
    weights datatype   &  [uint8; int8; float16; bfloat16]         \\
    output datatype   &  [uint8; int8; float16; bfloat16]          \\
operation & \multirow{2}{10cm}{[Convolution; CM Convolution; DW Convolution; Elementwise; MaxPool; AveragePool]}\\
&  \\
    execution mode & [Vector; Matrix; Vector FP16; Cuboid 16x16; Cuboid 8x16; Cuboid 4x16] \\
    activation function & [None; RELU; LRELU; Add; Sub; Mult] \\
    \bottomrule
  \end{tabular}
\end{table}
\FloatBarrier

\subsection{Error Bars }

\FloatBarrier
\pgfplotstableread[row sep=\\,col sep=&]{
    sample & median_ape \\
    1     & 0.81  \\
    2     & 0.6   \\
    3     & 0.75  \\
    4     & 0.68  \\
    5     & 0.71  \\
    6     & 0.74  \\
    7     & 0.7   \\
    8     & 0.68  \\
    9     & 0.77  \\
    10    & 0.73  \\
    }\experimentsdata

\begin{figure}[h]
\begin{tikzpicture}
    \begin{axis}[
            ybar,
            xtick=data,
            ymin=0,
            ymax=1,
	width=\textwidth,
	height=5cm,
        ]
        \addplot table[x=sample,y=median_ape]{\experimentsdata};
    \end{axis}
\end{tikzpicture}
  \caption{Median validation error bars (y) with respect to experiments with different seeds (x)}
  \label{experiments_table}
\end{figure}
\FloatBarrier

\subsection{Networks used for Assessment}

Details of the networks used for assessing VPUNN in comparison to the existing cost modelling within the VPU compiler. All networks were found to be either improved or unchanged in performance when VPUNN was used.

\FloatBarrier
\begin{table}[h]
  \caption{Networks analysed for performance improvement/degradation}
  \label{networks_list}
  \centering
  \begin{tabular}{lll}
    \toprule
    Network      & Precision & Hyperlink \\
    \midrule

\detokenize{por_tf_FP16-INT8_yolo-v2-tiny-ava-0001} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_yolo_v2_tiny_ava_0001.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16-INT8_yolo_v4} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/2021.2/omz_models_public_yolo_v4_tf_yolo_v4_tf.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16-INT8_yolo_v4_tiny} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_yolo_v4_tiny_tf.html}{OpenVINO Model Info}  \\

\detokenize{BDK_FP16-INT8_resnet-50-pytorch-sparse50} & Mixed FP16-INT8 & \multirow{2}{3cm}{Internal variant of  \href{https://docs.openvino.ai/latest/omz_models_model_resnet_50_pytorch.html}{OpenVino model}}\\
& &  \\

\detokenize{por_onnx_FP16-INT8_resnet-50-pytorch} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_resnet_50_pytorch.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16-INT8_googlenet-v1} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/2021.2/omz_models_public_googlenet_v1_tf_googlenet_v1_tf.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16_googlenet-v1} & FP16 & \href{https://docs.openvino.ai/2021.2/omz_models_public_googlenet_v1_tf_googlenet_v1_tf.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16-INT8_googlenet-v3} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_googlenet_v3.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16_googlenet-v3} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_googlenet_v3.html}{OpenVINO Model Info} \\
\detokenize{por_caffe2_FP16_squeezenet1.1} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_squeezenet1_1.html}{OpenVINO Model Info} \\
\detokenize{por_caffe2_FP16-INT8_squeezenet1.1} & Mixed FP16-INT8 &  \href{https://docs.openvino.ai/latest/omz_models_model_squeezenet1_1.html}{OpenVINO Model Info} \\
\detokenize{scale_caffe_FP16_mobilenet-v1-1.0-224} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v1_1_0_224_tf.html}{OpenVINO Model Info} \\
\detokenize{scale_caffe_FP16-INT8_mobilenet-v1-1.0-224} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v1_1_0_224_tf.html}{OpenVINO Model Info} \\
\detokenize{por_onnx_FP16-INT8_mobilenet-v2} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v2.html}{OpenVINO Model Info} \\
\detokenize{ea_caffe_FP16_age-gender-recognition-retail-0013_ww14} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_age_gender_recognition_retail_0013.html}{OpenVINO Model Info} \\
\detokenize{scale_caffe_FP16-INT8_age-gender-recognition-retail-0013} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_age_gender_recognition_retail_0013.html}{OpenVINO Model Info} \\
\detokenize{scale_onnx_FP16_resnet-18-pytorch} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_resnet_18_pytorch.html}{OpenVINO Model Info} \\
\detokenize{scale_onnx_FP16-INT8_resnet-18-pytorch} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_resnet_18_pytorch.html}{OpenVINO Model Info} \\
\detokenize{scale_caffe_FP16_mobilenet-v2} &  FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v2.html}{OpenVINO Model Info} \\
\detokenize{scale_caffe_FP16-INT8_mobilenet-v2} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v2.html}{OpenVINO Model Info} \\
\detokenize{FP16-INT8_InternalModel1} & Mixed FP16-INT8 & Internal Model \\
\detokenize{INT8_InternalModel2_V2} & INT8 & Internal Model \\
\detokenize{INT8_InternalMode2_V3} & INT8 & Internal Model\\
\detokenize{scale_tf_FP16_mobilenet-v3-large-1.0-224} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v3_large_1_0_224_tf.html}{OpenVINO Model Info} \\
\detokenize{scale_tf_FP16-INT8_mobilenet-v3-large-1.0-224} & Mixed FP16-INT8 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v3_large_1_0_224_tf.html}{OpenVINO Model Info} \\
\detokenize{por_tf_FP16_mobilenet-v3-small-1.0-224} & FP16 & \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v3_small_1_0_224_tf.html}{OpenVINO Model Info}  \\
\detokenize{por_tf_FP16-INT8_mobilenet-v3-small-1.0-224} & Mixed FP16-INT8 &  \href{https://docs.openvino.ai/latest/omz_models_model_mobilenet_v3_small_1_0_224_tf.html}{OpenVINO Model Info} \\
\detokenize{FP32_InternalModel3} & FP32 & Internal Model \\

    \bottomrule
  \end{tabular}
\end{table}
\FloatBarrier
\textbf{Note:} The networks referenced here may be using internal releases of OpenVINO or be non-published variants of the published network.

\textbf{Note:} weblinks accessed 2022-04-29


\section{Attention Flow}
\label{app:attention_flow}

Despite \emph{attention flow}, Grad-SHAP, and Deep-SHAP all (supposedly) being valid Shapley Value explanations \cite{ethayarajh-jurafsky-2021-attention}, Table \ref{tab:attn_flow_apx} shows that the agreement is low.

\input{tables/appendix/attn_flow}

\section{Result tables with standard deviation}
\label{app:results-with-std-dev}

In Tables \ref{app:distilbert} and \ref{app:BILSTM} we provide the same results as in Tables \ref{tab:distilbert_results} and \ref{tab:bilstm_results} of the main paper, but with the standard deviation for each experiment.

\input{tables/appendix/distilbert_correlation_table}
\input{tables/appendix/bilstm_agreement_table}

\section{Model Performance}
\label{app:performance}

We include a uniform activation baseline to contextualize the attention mechanism's utility. Table \ref{tab:accuracy} notes the gap between the performance of uniform and softmax attention in the BiLSTM is never higher than 1-2\%. This distinction aligns with the results of \citet{wiegreffe-pinter-2019-attention} and \citet{vashishth2019attention}, who argue claims of attention-based interpretability are stronger in situations in which models need the attention module to solve the underlying task. Of course, it is difficult to prove a causal effect in a deep network: the BiLSTM may solve the task differently depending on whether or not the attention mechanism is ablated, meaning it is still possible to make claims of interpretability whether or not softmax attention leads to higher task performance. We again emphasize that we do not wish to take a side in the ``attention explanation" debate, but this point further stresses the difficulty of proving anything with the \emph{agreement as evaluation} paradigm.

\input{tables/appendix/accuracy}

\section{Reproducibility}
\label{app:reproducibility}

Our code is publicly available at \url{https://github.com/sfschouten/court-of-xai}.
We conducted our experiments on Amazon Web Services \texttt{g4dn.xlarge} EC2 instances using an NVIDIA T4 GPU with 16GB of RAM. The version of PyTorch was \texttt{1.6.0+cu101}.
We refer to Table \ref{tab:runtime} for the average time to train each model on each dataset.

The DistilBERT model contained 66955779 trainable parameters and the BiLSTM model contained 12553519 trainable parameters, as reported by the AllenNLP library \cite{gardner-etal-2018-allennlp}.
Table \ref{tab:dataset_statistics} lists the number of instances in each split of each dataset and Table \ref{tab:val_accuracy} details the accuracy of our models on the validation sets during training.

Links to download versions of all datasets are included in our code repository. For posterity, links to all datasets are listed here:
\begin{itemize}
	\item \textbf{SST-2}: \url{https://github.com/successar/AttentionExplanation/tree/master/preprocess/SST}
	\item \textbf{IMDb}: \url{https://github.com/successar/AttentionExplanation/tree/master/preprocess/IMDB}
	\item \textbf{SNLI}: \url{https://nlp.stanford.edu/projects/snli/}
	\item \textbf{MNLI}: \url{https://cims.nyu.edu/~sbowman/multinli/}
	\item \textbf{XNLI}: \url{https://cims.nyu.edu/~sbowman/xnli/}
	\item \textbf{Quora Question Pair}: \url{https://drive.google.com/file/d/12b-cq6D45U5c-McPoq2wsFjzs6QduY_y/view?usp=sharing}
\end{itemize}

\input{tables/appendix/runtime}
\input{tables/appendix/dataset_statistics}
\input{tables/appendix/validation_accuracy}

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

As machine learning (ML) models are increasingly used in hybrid settings to make consequential decisions for humans, criteria for plausible and faithful explanations of their predictions remain speculative \citep{lipton2016mythos, jacovi-goldberg-2020-towards}.
Although there are many possible explanations for a model's decision, only those faithful to both the model's reasoning process and to human stakeholders are desirable \cite{jacovi-goldberg-2020-towards}.
The rest are irrelevant in the best case and harmful in the worst, particularly in critical domains such as law \cite{Kehl2017AlgorithmsIT}, finance \cite{grath_interpretable_2018}, and medicine \cite{caruana2015intelligble}.

Content moderation is a use case where explanations can help domain experts as part of a hybrid human-in-the-loop system.
Consider an ML model that predicts whether or not a post on a social media contains misinformation: when a post is automatically removed by the model and the user who created it appeals its removal, content moderators need to read through the entire post to identify why it was flagged as containing misinformation \cite{halevy_preserving_2022}.
Including explanations that identify which parts of the post are problematic can help content moderators decide if the model was correctly flagging the post.

ML practitioners frequently explain models by calculating each input's contribution toward an individual prediction.
Additivity --- treating all contributions as independent and quantifiable --- is a common simplifying assumption.
In this work, we focus on such additive explanation methods and refer to them as \emph{feature attribution methods}.
We denote the contribution of each input toward the model's decision as its \emph{importance}. We say that two different feature attribution methods \emph{agree} if there is a strong correlation between their computed rankings of input importance.

The attention mechanism \cite{bahdanau2014neural} in natural language processing (NLP) is a popular, albeit less rigorously motivated, way of obtaining explanations. Because the mechanism produces context vectors from which decoders can soft-search for prediction-relevant information, the weights assigned to inputs intuitively serve as proxies for their overall contribution towards a decision. Weights are often visualized as heatmaps over sequences \cite[e.g.,][]{pmlr-v37-xuc15,li-etal-2016-visualizing}, which can be particularly persuasive when examples are (unintentionally) cherry-picked to fit a narrative.
We define an attention-based explanation as a vector of attention weights that, similar to a feature attribution method, can be treated as a ranking of importance.

In their critique of attention-based explanations for NLP, \citet{jain-wallace-2019-attention} argue that faithful attention-based explanations must be highly \emph{agreeable}.\footnote{\citet{ethayarajh-jurafsky-2021-attention} use the term \emph{consistent}.} That is, their generated rankings of input importance must correlate with existing feature attribution methods.
Following \citet{jain-wallace-2019-attention} and their claim that ``attention is not explanation'', several recent papers have presented an increased agreement with a small set of feature attribution methods as evidence for their proposed method's ability to improve the faithfulness of the attention mechanism. For example, \citet{mohankumar-etal-2020-towards} show that minimizing hidden state conicity in a BiLSTM improves the Pearson correlation of attention weights with Integrated Gradients \cite{intgrad} attributions. As the popularity of \emph{agreement as evaluation} grows \citep{meister-etal-2021-sparse,abnar-zuidema-2020-quantifying}, we believe it is important to investigate diagnostic capacity of agreement as a metric by examining (i) more recent feature attribution methods, and (ii) more complex transformer-based models.

Under the paradigm of \emph{agreement as evaluation}, new explanation methods (i.e., attention-based) are compared to established explanation method(s) (i.e., feature attributions).
However, can any one explanation method act as the standard against which other explanation methods are evaluated?
Explanations are task-, model-, and context-specific \cite{doshi-velez_towards_2017}, and the performance of explanation methods depends on the particular diagnostic tests considered \citep{deyoung-etal-2020-eraser,Robnik-sikonja2018}.

In this work, our main research question is: \emph{How well do attention-based explanations correlate with recent feature attribution methods for NLP tasks?} Specifically, we want to investigate:
\begin{itemize}
	\item \textbf{RQ1:} Does the correlation depend on the model architecture (transformer- vs. LSTM-based)?
	\item \textbf{RQ2:} Does the correlation depend on the nature of the classification task (single- and pair-sequence)?
\end{itemize}

We investigate the following feature attribution methods: {LIME} \cite{lime}, {Integrated Gradients} \cite{intgrad}, {DeepLIFT} \cite{deeplift}, and two versions of SHAP \cite{shap}:
{Grad-SHAP} (based on Integrated Gradients) and {Deep-SHAP} (based on DeepLIFT).
We observe low agreement between attention-based explanations and feature attribution methods, across both models and both tasks.
We also observe low agreement across all explanation methods for the transformer-based model and for pair-sequence tasks.
We use this empirical evidence, along with our theoretical objections, to argue that practitioners should refrain from evaluating attention-based explanations based on their agreement with feature attribution methods.

\section{Related Work}

\citet{jain-wallace-2019-attention} introduced the \emph{agreement as evaluation} paradigm by comparing attention-based explanations with simple feature attribution methods using a bidirectional LSTM \cite{lstm} on single- and pair-sequence classification tasks.
They conclude that ``attention is not explanation'' due to the weak correlation between the rankings of input token importance obtained from attention weights and those obtained from two elementary feature attribution methods: (i) input $\times$ gradient \citep{inputxgradients1, inputxgradients2}, and (ii) leave-one-out \cite{li2017understanding}.
In our work, we test the generalizability of agreement as an evaluation metric by (i) testing on a more comprehensive set of feature attribution methods (see Section~\ref{section:xai-methods}), and (ii) testing on a transformer-based \cite{vaswani2017attention} model on the same types of tasks (see Section~\ref{section:models}).

The influential critique by \citet{jain-wallace-2019-attention} has sparked an ongoing debate about whether or not attention is explanation \citep{wiegreffe-pinter-2019-attention, serrano-smith-2019-attention, grimsley-etal-2020-attention}.
More recently, \citet{bastings-filippova-2020-elephant} have questioned the notion of ``attention as explanation'' as a whole, and suggest that in order to explain ML model predictions, the community should rely on methods that are explicitly created for this purpose (i.e., feature attributions), instead of seeking explanations from attention mechanisms.
In our work, we do not aim to take a position on the ``is attention explanation'' debate, but rather investigate the hypothesis that in order for attention mechanisms to be considered ``explanations'', they must correlate with existing feature attribution methods. This is an underlying assumption of not only the work by \citet{jain-wallace-2019-attention}, but also that of \citet{meister-etal-2021-sparse}, who show that inducing sparsity in the attention distribution decreases agreement with feature attribution methods, and \citet{abnar-zuidema-2020-quantifying}, who demonstrate their \emph{attention-flow} algorithm improves the correlation with attributions based on feature ablation.

\citet{atanasova-etal-2020-diagnostic} introduce a series of diagnostic tests to evaluate feature attribution methods for text classification.
They show that the performance of feature attribution methods, measured by using these diagnostic tests, largely depends on the model and task considered, but note that gradient-based methods tend to perform the best.
Similar to their work, we also compare and evaluate feature attribution methods, but only to investigate the suitability of \textit{agreement as evaluation}, not to determine a winning explanation method given several diagnostic tests.
In contrast, we evaluate five feature attribution methods on agreement as evaluation based on only the rank correlation using Kendall's-$\tau$, following \citet{jain-wallace-2019-attention}.

\citet{prasad-etal-2021-extent} define three alignment metrics that quantify how well human-annotated natural language explanations align with the explanations generated by the Integrated Gradients method \citep{intgrad}.
They find that the BERT \cite{devlin-etal-2019-bert} model has the highest alignment with human-annotated explanations.
However, unlike our work, this work focuses exclusively on using transformer-based models and Integrated Gradients, and does not question the \emph{agreement as evaluation} paradigm for feature attribution methods.

\citet{ding-koehn-2021-evaluating} introduce a human-annotated benchmark to evaluate feature attribution methods for NLP models.
They test two attribution methods on three types of NLP models and find that explanations from feature attribution methods are sensitive to changes in model configuration.
Similarly, we test five attribution methods on two types of NLP models, but our focus is on investigating the \emph{agreement as evaluation} paradigm, whereas \citet{ding-koehn-2021-evaluating} focus on investigating the correlation between feature attribution methods and a human-annotated benchmark.
We are unable to incorporate the benchmark proposed by \citet{ding-koehn-2021-evaluating} in our work because the human annotations are not present for every single token, and therefore we cannot turn them into a ranking.

This work builds off our prior ICML workshop paper \cite{neely2021order}, which several other papers have extended. \citet{feldhus-etal-2021-thermostat} introduce a software package to analyze instance-wise explanations for popular NLP models and tasks and, in doing so, partially reproduce one of our experiments. \citet{krishna_disagreement_2022} formally define and highlight the importance of the "disagreement problem" between feature attribution methods, which they find is a constant frustration for ML practitioners. They introduce metrics to capture the disagreement between top-k features and ``features of interest'' (e.g., those selected by an end-user) and find considerable disagreement between feature attribution methods for tabular, text, and image data on both real-world and research datasets.

\section{Explainability in NLP} \label{sec:method}
In this work we investigate two types of explanations for NLP tasks: (i) those from recent feature attribution methods in the \ac{XAI} literature and (ii) those based on attention scores. We evaluate on transformer- and LSTM-based models (\textbf{RQ1}), for both single- and pair-sequence tasks (\textbf{RQ2}).
Following \citet{jain-wallace-2019-attention}, we define an \emph{explanation} of an input token sequence as a vector comprised of an importance score for each token.
This vector can be used to rank the tokens by their importance score.

\subsection{Explanations from Feature Attributions} \label{ssec:recent_methods}
\label{section:xai-methods}

In our experiments, we focus on five recent feature attribution methods: {LIME} \cite{lime}, {Integrated Gradients} \cite{intgrad}, {DeepLIFT} \cite{deeplift}, and two versions of SHAP \cite{shap}:
{Grad-SHAP} and {Deep-SHAP}:

\begin{itemize}
	\item LIME \cite{lime} produces locally faithful explanations by learning an interpretable (e.g., linear) model from samples weighted by their proximity to the original instance.
	\item Integrated Gradients \cite{intgrad} calculates input feature attributions by accumulating the gradients obtained from the model along the straight-line path from a baseline to the original instance.
	\item DeepLIFT \cite{deeplift} also produces input feature attributions using the gradients, but it assigns scores based on the difference between a reference activation and the activation of the original instance. This allows the calculated contributions to remain non-zero even when the gradients are zero.
	\item SHAP \cite{shap} identifies a unique solution for the contribution of each input toward the prediction. Since this is computationally expensive, \citet{shap} propose approximations based on existing methods:
	\begin{itemize}
		\item[$\circ$] Grad-SHAP (based on Integrated Gradients).
		\item[$\circ$] Deep-SHAP (based on DeepLIFT).
	\end{itemize}
\end{itemize}

\subsection{Explanations from Attention Mechanisms} \label{ssec:attn_methods}

Given an input sequence of tokens $\textbf{S} = (t_{1}, ..., t_{n})$, we define an \emph{attention-based explanation} as an assignment of attention weights $\boldsymbol{\alpha} \in \mathbb{R}^{n}$ to the tokens in $\textbf{S}$ \citep{jain-wallace-2019-attention}.
Since the dimensionality of $\boldsymbol{\alpha}$ is architecture-dependent, it may be necessary to filter or aggregate the weights. In our experiments, this is only relevant for the self-attention mechanism in the transformer-based model we consider below (see Section~\ref{section:models}).

Previous analyses at the attention head level implicitly assume that contextual word embeddings remain tied to their corresponding tokens across self-attention layers \citep{baan2019understanding, clark-etal-2019-bert}.
This assumption may not hold in transformer-based models since information mixes across layers \cite{brunner2019identifiability}.
Therefore, we use the \emph{attention rollout} \cite{abnar-zuidema-2020-quantifying} method --- which assumes the identities of tokens are linearly combined through the self-attention layers based exclusively on attention weights --- to calculate post-hoc, token-level importance scores.%
\footnote{We also experimented with \textit{attention flow} \cite{abnar-zuidema-2020-quantifying}, see Appendix \ref{app:attention_flow}.}
Following \citet{abnar-zuidema-2020-quantifying}, we use the scores calculated for the last layer's [CLS] token, resulting in a final vector $\boldsymbol{\alpha} \in \mathbb{R}^{n}$ at the time of evaluation.

Recurrent models similarly suffer from issues of identifiability. In LSTM-based models, attention is computed over hidden representations across timesteps, which does not provide faithful token-level importance scores. Approaches that trace explanations back to individual timesteps \cite{bento_timeshap_2021} or input tokens \cite{tutek-snajder-2020-staying} are only just emerging.
Therefore, we analyze the raw attention weights for the LSTM-based model we consider below (see Section~\ref{section:models}).

\subsection{Agreement between Explanations}

Following \citet{jain-wallace-2019-attention}, we measure \emph{agreement} between the explanation methods as the mean Kendall-$\tau$ correlation \cite{kendall-tau} between the ranked importance scores of all input tokens, across all examples.
\footnote{We also calculated Spearman and Pearson correlations and obtained similar results.}
The Kendall-$\tau$ correlation is a widely used metric for comparing ranked lists; it measures the correlation between two ranked lists based on discordant pairs between the lists. Two items are considered \emph{discordant} if they are ranked differently on the two lists.
The Kendall-$\tau$ correlation can take on values in the $[-1, 1]$ interval, where negative values imply the rankings are negatively correlated while positive values imply the rankings are positively correlated \cite{kendall-tau}.

\section{Experimental Setup} \label{sec:experiments}

\subsection{Datasets} \label{sec:datasets}

We evaluate on two types of NLP classification tasks:
\begin{enumerate*}[label=(\roman*)]
	\item single-sequence, and
	\item pair-sequence.
\end{enumerate*}
For the single-sequence task, we perform binary sentiment classification on the Stanford Sentiment Treebank (SST-2) \cite{socher-etal-2013-parsing} and the IMDb Large Movie Reviews Corpus (IMDb) \cite{maas-etal-2011-learning}.
We use identical splits and pre-processing as \citet{jain-wallace-2019-attention}, but also remove sequences longer than 240 tokens for faster attribution calculation. This leaves us with roughly 78\% of the original instances in the IMDb dataset.

\begin{itemize}
	\item The SST-2 dataset consists of single sentences extracted from movie reviews and is used for binary sentiment analysis: predicting whether the review is positive or negative.
	\item The IMDb dataset consists of movie reviews as well, but contains longer sequences compared to SST-2. It is also used for binary sentiment analysis.
\end{itemize}

For the pair-sequence task, we examine natural language inference and understanding with the Multi Natural Language Inference (MNLI) corpus \cite{williams-etal-2020-predicting}, the Stanford Natural Language Inference (SNLI) corpus \cite{bowman-etal-2015-large}, and the Quora Question Pairs dataset \cite{quora}.

\begin{itemize}
	\item The MNLI dataset contains sentence pairs for a textual entailment task: given a pair of sentences, we want to predict whether or not one sentence implies the other. Since MNLI has no publicly available test set, we use the English subset of the XNLI \cite{conneau2018xnli} test set.
	\item The SNLI dataset~\cite{bowman-etal-2015-large} contains pairs of sentences and is used for the textual entailment task, similar to MNLI.
	\item The Quora Question Pairs dataset \cite{quora} contains pairs of questions from the Quora website, where the task is to classify if the two questions in a pair are duplicates or not. We use a custom split (80/10/10) for the Quora dataset, removing pairs with a combined count of 200 or more tokens (leaving 99.99\% of the original instances).
\end{itemize}

\subsection{Models}
\label{section:models}

We test two types of NLP models:
\begin{enumerate*}[label=(\roman*)]
	\item transformer-based, and
	\item LSTM-based.
\end{enumerate*}
The transformer-based model \cite{vaswani2017attention} relies on a self-attention mechanism to calculate representations for tokens.
For each layer in the transformer network, the representation of each token is updated by computing a weighted sum over all tokens represented in the entire sequence, together with a non-linear feedforward transformation.
The weight value of each token is determined using self-attention, which computes a similarity score for each pair of tokens in the sequence.
For every token in every layer of the network, an attention layer is used.
As a result, it is non-trivial to aggregate all the importance values per layer into one interpretable importance score per token.

In our work, we follow \citet{abnar-zuidema-2020-quantifying} by using the attention rollout scores for [CLS] tokens.
For the transformer-based model, we fine-tune the lighter, pre-trained DistilBERT variant \cite{sanh_distilbert_2019} instead of the full BERT model \cite{devlin-etal-2019-bert} to reduce the computational overhead and ecological footprint. For classification, we add a linear layer on top of the pooled output. We concatenate pair-sequences with a [SEP] token.

Unlike a transformer-based model, an LSTM-based model processes input tokens in sequential order.
For each token, $t \in \{1, \dots, n\}$, in the sequence, the global representation of the sequence, $\textbf{h}_{t}$, is updated using a non-linear transformation with an embedded representation of
\begin{inparaenum}[(i)]
	\item the current token in the sequence, $t_{n}$, and
	\item the previous global representation of the sequence, $\textbf{h}_{t-1}$, as input.
\end{inparaenum}
A bidirectional LSTM (Bi-LSTM) uses two stacked LSTMs to process tokens in both directions of the sequence.

For the Bi-LSTM, we use the same single-layered bidirectional encoder with the query-less additive $(tanh)$ attention and linear feedforward decoder as in the work of \citet{jain-wallace-2019-attention}.
As a result, the attention weight for each token in the sequence is solely based on its representation $\textbf{h}_{t}$, and not on a query item.
In pair-sequence tasks, we embed, encode, and induce attention over each sequence separately.
The decoder predicts the label from the concatenation of
\begin{inparaenum}[(i)]
	\item both context vectors $c_{1}$ and $c_{2}$,
	\item their absolute difference $|c_{1} - c_{2}|$, and
	\item their element-wise product $c_{1} \odot c_{2}$.
\end{inparaenum}

\subsection{Training the models}

We train three independently-seeded instances of each of the models described in Section~\ref{section:models} using the AllenNLP framework \cite{gardner-etal-2018-allennlp}, each for a maximum of 40 epochs. We use a patience value of 5 epochs for early stopping.
For DistilBERT, we fine-tune the standard ``base-uncased'' weights available in the HuggingFace library \cite{Wolf2019HuggingFacesTS} with the AdamW \cite{loshchilov2019decoupled} optimizer.
For the BiLSTM, we follow \citet{jain-wallace-2019-attention} and select a 128-dimensional encoder hidden state with a 300-dimensional embedding layer. We tune pre-trained FastText embeddings \cite{bojanowski-etal-2017-enriching} and optimize with the AMSGrad variant \cite{Tran_2019} of Adam \cite{kingma_dp_adam_2015}. Appendix \ref{app:performance} details model performance and indicates both the BiLSTM and DistillBERT are sufficiently accurate for our analysis.

\subsection{Explaining the models}

We leverage the Captum\footnote{\url{https://github.com/pytorch/captum}} implementations of LIME, Integrated Gradients, DeepLIFT, Grad-SHAP, and Deep-SHAP, and use the padding token as a baseline where applicable. For LIME, we mask tokens as features and use 1000 samples to train the interpretable models.
We apply our feature attribution methods to the predictions of each independently-seeded model for 500 instances randomly sampled from each test set.

Our code is publicly available.\footnote{\url{https://github.com/sfschouten/court-of-xai}} Refer to Appendix \ref{app:reproducibility} for more information on reproducing our experiments.

\section{Results}\label{sec:results}

In this section, we answer our main research question: \emph{How well do attention-based explanations correlate with recent feature attribution methods for NLP tasks?}
In general, we find that attention-based explanations do not correlate strongly with feature attribution methods, with some exceptions (see Section~\ref{sec:results2} and~\ref{sec:results3}).

\subsection{RQ1: Does the correlation depend on the model architecture?}
\label{sec:results1}
Tables~\ref{tab:distilbert_results} and~\ref{tab:bilstm_results} display the average\footnote{Across the 3 model instances, randomly selecting 500 instances from the test set using the training seed.}  Kendall-$\tau$ correlations between the explanation methods for the DistilBERT and BiLSTM models, respectively.\footnote{See Appendix \ref{app:results-with-std-dev} for the same tables including the standard deviation.}
A stronger correlation (i.e., agreement) is indicated by a darker blue colour in the table cell.
In general, we see that the agreement between explanation methods is substantially lower for the DistilBERT model than for the BiLSTM model.

\input{tables/results/DistlBERT}

For the DistilBERT model, we observe a weak correlation across all explanations -- almost none of the explanation methods we test, whether they are attention-based or feature attributions, correlate strongly with one another.
There are two exceptions: (i) Integrated Gradients moderately correlates with Grad-SHAP, and (ii) DeepLIFT moderately correlates with Deep-SHAP.
However, this is unsurprising since the implementation of Grad-SHAP is based on Integrated Gradients, and the implementation of Deep-SHAP is based on DeepLIFT.
In contrast, it is surprising to see the lack of correlation between Grad-SHAP and Deep-SHAP for the DistilBERT model, given that they are different implementations of the same algorithm, SHAP \cite{shap}.

\input{tables/results/BiLSTM}

For the BiLSTM model, we observe a weak correlation between attention-based explanations and feature attribution explanations (i.e., the first row of Table~\ref{tab:bilstm_results}).
Similar to the DistilBERT model, we see strong correlations for the methods that have similar underlying implementations.
We also see some strong correlation between feature attribution methods, especially for single-sequence tasks (see Section~\ref{sec:results2}).

Overall, we conclude that the correlation between explanation methods depends on the model architecture, which answers \textbf{RQ1}.
In general, the correlation is weaker for the DistilBERT model than for the BiLSTM.
However, the overall agreement between attention-based explanations and feature attribution explanations is weak, regardless of the model architecture.

\subsection{RQ2: Does the correlation depend on the nature of the classification task?}
\label{sec:results2}
To investigate \textbf{RQ2}, we examine five different datasets corresponding to two different tasks: IMDb and SST-2 are single-sequence tasks, while MNLI, Quora, and SNLI are pair-sequence tasks and show the results in Tables~\ref{tab:distilbert_results} and~\ref{tab:bilstm_results}.

For the BiLSTM model, the results are clear: there is a substantially stronger correlation among explanations for single-sequence tasks than for pair-sequence tasks.
For the DistilBERT model however, we do not see such a clear distinction -- there is, at best, a weak correlation across all explanation methods, regardless of the nature of the task.

Regarding \textbf{RQ2}, we conclude that explanation methods agree with each other more on single-sequence tasks than pair-sequence tasks. This relationship holds especially for the BiLSTM model compared to the DistilBERT model.

\subsection{Explanations that strongly correlate}
\label{sec:results3}

Given that some of the feature attribution methods are inherently related, we would expect to see stronger degrees of correlation between them:
\begin{enumerate*}[label=(\roman*)]
	\item  Grad-SHAP and Deep-SHAP are two different versions of the SHAP explanation method \cite{shap},
	\item  Grad-SHAP relies on the Integrated Gradients implementation to compute the token importances, and
	\item Deep-SHAP relies on DeepLIFT.
\end{enumerate*}
We see some strong correlations for (ii) and (iii) between these methods for both models.
These correlations are especially strong for the single-sequence tasks on the BiLSTM.
However, we see remarkably weak correlations for (i) on DistilBERT, and for the BiLSTM on the pair-sequence tasks.

\section{Discussion} \label{sec:discussion}
Based on the results in Section~\ref{sec:results}, we argue that rank correlation with existing feature attribution methods is not an appropriate measure for evaluating attention-based explanations.
In this section, we detail three main reasons for this:
\begin{enumerate*}[label=(\roman*)]
	\item the general lack of correlation across all explanation methods, especially for transformer-based models,
	\item the fact that similar explanations do not always result in correlated rankings, and
	\item the lack of justification for the existence of one ``ideal'' explanation, which is a fundamental assumption of the \emph{agreement as evaluation} paradigm.
\end{enumerate*}

\subsection{Lack of correlation between explanation methods.}
In Section~\ref{sec:results}, we have shown that there is a low degree of correlation between explanation methods, especially for the transformer-based model.\footnote{Although there are some examples of stronger correlations, these are either on a very specific combination of task and model (i.e., single sequence tasks on the BiLSTM), or due to similarities in the underlying implementation of the explanation methods. } Similar conclusions are observed in the work of \citep{feldhus-etal-2021-thermostat} and \citep{krishna_disagreement_2022}.
This makes it challenging to justify the expectation that in order for attention-based explanations to be valid, they should correlate with existing feature attribution methods.
If none of the existing feature attribution methods correlate with one another (as is the case for the transformer-based model), should we expect attention-based explanations to correlate with them?
Therefore, we argue against the use of rank correlation as an ``off-the-shelf'' tool to evaluate attention-based explanations.

A common critique of using Kendall-$\tau$ is that expecting the full rankings of token importance to correlate is unrealistic. In their original paper, \citet{jain-wallace-2019-attention} theorized that top-$k$ token comparison would reduce noise, improve correlation, and more closely align with the end user's interest in the most salient features. However, formulating a systematic method for calculating $k$ across various models, tasks, and datasets is difficult in practice and yields mixed results. \citet{treviso-martins-2020-explanation} show a dynamic selection of $k$ with a sparse attention mechanism more effectively conveys justifications of decisions than fixed values of $k$ (e.g., 5 or 10 tokens). In contrast, \citet{krishna_disagreement_2022} show a widespread disagreement of feature attribution methods across domains regardless of whether $k$ is fixed or dynamic (e.g., the top $X$\% of features or top $X$\% of tokens based on average sentence length).

\subsection{Similar explanations do not imply strongly correlated rankings}
In order to investigate what weak correlation looks like in practice, we apply a heatmap of importance scores for each of the tested explanation methods.
The darker the color, the more important the token is w.r.t. the model prediction.
Figure~\ref{fig:distill-ss} shows a sentiment analysis example using the transformer-based model on the SST-2 dataset.
Upon first glance, the visualizations in Figure~\ref{fig:distill-ss} look relatively similar: the methods are all highlighting roughly the same words -- almost all methods (except for Int-Grad) highlight the word `technically' as important.
However, the average Kendall-$\tau$ correlation across all methods for this example is very weak -- only $0.01$. This highlights the problem of evaluating explanations by measuring the correlation between two ranked lists of importance scores: two explanation methods can roughly indicate the same token(s) as important, but when similar tokens do not end up in similar positions in the two rankings, the overall agreement between the two methods can still be low.

\begin{figure}
	\centering
	\begin{subfigure}{.48\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{plots/figures/distill-ss-sample1074.pdf}
		\caption{Sentiment analysis example from the SST-2 dataset. The average Kendall-$\tau$ correlation across all methods for this example is $0.01$. }
		\label{fig:distill-ss}
	\end{subfigure}\hspace{1mm}
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{plots/figures/distill-ps-sample481.pdf}
		\caption{Textual entailment example from the MNLI dataset.The average Kendall-$\tau$ correlation across all methods for this example is $0.05$. }
		\label{fig:distill-ps}
	\end{subfigure}
	\caption{Examples of explanations for the transformer-based model. The brighter the color, the higher the attribution value}
	\label{fig:test}
\end{figure}

\subsection{Is there one ideal explanation?}
The \emph{agreement as evaluation} paradigm implicitly assumes the existence of a single ``ideal'' explanation (i.e. ranking of tokens) which all methods must uncover, and new feature attribution methods are evaluated based on how strongly they correlate with the ``ideal'' explanation. However, it is unclear whether this assumption holds.
For instance, input token importance rankings may only capture a narrow slice of the model's behavior, such that many plausible rankings exist.

Since many tasks may be too complex for humans to judge token-level importance, it can be unclear how to choose the ``ideal'' ranking.
While a handful of highly polar tokens are generally indicative of the label in binary sentiment classification \cite{sun-lu-2020-understanding}, annotators may be unsure how to rank the other tokens. The difficulty increases in the pair-sequence setting -- if two words indicate a contradiction, which is more important?
As a result, when agreement is measured in the presence of multiple faithful and plausible rankings, feature attribution methods may look deceptively problematic, even if they are not.

Figure~\ref{fig:distill-ps} shows an example of how different explanation methods can highlight different tokens as being important for the prediction, even though the underlying model is the same.
Figure~\ref{fig:distill-ps} shows a textual entailment example from the MNLI dataset, where the task is to predict whether or not the second sentence entails the first.
The first sentence comes after the [CLS] token: \textit{``But there are many more who still need our help.''}.
The second sentence comes after the first [SEP] token: \textit{``10,000 people still need our help''}.
Given that we see different explanations for the same example using the same model, it is unclear how to identify which of these explanations is the ``ideal'' one.

Instead of evaluating potential explanations based on how well they correlate with one another, we suggest that researchers and practitioners should test various explanation methods and employ a human-in-the-loop process to determine if the explanations align with human intuition for the particular use case at hand.
This can be used in combination with metrics such as those suggested by \citet{atanasova-etal-2020-diagnostic}.

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

In this work, we investigate how strongly attention-based explanations correlate with recent feature attribution methods for NLP tasks.
We compare attention-based explanations to five recent feature attribution methods, using both transformer- and LSTM-based models, on both single- and pair-sequence tasks.
Overall, we observe a low degree of correlation between the attention-based explanations and the feature attribution explanations.
For the transformer-based model, we find a weak correlation across all explanations for both task types.
For the BiLSTM, we observe some strong correlations between explanations, but only for simple single-sequence tasks.

Through our experiments, we discover
\begin{enumerate*}[label=(\roman*)]
	\item a general lack of correlation between explanation methods, especially for more complex settings (i.e., transformer-based model and pair-sequence tasks) which is corroborated by additional recent research on text, tabular, and image data \cite{feldhus-etal-2021-thermostat,krishna_disagreement_2022},
	\item that similar explanations do not always result in correlated rankings, and
	\item the existence of a single ``ideal'' explanation is questionable, which is a fundamental assumption of the \emph{agreement as evaluation} paradigm.
\end{enumerate*}
Without an external ground-truth explanation, all that rank correlation tells us is whether or not two rankings are similar.
For this reason, we recommend practitioners stop using \textit{agreement as evaluation} for attention-based explanations.

In future work, we plan to approach this problem from first principles by formulating toy data to guarantee a single, correct top-$k$ or full ranking for each instance to see if XAI methods are consistently capable of recovering the ideal explanation under this setting.

\section*{Acknowledgements}
We thank the anonymous reviewers, Bilal Alsallakh and Maarten de Rijke for their helpful discussions and suggestions.
This research was supported by Ahold Delhaize and the Netherlands Organisation for Scientific Research under project nr. 652.001.003 and the Nationale Politie.
All content represents the opinion of the authors, which is not necessarily shared or endorsed by their respective employers and/or sponsors.