\section*{Introduction}
Rank one perturbations of random Hermitian matrices appear naturally in a variety of contexts. An important example is the celebrated phenomenon of \textit{BBP transition} (for Baik-Ben Arous-Péché, see \cites{BBP_paper, Peche_BBP_GUE}), which arises when the perturbation is a positive rank-one Hermitian matrix, and so slightly `pushes' the spectrum to the right, to the point where one outlier is clearly separated from the bulk. Another example is that of anti-Hermitian perturbations, which play an important role in scattering theory (for a general presentation of this application, see Chapter 34 of \cite{RMT_handbook} and \cite{FyodorovSommers1997}). The questions we consider here have much in common with these two examples, and could be summarized as follows: what can be rigorously established about rank-one anti-Hermitian perturbations of a random Hermitian matrix considered dynamically, that is, when the coupling parameter is interpreted as time? And more specifically, what can be said about the eventual outlier? \medskip

Throughout the paper, the essential assumption is that $H$ be a random Hermitian matrix for which the uniform isotropic local law (Theorem \ref{thm:UILL}) is known to hold. For the sake of definiteness, say that we consider the Wigner ensemble with the following standard assumptions: entries $h_{ij}, i \leq j$ are independent, off-diagonal (resp. diagonal) entries are identically distributed with continuous distribution on $\mathbb{C}$ (resp. $\mathbb{R}$) such that $\mathbb{E} h_{ij}=0$, $\mathbb{E} |h_{ij}|^2=1/N$ and finite moments, i.e. $\mathbb{E} |\sqrt{N} h_{ij}|^p\le C_p$ for all $p$. These are the Wigner matrices with which we work by default -- although the method and results also hold under more general conditions, such as those of \cite{ErdoesKruegerSchroeder}. As the entries of $H$ are assumed to have a continuous distribution, it holds almost surely that $H$ has $N$ distinct real eigenvalues $\mu_1 < \dots < \mu_N$; and we denote by $(u_i)_{i=1}^N$ a choice of associated  normalized  eigenvectors. \medskip

We consider the following process, which is a rank one perturbation of $H$:
\begin{equation}\label{def_system}
G = G_t := H + i t vv^*, \qquad t \in \mathbb{R},
\end{equation}
where $v$ is a random unit vector, chosen uniformly on the sphere and independent of $H$. However, the randomness of $v$ is not a very relevant feature (as long as independence holds). The main results, indeed, are proved for any fixed $v$; only the proofs of some preliminary facts are greatly simplified when stating them with respect to the randomness in $v$. \medskip

It is straightforward to check that the eigenvalues of $G_t$ lie in the upper half-plane for $t>0$, and that $G_{-t}=G_t^*$, so that the eigenvalue trajectories for $t<0$ and $t>0$ are symmetric to the real axis.
Another deterministic property is that, as $t \rightarrow \infty$, the spectrum is composed of one outlier that diverges ($\lambda_{j_{\text{out}}}(t) \approx it$) and $N-1$ eigenvalues that converge to specific locations on the real line.
Note that the distribution of eigenvalues at any fixed $t$ is known when $H$ is taken from an integrable ensemble such as GUE or GOE (see for instance \cites{FyodorovKhoruzhenko,FyodorovSommers}). However, the questions we ask here are of a dynamic nature: we are interested in the evolution of the spectrum $\{ \lambda_1(t), \dots, \lambda_N(t) \}$ of $G_t$ when $t \in [0, +\infty) $. The fact that indices can be given consistently to form $N$ continuous trajectories is a consequence of the non-intersection of trajectories (Theorem \ref{thm:non_crossing}); we choose these indices so that $\lambda_i(0)=\mu_i$.
\medskip

Section \ref{first_prop_sec} introduces some basic properties of these dynamics, either deterministic or probabilistic. In particular, it is established that trajectories are almost surely non-crossing (Theorem \ref{thm:non_crossing}); moreover, they are everywhere differentiable and follow an ordinary differential equation (Theorem \ref{thm:ode}), which brings to light a link between these trajectories and the eigenvector overlaps of the same system, studied in \cite{FyodorovMehlig}. Seen from this angle, the system appears to be extremely unstable, and therefore unpredictable. \medskip

Section \ref{subtle_prop_sec} relies on the isotropic local law (Theorem \ref{thm:UILL}), borrowed from the existing literature on Hermitian random matrices, to give more precise high-probability estimates on these trajectories. Most importantly, we establish in Theorem \ref{thm:the_outlier} that the outlier is distinctly separated from the rest of the spectrum at all times $t > 1+N^{-1/3+\varepsilon}$, with $\varepsilon>0$. This is a first step towards proving that the relevant timescale is exactly $1+ \mathrm{O}( N^{-1/3})$, which happens to be the same timescale as for the BBP transition. Inspired by the current manuscript, Yan Fyodorov \cite{FyodorovPrivateCom} provided further heuristics for this timescale, based upon an explicit formula for the density of the eigenvalues in the case when $H$ is GUE \cite{FyodorovSommers1996}. \medskip

A natural question that is left open is that of the origin of the outlier: from which eigenvalue $\mu_j$ of $H$ is this particular trajectory more likely to originate? In the context of a Hermitian perturbation, the answer is  trivial;
 for an anti-Hermitian perturbation it becomes very subtle. Heuristic arguments as well as numerical simulations seem to imply that the eigenvalues closer to the origin are much more likely to become the outlier when $t$ increases. However, the absence of a local law very near the spectrum
prevents us from turning this phenomenology into a rigorous statement.

\subsection*{Notations and conventions}

We introduce the following standard definition.

\begin{definition}[High Probability]
A sequence of events $(A_N)_{N \geq  1}$ is said to happen with high probability if for any $D>0$ the inequality
$$
\mathbb{P} \left( A_N^c \right) < N^{-D}
$$
holds for sufficiently large $N$.
\end{definition}

It is customary, when working with Wigner  matrices, to define the function $m_{\text{sc}}$ on $\mathbb{C} \backslash [-2,2]$,
the Stieltjes transform of the Wigner semicircle distribution $\rho_{\text{sc}}(x) =\frac{1}{2\pi}\sqrt{4-x^2}$ on $[-2,2]$.
This function is the natural approximation of the resolvent on both the upper and lower half-planes.
Note that $m_{sc}$ has a jump discontinuity on $[-2,2]$
and $m_{\text{sc}}(\overline{z}) = \overline{m_{\text{sc}}(z)}$.
In this paper we  need  a slight modification of this function on the lower half-plane that is holomorphic through $[-2,2]$,
i.e. we define
\begin{equation}\label{Stieltjes_msc}
\mathfrak{m} (z) =\frac{-z + \sqrt{z^2-4}}{2}
\end{equation}
which is holomorphic on $\mathbb{C} \backslash (-\infty,-2]\cup[2,+\infty)$ with the appropriate choice of branch-cut for the square root, such that $\Im ( \sqrt{z^2-4} ) >0$ for every $z$ in this domain. In particular, $\mathfrak{m}=m_{\text{sc}}$ on the upper-half plane. It is
a solution to the equation $\mathfrak{m}(z)^2 + z \mathfrak{m}(z)+1 =0$, so that $\mathfrak{m}$ defines a bijection from $\mathbb{C} \backslash (-\infty,-2]\cup[2,+\infty)$ to its image with
\begin{equation}\label{inverse_m}
z= -\frac{1+\mathfrak{m}(z)^2}{\mathfrak{m}(z)}.
\end{equation}
We also define, for any $t>0$,
\begin{equation}
t^* := t - \frac{1}{t},
\quad \text{such that} \quad
\mathfrak{m}(it^*) = \frac{i}{t}.
\end{equation}
It is a very important fact (for most results in Section \ref{subtle_prop_sec}) that the holomorphic function $\mathfrak{m}(z)-i/t$ has only one zero at $z=it^*$, with multiplicity one; this zero being in the upper-half plane if and only if $t\geq 1$. The role of the quantity $t^*$ as an approximation to the resonant eigenvalue  was already noted, see for instance \cite[p.1950]{FyodorovSommers1997}.

\newpage
\section{First properties of trajectories}\label{first_prop_sec}

In this section, we consider a fixed matrix $H$. The results are either deterministic, or stated with respect to $\mathbb{P}_v$, the randomness in $v$, a uniform unit vector.

\subsection{Weighted resolvent and non-intersection of trajectories}

For now, one can assume that the indices of the eigenvalues $\lambda_j(t)$ of $G_t$ are given arbitrarily for each $t$. One of the goals of the following results is to establish that the indices can be given in a consistant way, with $\lambda_j(t)$ being `the' trajectory such that $\lambda_j(0)=\mu_j$.

\begin{definition}
We define the weighted resolvent associated to $H$ and the unit vector $v$ by
\begin{equation}\label{weighted_resolvent}
    \mathscr{W} (z) := \sum_{j=1}^N \frac{|\langle u_j | v \rangle|^2}{\mu_j - z} = \langle v | (H-z)^{-1} v \rangle.
\end{equation}
\end{definition}
The name \textit{weighted resolvent} refers to the fact that  $\mathscr{W}(z)$ can be considered as
a weighted sum with weights $|\langle u_j | v \rangle|^2$ that sum up to one. Since $\mathbb{E}_v |\langle u_j | v \rangle|^2 = \frac{1}{N}$,
where expectation is with respect to the uniform unit vector $v$,
so $ \mathbb{E}_v \mathscr{W}(z) $ is the usual trace of the resolvent of $H$.

\begin{proposition}\label{level_lines}
For any $t \neq 0$,
\begin{equation}
z \in \mathrm{Sp} \left( G_t \right)
\quad
\Leftrightarrow \quad \mathscr{W}(z)= \frac{ i }{ t }.
\end{equation}
As a consequence, the trajectories of eigenvalues for the system (\ref{def_system}) are given by the zero level lines of the real part of the weighted resolvent $\mathscr{W}$:
\begin{equation}
\bigcup_{j=1}^N
\{
\lambda_j(t) \ : \ \ t \in \mathbb{R}^*
\}
=
\{
z \in \mathbb{C} \backslash \mathrm{Sp} H \ : \ \Re \mathscr{W}(z) = 0
\}.
\end{equation}
In particular,
\begin{equation}\label{real_part_bound}
\forall j, t, \quad \mu_1 \leq \Re \lambda_j(t) \leq \mu_N
\end{equation}
with equality happening only for $t=0$, or if $|\langle u_1 | v \rangle|^2=1$ (resp. $|u_N^* v|^2=1$).
\end{proposition}

\begin{proof}
For any $z \notin \mathbb{R}$, the matrix $(H-z)^{-1} vv^*$ has rank $1$, and in particular its trace $\mathscr{W}(z)$ is its only non-zero eigenvalue. We write:
\begin{equation}
\det (G_t -z )
=
\det (H-z) \det\left(I_N + it (H-z)^{-1} vv^* \right)
=
\det (H-z) \left(1 + it \mathscr{W}(z) \right).
\end{equation}
The result follows, as $ z \in \mathrm{Sp} (G_t)$ is  equivalent to $1 + it \mathscr{W}(z)=0$. This is equivalent to $z = E+i \eta$ being a solution to the equation
\begin{equation}\label{algebraic_variety}
    \sum_{j=1}^N |\langle u_j | v \rangle|^2 \frac{\mu_j - E}{(\mu_j -E)^2 + \eta^2} = 0.
\end{equation}
and \eqref{real_part_bound} follows by inspection, all terms having the same sign outside the vertical strip $\{\mu_1 \leq \Re z \leq \mu_N\}$.
\end{proof}

We now prove that the trajectories almost surely do not cross. This is in fact true for any $N$; we give below a concise argument that requires $N \geq 5$, which is enough for our purpose, as subsequent results concern large values of $N$.

\begin{theorem}\label{thm:non_crossing} We assume $N \geq 5$, and that $H$ has distinct real eigenvalues $\mu_1 < \dots < \mu_N$. Almost surely with respect to $\mathbb{P}_v$, the randomness in $v$, the trajectories of the system \eqref{def_system} do not intersect, nor do they self-intersect: that is, for all $i \neq j$,
\begin{equation}
\mathbb{P}_v \left(
\exists t,s \geq 0, \ \lambda_i (t) = \lambda_j (s)
\right)
= 0
\quad \& \quad
\mathbb{P}_v \left(
\exists t,s \geq 0, t \neq s, \ \lambda_j (t) = \lambda_j (s)
\right)
= 0 .
\end{equation}
\end{theorem}

\begin{proof}
If for some $j_1, j_2$, $\lambda_{j_1} (t) = \lambda_{j_2} (s) = z$, then by Proposition \ref{level_lines},
\begin{equation}\label{equality}
\mathscr{W}(z) = \frac{ i }{ t } = \frac{ i }{ s }
\end{equation}
and so $t=s$ (if $t$ or $s=0$, \eqref{equality} is understood as meaning that $z$ is a pole of $\mathscr{W}$). This rules out self-intersection, as well as intersection of two distinct trajectories for different times. Intersection of two trajectories at $t=0$ is ruled out by the eigenvalues of $H$ being distinct. The remaining possibility is that two trajectories intersect at the same time $t>0$, which implies that $\mathscr{W} - i/t$ vanishes at $z \notin \mathbb{R}$ with multiplicity at least $2$; in particular $\mathscr{W}'(z)=0$. This is to say that an intersection point $z$ is such that the two conditions
\begin{equation}\label{two_conditions}
\Re \mathscr{W} (z) = 0
\quad \& \quad
\mathscr{W} ' (z) = 0
\end{equation}
are met. We will prove that almost surely, there is no $z$ that checks both conditions. First, note that
\begin{equation}
\mathscr{W}'(z) = \sum_{j=1}^{N} \frac{ |\langle u_j | v \rangle |^2 }{ (\mu_j - z)^2 }
\end{equation}
and so the condition $\mathscr{W}'(z) = 0$ is equivalent to $z$ being the root of a real polynomial of degree $2N-2$ (almost surely) that does not vanish on the real line. There is a finite number of such points: namely, almost surely $N-1$ conjugated pairs that we denote $(Z_k,\overline{Z_k})_{k =1}^{N-1}$, with $\Im Z_k >0$, counted with multiplicity.
For any $z \notin \mathbb{R}$, let us define the real vectors:
\begin{equation}\label{three_real_vectors}
Y_1(z) := \left( \Re \frac{1}{ \mu_j -z } \right)_{j=1}^N, \quad
Y_2 (z) := \left( \Re \frac{1}{ (\mu_j -z)^2 } \right)_{j=1}^N, \quad
Y_3 (z) := \left( \Im \frac{1}{ (\mu_j -z)^2 } \right)_{j=1}^N,
\end{equation}
and notice that the first condition of (\ref{two_conditions}) can be written as
\begin{equation}\label{first_condition}
\langle X \ | \ Y_1(z) \rangle = 0
\end{equation}
where $X := (|\langle u_j | v \rangle |^2)_{j=1}^N$, and the second one similarly as
\begin{equation}\label{second_condition}
\langle X \ | \ Y_2(z) \rangle =0, \qquad  \langle X \ | \ Y_3(z) \rangle =0.
\end{equation}
We will rely on the (deterministic) fact that $Y_1(z)$ is not in the $\mathbb{R}$-span of $Y_2(z), Y_3(z)$ for any $z \in \mathbb{C} \backslash \mathbb{R}$.

\begin{lemma}
Assuming $N \geq 5$, for any $z = E+i \eta \notin \mathbb{R}$,
$
Y_1(z) \notin \mathrm{Span}_{\mathbb{R}} (Y_2(z), Y_3(z)).
$
\end{lemma}

\begin{proof}
For any $z = E+i \eta$ with $\eta \neq 0$, we denote
\begin{equation}
a_j + i b_j = \frac{1}{\mu_j - z}.
\end{equation}
In particular, every pair $(a_j, b_j)$ solves
\begin{equation}\label{circle_eta}
b_j = \eta (a_j^2 + b_j^2)
\end{equation}
which is the equation of the circle
$\mathscr{C}_{\eta}$ with center $\frac{i}{2\eta}$ and radius $\frac{1}{2\eta}$, in the $(a,b)$ plane. Moreover, the $\mu_j$'s being distinct, $(a_j,b_j)$ are $N$ distinct points on $\mathscr{C}_\eta$. \\ \\
Moreover, assuming that $Y_1(z) = \alpha Y_2(z)+\beta Y_3(z)$ for some $\alpha, \beta \in \mathbb{R}$, we have
\begin{equation}\label{hyperbola}
a_j = \alpha (a_j^2 - b_j^2) + 2 \beta a_j b_j ,
\end{equation}
which is the equation of a hyperbola $\mathscr{H}_{\alpha, \beta}$ (including the degenerate case, that yields a union of two lines). By general theory (e.g. Bézout's theorem for curves), we have
\begin{equation}\label{Bezout_curve}
N \leq  | \mathscr{C}_{\eta} \cap \mathscr{H}_{\alpha,\beta} | \leq 4,
\end{equation}
which is a contradiction.
\end{proof}

For any fixed $z \in \mathbb{C} \backslash \mathbb{R}$, conditionally on (\ref{second_condition}), the probability that (\ref{first_condition}) holds as well is zero, as $Y_1(z)$ is linearly independent of $Y_2(z), Y_3(z)$ and $X$ has a continuous distribution on the simplex
$$ \left\{ (x_1, \dots, x_N) \ : \ \ x_j \geq 0 , \ \sum_{j=1}^N x_j =1 \  \right\}.$$
So what the above argument allows to conclude is that
\begin{equation}
\forall z \in \mathbb{C} \backslash \mathbb{R}, \
\mathbb{P}_v \left( \Re \mathscr{W}(z) = 0 \ | \ \mathscr{W}'(z)=0 \right) = 0.
\end{equation}
We want to deduce from this that
\begin{equation}\label{no_z_exists}
\mathbb{P}_v \left( \exists z \in \mathbb{C} \backslash \mathbb{R}, \
\Re \mathscr{W} (z) = 0
\ \& \
\mathscr{W} ' (z) = 0  \right) = 0.
\end{equation}
This follows for the following reasons, which essentially rely on the roots $\{Z_k\}_{k=1}^{N-1}$ of $\mathscr{W}'$ being finitely many. Let $\pi \in \mathfrak{S}_{N-1}$ be a uniform random permutation (in the rest of the proof, $\mathbb{P}$ denotes the probability with respect to $\pi$ as well as the previous randomness in $v$), and $N_z$ denote the order of vanishing of $\mathscr{W}'$ at $z$, that is
$$
N_z := \# \{ k \ | \ Z_k=z \} .
$$
Note the equality of events $\{ \mathscr{W}'(z) = 0\} = \{ N_z \geq 1\}.$ The randomness in $\pi$ ensures
$$
\mathbb{P} \left( Z_{\pi(k)}=z \ | \ N_z =m \right) = \frac{m}{N-1} \geq \frac{1}{N-1},
$$
and so, by partitioning the event $\{ N_z \geq 1 \}$,
$$
\mathbb{P} \left( Z_{\pi(k)}=z \ | \ N_z \geq 1 \right)
\geq
\frac{1}{N-1}.
$$
We then write, for any $k$,
\begin{align*}
0
 = \mathbb{P} \left( \Re \mathscr{W}(z) = 0 \ | \ N_z \geq 1 \right)
& \geq \mathbb{P} \left( Z_{\pi(k)} =z \ | \ N_z \geq 1 \right)
\mathbb{P} \left( \Re \mathscr{W}(z) = 0 \ | \  Z_{\pi(k)}=z \right) \\
& \geq \frac{1}{N-1}
\mathbb{P} \left( \Re \mathscr{W}(Z_{\pi(k)}) = 0 \ | \  Z_{\pi(k)} =z \right),
\end{align*}
and deduce that the right hand side is zero. We then integrate over the upper half-plane and conclude that for every $k$,
$$\mathbb{P} \left( \Re \mathscr{W}(Z_{\pi(k)}) =0 \right) = 0.$$
As there are finitely many $Z_k$'s, we can use the union bound,
\begin{equation}
\mathbb{P}_v \left( \exists k, \ \Re \mathscr{W}(Z_k) =0 \right)
= \mathbb{P} \left( \exists k, \ \Re \mathscr{W}( Z_{\pi(k)} ) =0 \right)
\leq \sum_{k=1}^{N-1} \mathbb{P} \left( \Re \mathscr{W}(Z_{\pi(k)}) =0 \right) = 0.
\end{equation}
This is to say that, almost surely, no root of $\mathscr{W}'$ satisfies \eqref{first_condition}, and the claim \eqref{no_z_exists} follows.
\end{proof}

An important consequence of Theorem \ref{thm:non_crossing} is the possibility of choosing an coherent indexation such that each $\lambda_j(t)$ is a (uniquely defined) distinct continuous trajectory. We now study the deterministic behavior of these trajectories.

\subsection{Deterministic evolution}

The main features of these $N$ almost surely non-crossing continuous trajectories are as follows:
\begin{enumerate}
\item[(i)] At $t=0$, all eigenvalues are real.
\item[(ii)] For $t>0$, all eigenvalues are in the upper-half plane.
\item[(iii)] When $t \rightarrow \infty$, one eigenvalue (`the outlier' $\lambda_{j_{\text{out}}}$) diverges with $\Im \lambda_{j_{\text{out}}} \xrightarrow[t \rightarrow \infty]{} \infty$ and bounded real part, and the rest of the spectrum converges to $N-1$ distinct points on the real line.
\end{enumerate}
This last property can be established by the Schur complement identity, that allows to identify the limit points as the eigenvalues of the projection of the operator $H$ on the space orthogonal to $v$. \medskip

A remarkable fact is that this evolution of eigenvalues can be described by a closed system of differential equations when the perturbation has rank one.

\begin{theorem}[Differential Equation]\label{thm:ode} Let $H$ be a Hermitian matrix with simple eigenvalues $(\mu_j)_{j=1}^N$ and associated unit eigenvectors $(u_j)_{j=1}^N$. The evolution of the eigenvalues $(\lambda_j(t))_{j=1}^N$ of $G_t = H+i t vv^*$ can be described by the following closed system of equations, as long as the eigenvalues are distinct\footnote{This is almost surely the case when $v$ is a random unit vector, by Theorem \ref{thm:non_crossing}.}. The initial condition is $\lambda_j(0) = \mu_j $. For $t=0$, one has
\begin{equation}\label{initial_push}
\lambda_j '(0) = i |v^* u_j |^2 ,
\end{equation}
and for $t>0$, assuming the eigenvalues to be distinct at that time,
\begin{equation}\label{evolution_eq}
\lambda_j '(t)
=  \frac{i \Im \lambda_j(t)}{t} \prod_{k \neq j} \left( 1+ \frac{2i \Im \lambda_k(t)}{\lambda_j(t) - \lambda_k(t)} \right).
\end{equation}
\end{theorem}

Note that the product in the right hand side of \eqref{evolution_eq} is similar to the formula for diagonal overlaps $\mathscr{O}_{jj}$ in \cite{FyodorovMehlig}, thus establishing a particularly straighforward link between eigenvalue variability and diagonal overlaps in that setting.

\begin{proof}
Let $(L_i,R_j)$ be a bi-orthogonal system of left and right eigenvectors associated to $G=G_t$.
 By standard techniques
 (see for instance Appendix A of \cites{BourgadeDubach}), we find that the derivative of the eigenvalue $\lambda_j$ is given by:
\begin{equation}\label{first_derivative}
\lambda'_j(t) = \langle L_j | G'_t | R_j \rangle
\end{equation}
In the present case, $G'_t = i v v^*$, and so
\begin{equation}\label{first_derivative1}
\lambda_j ' (t)
= i \langle L_j | v \rangle \langle v | R_j \rangle .
\end{equation}
For $t=0$, $G=H$ and $R_j=L_j=u_j$, so that equation (\ref{initial_push}) follows. For $t>0$, we first notice that the quantity $\langle L_1 | v \rangle \langle v | R_1 \rangle $ is invariant under a unitary change of basis, so we can compute it for a Schur form of $G$. As this Schur form $T$ is conjugated to $G$ by a unitary change of variable, we have
\begin{equation}
T = \widetilde{H} + it \tilde{v} \tilde{v}^*
\end{equation}
with $\widetilde{H} = U H U^*$, $\tilde{v} = Uv$. We will simply continue to denote these by $H$ and $v$ in order to not overload notations.
As $T$ is upper-triangular, we have
\begin{equation}
\forall a<b, \quad T_{ba} =  H_{ba} + i t v_b \overline{v_a} =0,
\end{equation}
which implies
\begin{equation}
\forall a<b, \quad  T_{ab}
= H_{ab} + it v_a \overline{v_b}
= \overline{H}_{ba} + it v_a \overline{v_b}
= 2it v_a \overline{v_b},
\end{equation}
and on the diagonal, $T_{aa} = \lambda_a = H_{aa} + it |v_a|^2$, which implies
\begin{equation}\label{real_im_part}
\Re \lambda_a = H_{aa}, \qquad \Im \lambda_a = t |v_a|^2.
\end{equation}
This Schur form can be chosen so that any given eigenvalue is the first on the diagonal, and so we work now with $\lambda_1$
(i.e. we prove \eqref{evolution_eq} for $j=1$) without loss of generality. Another consequence of $T$ being triangular is that $R_1=e_1$ and so $\langle v | R_1 \rangle = \overline{v_1}$. We now compute $\langle L_1 | v \rangle$; in the following argument, we denote $L_1 = (1,\ell_2, \dots,\ell_N)$, and $L_1^{(d)}=(1,\ell_2, \dots, \ell_d)$ so that $L_1=L_1^{(N)}$. This vector satisfies the following recursion, which follows from the definition of $L_1$ and $T$ being triangular:
\begin{equation}\label{Schur_rec}
\ell_{k+1} = \frac{1}{\lambda_1 - \lambda_{k+1}} \langle L_1^{(k)} | \tau_{k+1} \rangle = \frac{2 i t \overline{v_{k+1}}}{\lambda_1 - \lambda_{k+1}}  \langle L_1^{(k)} | v^{(k)} \rangle
\end{equation}
where $\tau_{k+1}$ is the column vector of the first $k$ entries of the $(k+1)$th column of $T$. The recursion for $\langle L_1 | v \rangle$ is then initiated by
$$
\langle L_1^{(1)} | v^{(1)} \rangle = v_1,
$$
and continued in the following way:
\begin{align*}
\langle L_1^{(k+1)} | v^{(k+1)} \rangle & = \langle L_1^{(k)} | v^{(k)} \rangle +\ell_{k+1} v_{k+1} \\
& = \langle L_1^{(k)} | v^{(k)} \rangle + \frac{2 i t |v_{k+1}|^2}{\lambda_1 - \lambda_{k+1}} L_1^{(k)}v^{(k)} \\
& = \langle L_1^{(k)} | v^{(k)} \rangle \left( 1 + \frac{2 i t |v_{k+1}|^2}{\lambda_1 - \lambda_{k+1}} \right)
\end{align*}
where we replaced $\ell_{k+1}$ using (\ref{Schur_rec}).
Eq. (\ref{real_im_part}) gives us $2 i t |v_{k+1}|^2 = \lambda_{k+1} - \overline{\lambda_{k+1}}$, so that
\begin{equation}
    \langle L_1^{(k+1)} | v^{(k+1)} \rangle
    = \langle L_1^{(k)} | v^{(k)} \rangle \left( 1 + \frac{\lambda_{k+1} - \overline{\lambda_{k+1}}}{\lambda_1 - \lambda_{k+1}} \right)
    = \langle L_1^{(k)} | v^{(k)} \rangle  \frac{\lambda_{1} - \overline{\lambda_{k+1}}}{\lambda_1 - \lambda_{k+1}},
\end{equation}
and finally
$$
\langle L_1 | v \rangle = \langle L_1^{(N)} | v \rangle = v_1 \prod_{k=2}^N \frac{\lambda_{1} - \overline{\lambda_{k}}}{\lambda_1 - \lambda_{k}}.
$$
It now follows from~\eqref{first_derivative1} that
$$
\lambda_1 ' (t) = i \langle L_1 | v \rangle \langle v | R_1 \rangle = i |v_1|^2  \prod_{k=2}^N \frac{\lambda_{1} - \overline{\lambda_{k}}}{\lambda_1 - \lambda_{k}},
$$
and eq. (\ref{real_im_part}) allows us to obtain the expression (\ref{evolution_eq}), that is a function of eigenvalues only, valid for any $t>0$.
\end{proof}

\newpage
\section{Properties of trajectories via isotropic local law}\label{subtle_prop_sec}

In this section we state a few estimates on the trajectories of the system \eqref{def_system}, based on the approximation of the weighted resolvent $\mathscr{W}(z)$. The unit vector $v$ is assumed to be fixed in this entire section, and all results are stated with respect to $\mathbb{P}_H$, the randomness of $H$, a Wigner matrix.

We let $T \geq 2$ be an arbitrary fixed constant; the small time ($t \leq T$) and large time ($t>T$) behavior will be analysed somewhat differently.

For any $\zeta, L>0$, we consider the following spectral domains:
\begin{equation}\label{S_strip}
\mathscr{S}_{\zeta} := \{ z=E+i\eta \in \mathbb{C} \ : \ |E| < 3, \  N^{-1+\zeta} \leq \eta < N^{100} \},
\end{equation}
\begin{equation}\label{S_L_strip}
\mathscr{S}_{\zeta,L} := \{ z=E+i\eta \in \mathbb{C} \ : \ |E| < 3, \  N^{-1+\zeta} \leq \eta < L \},
\end{equation}
and
\begin{equation}\label{R_rectangle}
\mathscr{R}_{ \zeta} :=  \{ z=E+i\eta \in \mathbb{C} \ : \ |E| < 3, \ 0 \leq \eta < N^{-1+\zeta}\}.
\end{equation}

The essential input in this section is the uniform isotropic local law, taken over from the existing literature. Local laws in general aim at approximating the resolvent by some deterministic quantity ($m_{\text{sc}}$ times the identity matrix for Wigner matrices); `isotropic' refers  to scalar products $\langle v, (H-z)^{-1}v\rangle$ with some fixed vector $v$, i.e. to the weighted resolvent $\mathscr{W}$, and `uniform' refers to uniformity in the parameter $z$. An isotropic local law, for a given $z$ in a bounded domain, was first given in \cite{KnowlesYin2013}.  We need its following version:

\begin{theorem}[Uniform Isotropic Local Law]\label{thm:UILL}
For any $\zeta, \varepsilon, D>0$ and fixed unit vector $v$,
\begin{equation}\label{UILL_bound}
\mathbb{P}_H \left(
\exists z \in \mathscr{S}_\zeta, \
|\mathscr{W}(z) - \mathfrak{m}(z)| > \frac{N^{\varepsilon}}{\sqrt{N\eta} (1+\eta^2)^{3/4}}
\right)
< N^{-D}.
\end{equation}
\end{theorem}

\begin{proof}
In the bulk spectrum, $|\Re z|\le 2-\epsilon$, this result was stated in Thm 2.1,  eq. (2.6a) and (2.7a) of \cite{ErdoesKruegerSchroeder} even for much more general Hermitian random matrices with possibly correlated entries. The edge regime was settled in Eq. (2.6a) of \cite{AltErdoesKruegerSchroeder}, where the optimal bound is in fact slightly better
than~\eqref{UILL_bound}.
Together, these references provide an isotropic local law for any fixed $z \in \mathscr{S}_{\zeta}$. Uniformity in $z$ can be achieved by
Lipschitz continuity of the functions at stake and using a dense grid of fixed spectral parameters.
\end{proof}

We will use this isotropic local law together with the following classical theorem.

\paragraph{Rouché's Theorem:}
\textit{
Let $f$  and $g$ be two holomorphic functions on a domain $\Omega\subset\mathbb{C}$ with closed and simple boundary $\partial \Omega$. If $|f(z)-g(z)|<|g(z)|$ on $\partial \Omega$, then $f$ and $g$ have the same number of zeros in $\Omega$, counted with multiplicity.
} \medskip

The following proofs all have in common that we determine a domain on which the inequality $|f-g| < |g|$ holds, with $f= \mathscr{W} - i/t$ and $g=\mathfrak{m}-i/t$. The conclusion then either follows immediately, as this strict inequality clearly prevents $f=0$ on the relevant domain, or by applying Rouché's theorem on a Jordan curve: the (deterministic) zeros of $g$ being known, this allows us to draw some conclusions as to the (random) zeros of $f$.

\begin{theorem}\label{thm:elliptic_region}
For any $\varepsilon, \zeta >0$, it holds with $\mathbb{P}_H$-high probability that all eigenvalues $\lambda_i(t)$, for any time $t \in (0,T)$, lie in the domain $\mathscr{E}_{t, \varepsilon} \cup \mathscr{R}_{\zeta}$, where
\begin{equation}\label{elliptic_curve}
\mathscr{E}_{t, \varepsilon} := \{ z=E+i\eta \in \mathbb{C} \ : \
E^2 + \left(\eta - t^* \right)^2 < \frac{N^{\varepsilon}}{N \eta}, \ |E| < 3, 0 \leq \eta < T \}.
\end{equation}
\end{theorem}

\begin{corollary}\label{thm:hyperbolic_region}
For any $\varepsilon >0$, with $\mathbb{P}_H$-high probability, all trajectories up to time $T$ are in the domain
\begin{equation}
\mathscr{H}_{\varepsilon} :=
\left\{ z = E+i\eta \in \mathbb{C} \ : \ \eta E^2 < N^{-1+\varepsilon}, \ |E| <3, \ 0 \leq \eta < T \right\}
.
\end{equation}
\end{corollary}

\begin{proof}[Proof of Theorem \ref{thm:elliptic_region}]
It was proved in~\eqref{real_part_bound} that $\mu_1 \leq \Re \lambda_j(t) \leq \mu_N$ for all $j$ and $t>0$, and it holds with $\mathbb{P}_H$-high probability that $-3< \mu_1 < \mu_N < 3$ (any number larger than $2$ would do), so that the bound $| \Re \lambda_j(t)| <3$ follows immediately. Similarly, with probability one, $0< \Im \lambda_j (t) < t < T$ for all $j$. \medskip

We now justify the main inequality, which results from the following observations. \\
First, by the uniform isotropic law (Theorem \ref{thm:UILL}) for any $\varepsilon, \zeta>0$, with $\mathbb{P}_H$-high probability
\begin{equation}\label{up}
\forall z \in \mathscr{S}_\zeta, \quad
 |\mathscr{W}(z) - \mathfrak{m}(z)| < \frac{N^{\varepsilon/4}}{\sqrt{N \eta}}.
\end{equation}
Second, using \eqref{inverse_m}, we write:
\begin{equation}\label{lipstick}
|z_1 - z_2|
 = \left|\frac{1+\mathfrak{m}(z_1)^2}{\mathfrak{m}(z_1)} - \frac{1+\mathfrak{m}(z_2)^2}{\mathfrak{m}(z_2)} \right|
\leq |\mathfrak{m}(z_1) - \mathfrak{m}(z_2)| \left(1 + \frac{1}{|\mathfrak{m}(z_1) \mathfrak{m}(z_2)|} \right) .
\end{equation}
We apply this to $z_2 = it^*$, $z_1=E+i\eta $ with $|E|<3, t<T, \eta < 2T$ and conclude that in the bounded domain $\mathscr{S}_{\zeta, 2T}$,
\begin{equation}\label{mi}
|z-it^*| < C_T \left| \mathfrak{m}(z) - \frac{i}{t} \right|
\end{equation}
with a constant $C_T$ that only depends on $T$.
Finally, for any $z \in \mathscr{S}_{\zeta,2T} \backslash \mathscr{E}_{t,\varepsilon}$, we have (by definition of $ \mathscr{E}_{t,\varepsilon}$)
\begin{equation}\label{last}
|z-it^*| > \frac{N^{\varepsilon/2}}{\sqrt{N \eta}},
\end{equation}
and therefore the following sequence of inequalities holds for $z \in \mathscr{S}_{\zeta,2T} \backslash \mathscr{E}_{t,\varepsilon}$, bringing together \eqref{up}, \eqref{last}, and \eqref{mi}.
\begin{equation}\label{put_together}
|\mathscr{W}(z) - \mathfrak{m}(z)|
< \frac{N^{\varepsilon/4}}{\sqrt{N \eta}}
< N^{-\varepsilon/4} |z-it^*|
< \left| \mathfrak{m}(z) - \frac{i}{t} \right|,
\end{equation}
using $C_T\le N^{\varepsilon/4}$.
As noticed above (see \eqref{inverse_m}), $\mathfrak{m}(z)-i/t$ has only one zero at $z=it^*$, which is trivially outside $ \mathscr{S}_\zeta \backslash \mathscr{E}_{t,\varepsilon}$. The conclusion is that $\mathscr{W}(z)$ does not take the value $i/t$ on $ \mathscr{S}_\zeta \backslash \mathscr{E}_{t,\varepsilon}$, which is to say that all eigenvalues at time $t$ are in $\mathscr{E}_{t,\varepsilon} \cup \mathscr{R}_{\zeta}$.
\end{proof}

\begin{proof}[Proof of Corollary \ref{thm:hyperbolic_region}]
We apply Theorem \ref{thm:elliptic_region} and note that
\begin{equation}
\bigcup_{0<t<T} \mathscr{E}_{t,\varepsilon} \subset \mathscr{H}_{\varepsilon},
\end{equation}
which proves that for any $\varepsilon, \zeta >0$ all trajectories lie in $\mathscr{H}_{\varepsilon} \cup \mathscr{R}_{\zeta}$; choosing $\zeta < \varepsilon$ ensures that $\mathscr{R}_{\zeta} \subset \mathscr{H}_{\varepsilon}$, so that with $\mathbb{P}_H$-high probability all trajectories lie in $\mathscr{H}_{\varepsilon}$.
\end{proof}

\begin{figure}[h!]
\includegraphics[width=\textwidth]{elliptic_fig.png}
\caption{Shape of the domain \eqref{elliptic_curve} when $t$ increases (picture not to scale). The first two images correspond to some $t_1 < t_2 <1+N^{-1/3-\varepsilon}$; the last one corresponds to $t_4 > 1+N^{-1/3+\varepsilon}$, when the outlier can be isolated by applying Rouché's theorem around the upper connected component.}
\label{elliptic_fig}
\end{figure}

The above argument relied directly on the strict inequality \eqref{put_together}. For the next result, we rely on Rouché's theorem in order to isolate one particular eigenvalue, which we can do as soon as the relevant domain has a bounded connected component in $\mathscr{S}_{\zeta,2T}$, as illustrated on Fig.~\ref{elliptic_fig}.

\begin{theorem}[Emergence of an outlier]\label{thm:the_outlier}
For any $\varepsilon >0$, with $\mathbb{P}_H$-high probability, at all times $t \in ( 1+ N^{-1/3+\varepsilon}, T)$ the outlier is in the disk $D \left( it^*,
\frac{N^{\varepsilon/4}}{\sqrt{N t^*}} \right)$, whereas all other eigenvalues  satisfy $\Im \lambda_j (t)<\frac{N^{\varepsilon}}{N(t^*)^2}$. In particular they are well separated from the outlier\footnote{Recall that $t^*=t-1/t$, and so $t^* \sim t-1$ in any regime s.t. $t \sim 1$.}.
\end{theorem}

\begin{proof}
We consider some $\varepsilon'>0$ such that $\varepsilon'< \varepsilon/2$, and the domain $\mathscr{E}_{t,\varepsilon'}$ similarly as in \eqref{elliptic_curve}.
We go through the same steps as in the proof of Theorem \ref{thm:elliptic_region} and note that the inequalities \eqref{up}, \eqref{mi}, \eqref{last}, and therefore also \eqref{put_together} hold on $\mathscr{S}_{\zeta,2T} \backslash \mathscr{E}_{t, \varepsilon'}$, allowing us to invoke Rouché on any Jordan curve inside this domain. \\ \\
We further note that, for $t> 1+ N^{-1/3+\varepsilon}$ and $N$ large enough, the domain $\mathscr{E}_{t, \varepsilon'}$ has two connected components, which is a direct calculation. Let us prove that one connected component lies in the disk $D \left( it^*,
\frac{N^{\varepsilon/2}}{\sqrt{N t^*}} \right)$. First, $it^* \in \mathscr{E}_{t, \varepsilon'}$ by inspection. Then, note that $t>1+N^{-1/3+\varepsilon}$ implies $t^* > N^{-1/3+\varepsilon}$, so that for $N$ large enough,
\begin{equation}
\frac{N^{\varepsilon/4}}{\sqrt{Nt^*}} < \frac13 {t^*},
\end{equation}
which implies that for any $z=E+i\eta \in \partial D(it^*, \frac{N^{\varepsilon/4}}{\sqrt{N t^*}})$,
\begin{equation}
\eta \geq t^* - \frac{N^{\varepsilon/4}}{\sqrt{Nt^*}} > \frac23 t^*
\end{equation}
and so, for any point on that disk,
\begin{equation}
E^2 + (\eta-t^*)^2 = \frac{N^{\varepsilon/2}}{Nt^*} > \frac{N^{\varepsilon'}}{N \eta}.
\end{equation}
This and other direct considerations show that $\partial D(it^*, \frac{N^{\varepsilon/4}}{\sqrt{N t^*}}) \subset \mathscr{S}_{\zeta,2T} \backslash \mathscr{E}_{t, \varepsilon'}$. We can apply Rouché on this circle. Owing to the fact that $\mathfrak{m} - i/t$ has only one zero at $it^*$, this proves the first statement about the outlier. \\ \\
The second statement follows from checking that the second connected component is below height $\frac{N^{\varepsilon}}{N (t^*)^2}$. First, as $t^*> N^{-1/3+\varepsilon}$, we have
\begin{equation}
t^* - \frac{N^{\varepsilon}}{N (t^*)^2} > \frac23 t^*.
\end{equation}
Therefore, for any $z=E+i \frac{N^{\varepsilon}}{N (t^*)^2}$,
\begin{equation}
\eta E^2 + \eta (\eta -t^*)^2 > \frac49 (t^*)^2 \frac{N^{\varepsilon}}{N (t^*)^2} > \frac{N^{\varepsilon'}}{N}
\end{equation}
for any $\varepsilon'<\varepsilon$. This is valid outside $\mathscr{R}_{\zeta}$ for some $\zeta>0$; the second statement follows by choosing $\zeta< \varepsilon$.
\end{proof}

We finally present two complements of our main theorem, before and after the timescale at which the outlier can be isolated. The bounds we obtain for small $t$, up to slightly below the relevant $1+N^{-1/3}$ timescale, are given in Proposition \ref{prop:small_t}, whereas the bounds for $t \geq T$
are given in Proposition \ref{prop:large_t}.

\begin{proposition}[Small $t$ bounds]\label{prop:small_t}
For any $\zeta,\varepsilon>0$, with $\mathbb{P}_H$-high probability, at all times $t < 1+N^{-1/3-\varepsilon} $, all eigenvalues satisfy $\Im \lambda_j(t) < N^{-1/3 + \varepsilon}$. Moreover, for $t < 1- N^{-1/3 + \varepsilon}$,
\begin{equation}\label{bound2}
\forall j \quad \Im \lambda_j(t) < \max \left(\frac{N^{\varepsilon}}{N (t^*)^2},\frac{N^{\zeta}}{N} \right).
\end{equation}
\end{proposition}

\begin{proof}
These bounds are direct consequences of Theorem \ref{thm:elliptic_region} applied for a well chosen $\varepsilon'$. For instance, if $t < 1$, then $t^* <0$ and the inequality in \eqref{elliptic_curve} implies
$ \eta^3 < \eta (\eta-t^*)^2 < N^{-1+\varepsilon'},$
from which the bound $\Im \lambda_j < N^{-1/3+\varepsilon}$ follows if $\varepsilon'<\varepsilon$. It is a calculus exercise to check that this inequality still holds as long as $1 \leq t < 1+N^{-1/3-\varepsilon}$, when choosing an appropriate $\varepsilon'$. \\

If $t< 1-N^{-1/3+\varepsilon}$, it can be directly checked that the domain \eqref{elliptic_curve} is connected and contains the origin, so that any line that is not contained in it actually bounds it. Together with the fact that in that regime,
$$
\frac{N^{\varepsilon}}{N (t^*)^2} \ll |t^*|
$$
one can check that the line $\{ \eta = \frac{N^{\varepsilon}}{N (t^*)^2} \}$ is not in $\mathscr{E}_{t, \varepsilon}$, and therefore all eigevalues are below this threshold, or in $\mathscr{R}_{\zeta}$.
\end{proof}

\begin{proposition}[Large $t$ bounds]\label{prop:large_t}
For any $\zeta>\varepsilon>0$, with $\mathbb{P}_H$-high probability, at all times $t \in [T,N^{99}]$, the outlier is in the
small disk $D(it^*, N^{-1/2+\varepsilon})$, while all other eigenvalues are in $\mathscr{R}_{\zeta}$.
\end{proposition}

\begin{proof}
The only difference with the previous proofs is the change of domain, from a bounded one close to the real line,
to a domain
far away from the real line. There are two
consequences of this change: on the one hand, the isotropic law gives a better bound, but on the other hand the inequality \eqref{mi}
has to be replaced by a weaker one. For this, we apply the inequality \eqref{lipstick} this to $z_2 = it^*$, $z_1=E+i\eta$ with $\eta \in [t/2,2t]$,
\begin{equation}
|z - it^*| < \left|\mathfrak{m}(z) - \frac{i}{t}\right| \left(1 + \frac{t}{|\mathfrak{m}(z)|} \right)
< N^{\varepsilon/2} \eta^2 \left|\mathfrak{m}(z) - \frac{i}{t}\right|.
\end{equation}
In the domain
\begin{equation}
\mathscr{F}_{t, \varepsilon} := \{ |E|<3, \ t/2 < \eta < 2t, \ |z-it^*| \geq  N^{-1/2+\varepsilon} \},
\end{equation}
which is a rectangle with a small disk removed, we have the sequence of inequalities:
\begin{equation}
|\mathscr{W}(z) - \mathfrak{m}(z)|
< \frac{N^{\varepsilon/2}}{\eta^2 \sqrt{N}}
\leq N^{-\varepsilon/2} \frac{|z-it^*|}{\eta^2}
< \left| \mathfrak{m}(z) - \frac{i}{t} \right|
\end{equation}
and so we can apply Rouché on $\partial D(it^*, N^{-1/2+\varepsilon})$. We recall that the comparison function $\mathfrak{m} - i/t$ has exactly one root, which is at $it^*$. This proves that $\mathscr{W}-i/t$ also has exactly one root inside the disk, which is the outlier. As for the other eigenvalues, the argument from Theorem \ref{thm:the_outlier} works the same; as we assume $\zeta> \varepsilon$, the region $\mathscr{R}_{\zeta}$ contains the second connected component of \eqref{elliptic_curve}. \color{black}.
\end{proof}

\vspace{-.3in}

\section*{Acknowledgment}

We would like to thank Paul Bourgade, Victor Dubach and Yan Fyodorov for many useful remarks. \medskip

\noindent G. Dubach gratefully acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk{\l}odowska-Curie Grant Agreement No. 754411.

\begin{bibdiv}
\begin{biblist}

\bib{RMT_handbook}{book}{
  title={The Oxford Handbook of Random Matrix Theory},
  author={Akemann, G.},
  author={Baik, J.},
  author={Di Francesco, P.},
  year={2011},
  publisher={Oxford University Press}
}

\bib{AltErdoesKruegerSchroeder}{article}{
  title={Correlated random matrices: band rigidity and edge universality},
  author={Alt, J.},
  author={Erd{\H{o}}s, L.},
  author={Kr{\"u}ger, T.},
  author={Schr{\"o}der, D.},
  journal={The Annals of Probability},
  volume={48},
  number={2},
  pages={963--1001},
  year={2020},
  publisher={Institute of Mathematical Statistics}
}

\bib{BBP_paper}{article}{
  title={Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices},
  author={Baik, J.},
  author={Ben Arous, G.},
  author={P{\'e}ch{\'e}, S.},
  journal={The Annals of Probability},
  volume={33},
  number={5},
  pages={1643--1697},
  year={2005},
  publisher={Institute of Mathematical Statistics}
}

\bib{BourgadeDubach}{article}{
  title={The distribution of overlaps between eigenvectors of Ginibre matrices},
  author={Bourgade, P.},
  author={Dubach, G.},
  journal={Probability Theory and Related Fields},
  pages={1--68},
  year={2019},
  publisher={Springer}
}

\bib{Cipolloni2020}{article}{
  title={Functional central limit theorems for Wigner matrices},
  author={Cipolloni, G.},
  author={Erd{\H{o}}s, L.},
  author={Schr{\"o}der, D.},
  eprint={arXiv:2012.13218},
  year={2020}
}

\bib{ErdoesKruegerSchroeder}{article}{
  title={Random matrices with slow correlation decay},
  author={Erd{\H{o}}s, L.},
  author={Kr{\"u}ger, T.},
  author={Schr{\"o}der, D.},
  booktitle={Forum of Mathematics, Sigma},
  volume={7},
  year={2019},
  organization={Cambridge University Press}
}

\bib{ErdoesSchroeder2016}{article}{
  title={Fluctuations of functions of Wigner matrices},
  author={Erd{\H{o}}s, L.},
  author={Schr{\"o}der, D.},
  journal={Electronic Communications in Probability},
  volume={21},
  pages={1--15},
  year={2016},
  publisher={Institute of Mathematical Statistics and Bernoulli Society}
}

\bib{FyodorovPrivateCom}{article}{
  author={Fyodorov, Y. V.},
  journal={private communication},
  year={2021}
}

\bib{FyodorovKhoruzhenko}{article}{
  title={Systematic analytical approach to correlation functions of resonances in quantum chaotic scattering},
  author={Fyodorov, Y. V.},
  author={Khoruzhenko, B. A.},
  journal={Physical Review Letters},
  volume={83},
  number={1},
  pages={65},
  year={1999},
  publisher={APS}
}

\bib{FyodorovMehlig}{article}{
  title={Statistics of resonances and nonorthogonal eigenfunctions in a model for single-channel chaotic scattering},
  author={Fyodorov, Y. V.},
  author={Mehlig, B.},
  journal={Physical Review E},
  volume={66},
  number={4},
  year={2002},
  publisher={APS}
}

\bib{FyodorovSommers1996}{article}{
  title={Statistics of S-matrix poles in few-channel chaotic scattering: Crossover from isolated to overlapping resonances},
  author={Fyodorov, Y.V.},
  author={Sommers, H.J.},
  journal={Journal of Experimental and Theoretical Physics Letters},
  volume={63},
  number={12},
  pages={1026--1030},
  year={1996},
  publisher={Springer}
}

\bib{FyodorovSommers1997}{article}{
  title={Statistics of resonance poles, phase shifts and time delays in quantum chaotic scattering: Random matrix approach for systems with broken time-reversal invariance},
  author={Fyodorov, Y.V.},
  author={Sommers, H.J.},
  journal={Journal of Mathematical Physics},
  volume={38},
  number={4},
  pages={1918--1981},
  year={1997},
  publisher={American Institute of Physics}
}

\bib{FyodorovSommers}{article}{
  title={Random matrices close to Hermitian or unitary: overview of methods and results},
  author={Fyodorov, Y.V.},
  author={Sommers, H.J.},
  journal={Journal of Physics A: Mathematical and General},
  volume={36},
  number={12},
  pages={3303},
  year={2003},
  publisher={IOP Publishing}
}

\bib{HeKnowles2016}{article}{
  title={Mesoscopic eigenvalue statistics of Wigner matrices},
  author={He, Y.},
  author={Knowles, A.},
  journal={Ann. Appl. Probab.},
  volume={27},
  number={3},
  pages={1510-1550},
  year={2017}
}

\bib{KnowlesYin2013}{article}{
  title={The isotropic semicircle law and deformation of Wigner matrices},
  author={Knowles, A.},
  author={Yin, J.},
  journal={Communications on Pure and Applied Mathematics},
  volume={66},
  number={11},
  pages={1663--1749},
  year={2013},
  publisher={Wiley Online Library}
}

\bib{Peche_BBP_GUE}{article}{
  title={The largest eigenvalue of small rank perturbations of Hermitian random matrices},
  author={P{\'e}ch{\'e}, S.},
  journal={Probability Theory and Related Fields},
  volume={134},
  number={1},
  pages={127--173},
  year={2006},
  publisher={Springer}
}

\end{biblist}
\end{bibdiv}

\end{document}


\section{Method}
In this section, we first describe our proposed Discretized integrated gradients (DIG{}) and the desirable explanation axioms satisfied by it. Then we describe an interpolation algorithm that leverages our DIG{} in discrete textual domains.
Please refer to Appendix \ref{sec:background} for a brief introduction of the attribution-based explanation setup and the integrated gradients method.

\subsection{Discretized integrated gradients}
\label{sec:dig}
Below, we define our DIG{} formulation that allows interpolations along non-linear paths:
\begin{equation}
\label{eq:dig_intergral}
\textrm{ DIG{} }_{i}(x) = \int_{x^k_i=x^{\prime}_i}^{x_i} \frac{\partial F\left(x^k\right)}{\partial x_{i}} d x^k_i.
\end{equation}
Here, $x^k_i$ refers to the $i^{th}$ dimension of the $k^{th}$ interpolated point between input $x$ and  baseline $x^{\prime}$. The only constraint on $x^k_i$'s is that each interpolation should be monotonic between $x_i$ and $x^{\prime}_i$, i.e., $\forall j,k \in \{1, ..., m\}$; j < k,
\begin{equation}
\begin{split}
\label{eq:monotonic_constraints}
&x^{\prime}_i \leq x^j_i \leq x^k_i \leq x_i \textrm{~~if~~} x^{\prime}_i \leq x_i, \\
&x^{\prime}_i \geq x^j_i \geq x^k_i \geq x_i \textrm{~~otherwise}.
\end{split}
\end{equation}
This is essential because it allows approximating the integral in Eq. \ref{eq:dig_intergral} using Riemann summation\footnote{\url{https://en.wikipedia.org/wiki/Riemann\_sum}} which requires monotonic paths. We note that the interpolation points used by IG naturally satisfy this constraint since they lie along a straight line joining $x$ and $x^{\prime}$. The key distinction of our formulation from IG is that DIG{} is agnostic of any fixed step size parameter $\alpha$ and thus allows non-linear interpolation paths in the embedding space.
The integral approximation of DIG{} is defined as follows:
\begin{equation}
\label{eq:dig_sum}
\textrm{DIG{}}^{\textrm{approx}}_{i}(x) = \Sigma_{k=1}^{m} \frac{\partial F\left(x^k\right)}{\partial x_{i}} \times \big(x^{k+1}_i - x^k_i\big),
\end{equation}
where $m$ is the total number of steps considered for the approximation.

\subsection{Axioms satisfied by DIG{}}
\label{sec:axioms}
As described in prior works \cite{sundararajan2017axiomatic,shrikumar2017learning}, a good explanation algorithm should satisfy certain desirable axioms which justify the use of the algorithm for generating model explanations. Similar to IG, DIG{} also satisfies many such desirable axioms. First, DIG{} satisfies \textit{Implementation Invariance} which states that attributions should be identical for two \textit{functionally equivalent} models. Two models are functionally equivalent if they have the same output for the same input, irrespective of any differences in the model's internal implementation design. Further, DIG{} satisfies  \textit{Completeness} which states that the sum of the attributions for an input should add up to the difference between the output of the model at the input and the baseline, i.e., $\sum_i \textrm{ DIG{} }_{i}(x) = F(x) - F(x^{\prime})$. This ensures that if $F(x^{\prime}) \approx 0$ then the output is completely attributed to the inputs. Thirdly, DIG{} satisfies \textit{Sensitivity} which states that attributions of inputs should be zero if the model does not depend (mathematically) on the input.
Please refer to Appendix \ref{sec:path} for further comparisons of DIG{} with IG.

\subsection{Interpolation algorithm}
Here, we describe our proposed interpolation algorithm that searches for intermediate interpolation points between the input word embedding and the baseline embedding. Once we have the desired interpolation points, we can use Equation \ref{eq:dig_sum} to compute the word attributions similar to the IG algorithm. Please refer to Section \ref{sec:background_ig} for more details about application of IG to text.

\paragraph{Design Consideration.} First, we discuss the key design considerations we need to consider of our interpolation algorithm. Clearly, our interpolation points need to satisfy the monotonicity constraints defined in Equation \ref{eq:monotonic_constraints} so that we can use the Riemann sum approximation of DIG{}. Hence, we need to ensure that every intermediate point lies in a monotonic path. Also, the interpolation points should lie close to the original words in the embedding space to ensure that the model gradients faithfully define the model behavior.
\\

Now, we define the notion of \textit{closeness} for our specific use-case of explaining textual models. To calculate how far the interpolated words are from some true word embedding in the vocabulary, we can compute the distance of the interpolated point from the nearest word in the vocabulary. We define this as the word-approximation error (WAE). More specifically, if $w^k$ denotes the $k^{th}$ interpolation point for a word $w$, then its word-approximation error along the interpolated path is defined as:
\begin{equation}
\label{eq:word_approx_error}
\textrm{WAE}_w = \frac{1}{m}\sum_{k=1}^{m} \min_{x \in V}dist(w^k - x),
\end{equation}
where $V$ is the embedding matrix of all the words in the vocabulary. WAE of a sentence is the average WAE of all words in the sentence. Intuitively, minimizing WAE will ensure that the interpolated points are close to some real word embedding in the vocabulary which in turn ensures that output gradients of $F$ are not computed for some out-of-distribution unseen embedding points.

We observe that to minimize WAE without the monotonic constraints defined in Section \ref{sec:dig}, one can define some heuristic to search for interpolation points that belong to the set $V$ (i.e., select words from the vocabulary as interpolation points), leading to a zero WAE. Motivated by this, for a given input word embedding, we first search for an anchor word from the vocabulary that can be considered as the next interpolation point. Since the anchor point need not be monotonic w.r.t. the given input, we then optimally perturb the dimensions of the anchor word so that they satisfy the monotonicity constraints in Equation \ref{eq:monotonic_constraints}. This perturbed point becomes our first interpolation. For subsequent interpolation points, we repeat the above steps using the previous anchor and perturbed points. Formally, we break our interpolation algorithm into two parts:
\begin{enumerate}[label=(\roman*)]
	\item \textsc{AnchorSearch}: In this step, given the initial word embedding $w$, we search for an anchor word embedding $a \in V$.

	\item \textsc{Monotonize}: This step takes the anchor embedding $a$ and modifies its dimensions to create a new embedding $c$ such that all dimensions of $c$ are monotonic between the input $w$ and the baseline $w^{\prime}$.
\end{enumerate}
Overall, given an initial input word embedding $w$ and a baseline embedding $w^{\prime}$, our interpolation algorithm interpolates points from $w$ to $w^{\prime}$ (which is in decreasing order of $k$ in Eq. \ref{eq:dig_sum}). It proceeds by calling \textsc{AnchorSearch} on $w$ to get an anchor word $a$. Then, it applies \textsc{Monotonize} on $a$ to get the monotonic embedding $c$. This is our first interpolated point (in reverse order), i.e., $c = w^{m-1}$. Now, the $a$ becomes the new $w$ for the next iteration and the process continues till $m$ steps. Next, we describe in detail our specific formulations of the \textsc{Monotonize} and \textsc{AnchorSearch} algorithms.
\\

\noindent\textsc{\textbf{Monotonize}}: In this step, given an anchor word embedding $a$, we modify the non-monotonic dimensions of $a$ such that they become monotonic w.r.t. $w$ and $w^{\prime}$. The monotonic dimensions of a vector $a$ is given by:
\begin{align*}
M_a &= \{j~|~ w^{\prime}_j \leq a_j \leq w_j, j \in \{1, ..., D\} \}\\
&\cup\{j ~|~ w^{\prime}_j \geq a_j \geq w_j, j \in \{1, ..., D\} \},
\end{align*}
where $D$ is the word embedding dimension. The number of monotonic dimensions is given by the size of the set defined as $|M_a|$. Thus, the non-monotonic dimensions $\overline{M_a}$ is the set complement of the monotonic dimensions, i.e., $\overline{M_a} = \{1, ..., D\} - M_a$, where the subtraction is the set-diff operation. Let the final monotonic vector be $c$. We define the \textsc{Monotonize} operations as follows:
\begin{align*}
c[M_a] &\leftarrow a[M_a],\\
c[\overline{M_a}] &\leftarrow w[\overline{M_a}] - \frac{1}{m} \times (w[\overline{M_a}] - w^{\prime}[\overline{M_a}]),
\end{align*}
where $m$ is the total number of interpolation points we want to select in the path. It can be easily seen that $c$ is monotonic w.r.t. $w$ and $w^{\prime}$ according to the definition in Equation \ref{eq:monotonic_constraints}.
\\

\noindent\textsc{\textbf{AnchorSearch}}: First, we preprocess the word embedding in $V$ to find the top-$K$ nearest neighbor for each word. We consider this neighborhood for candidate anchor selection. Let us denote the $K$-neighbors for a word $w$ by $KNN_V(w)$.
We define two heuristics to search for the next anchor word: \textsc{Greedy}{} and \textsc{MaxCount}{}.

In the \textsc{Greedy}{} heuristic, we first compute the monotonic embedding corresponding to each word in the neighborhood $KNN_V(w)$ using the \textsc{Monotonize} step. Then, we select the anchor word $a$ that is closest to its corresponding monotonic embedding obtained from the above step. This can be thought of as minimizing the WAE metric for a single interpolated word. The key intuition here is to locally optimize for smallest perturbations at each iterative selection step. This heuristic is depicted in Figure \ref{fig:dig_greedy} and the algorithm is presented in Algorithm \ref{algo:dig_greedy} in Appendix.

In the \textsc{MaxCount}{} heuristic, we select the anchor $a$ as the word in $KNN_V(w)$ with the highest number of monotonic dimensions. Precisely, the anchor is given by:
\begin{equation*}
a = \argmax_{a' \in KNN_V(w)} {|M_{a'}|}.
\end{equation*}
The intuition of this heuristic is that the vector with highest number of monotonic dimensions would require the minimum number of dimensions being perturbed in the \textsc{Monotonize} step and hence, would be close to a word in the vocabulary. This heuristic is depicted in Figure \ref{fig:dig_count} and the algorithm is presented in Algorithm \ref{algo:dig_count} in Appendix.
\section{Introduction}
\label{sec:intro}

In the past few years, natural language processing has seen tremendous progress, largely due to strong performances yielded by pre-trained language models \cite{devlin2019bert,radford2019language,brown2020language}. But even with this impressive performance, it can still be difficult to understand the underlying reasoning for the preferred predictions leading to distrust among end-users \cite{lipton2018themythos}.
Hence, improving model interpretability has become a central focus in the community with an increasing effort in developing methods that can explain model behaviors \cite{lime,binder2016layer,li2016understanding,sundararajan2017axiomatic,shrikumar2017learning,lundberg2017unified,murdoch2018beyond}.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.85\columnwidth]{figures/intro_fig.pdf}
	\caption{\textbf{An illustration of paths used in IG and DIG{}}. IG uses a straight line interpolation with points as depicted by green squares. In contrast, DIG{} uses a non-linear path (shown in blue) with interpolation points (red stars) lying close to words in the embedding space.}
	\label{fig:intro}
	\vspace{-0.2cm}
\end{figure}

Explanations in NLP are typically represented at a word-level or phrase-level by quantifying the contributions of the words or phrases to the model's prediction by a scalar score. These explanation methods are commonly referred as \textit{attribution-based methods} \cite{murdoch2018beyond,ancona2017towards}.
Integrated Gradients (IG) \cite{sundararajan2017axiomatic} is a prominent attribution-based explanation method used due to the many desirable explanation axioms and ease of gradient computation. It computes the partial derivatives of the model output with respect to each input feature as the features are interpolated along a straight-line path from the given input to a baseline value. For example, say we want to compute the attribution for the word ``\textit{good}'' in the sentence ``\textit{the movie was good!}'' using IG. The straight-line interpolation path used by IG is depicted in green in Figure \ref{fig:intro}. Here, the baseline word is defined as the ``\textit{<pad>}'' embedding and the green squares are the intermediate interpolation points in the embedding space.

While this method can be used for attributing inputs in both continuous (e.g., image, audio, etc.) and discrete (e.g., text, molecules, etc.) domains \cite{sundararajan2017axiomatic}, their usage in the discrete domain has some limitations. Since the interpolation is done along a straight-line path joining the input word embedding and the baseline embedding (``<pad>'' in Figure \ref{fig:intro}), the interpolated points are not necessarily representative of the discrete word embedding distribution. Specifically, let a dummy word embedding space be defined by the words represented by black dots in Figure \ref{fig:intro}. Then we can see that some of the green squares can be very far-off from any original word in the embedding space. Since the underlying language model is trained to effectively work with the specific word embedding space as input, using these out-of-distribution green interpolated samples as intermediate inputs to calculate gradients can lead to sub-optimal attributions.

\begin{figure*}
	\centering
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=0.8\columnwidth]{figures/dig_greedy.pdf}
		\caption{DIG{}-\textsc{Greedy}{}}
		\label{fig:dig_greedy}
	\end{subfigure}%
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=0.8\columnwidth]{figures/dig_count.pdf}
		\caption{DIG{}-\textsc{MaxCount}{}}
		\label{fig:dig_count}
	\end{subfigure}
	\caption{
	\textbf{Overview of paths used in DIG{} and IG}. The gray region is the neighborhood of $w$. Green line depicts the straight-line path used by IG. \textbf{Left}: In DIG{}-\textsc{Greedy}{}, we first monotonize each word in the neighborhood (red arrow) and the word closest to its corresponding monotonic point is selected as the anchor ($w_5$ since the red arrow of $w_5$ has the smallest magnitude). \textbf{Right}: In DIG{}-\textsc{MaxCount}{} we select the word with the highest number of monotonic dimensions (count shown in $[.]$) as the anchor word ($w_4$), followed by changing the non-monotonic dimensions of $w_4$ (red arrow to $c$). Repeating this iteratively gives the non-linear blue path for DIG{} with the red stars as interpolation points. Please refer to Section \ref{sec:dig} for more details. Figure best viewed in color.}
	\label{fig:dig}
\end{figure*}

To mitigate these limitations, we propose a Discretized integrated gradients (DIG{}) formulation by relaxing the constraints of searching for interpolation points along a straight-line path. Relaxing this linear-path constraint leads to a new constraint on the interpolation paths in DIG{} that points along the path should be monotonically situated between the input word embedding and the baseline embedding. Hence, in DIG{}, our main objective is to monotonically interpolate between the input word embedding and baseline such that the intermediate points are close to real data samples. This would ensure that the interpolated points are more representative of the word embedding distribution, enabling more faithful model gradient computations. To this end, we propose two interpolation strategies that search for an optimal anchor word embedding in the real data space and then modify it such that it lies monotonically between the input word and baseline (see Fig.~\ref{fig:intro} for an illustration).

We apply DIG{} using our proposed interpolation algorithms to generate attributions for three pre-trained language models - BERT \cite{devlin2019bert}, DistilBERT \cite{sanh2020distilbert}, and RoBERTa \cite{liu2019roberta}, each fine-tuned separately on three sentiment classification datasets - SST2 \cite{socher2013recursive}, IMDB \cite{mass2011learning}, and Rotten Tomatoes \cite{rotten_tomatoes}. We find that our proposed interpolation strategies achieve a superior performance compared to integrated gradients and other gradient-based baselines on eight out of the nine settings across different metrics. Further, we also observe that on average, end-users find explanations provided by DIG{} to be more plausible justifications of model behavior than the explanations from other baselines.
\section{Acknowledgments}
Mention scholarship
\reminder{Fill this}
\section{Conclusion}
In this paper, we proposed Discretized integrated gradients (DIG{}) which is effective in explaining models working with discrete text data. Further, we proposed two interpolation strategies - DIG{}-\textsc{Greedy}{} and DIG{}-\textsc{MaxCount}{} that generate non-linear interpolation paths for word embedding space. Finally, we established the effectiveness of DIG{} over integrated gradients and other gradient-based baselines through experiments on multiple language models and datasets. We also conduct human evaluations and find that DIG{} enhances human trust on model predictions.
\section{Experimental Setup}
In this section, we describe the datasets and models used for evaluating our proposed algorithm.

\paragraph{Datasets.} The SST2 \cite{socher2013recursive} dataset has 6920/872/1821 example sentences in the train/dev/test sets. The task is binary classification into positive/negative sentiment. The IMDB \cite{mass2011learning} dataset has 25000/25000 example reviews in the train/test sets with similar binary labels for positive and negative sentiment. Similarly, the Rotten Tomatoes (RT) \cite{rotten_tomatoes} dataset has 5331 positive and 5331 negative review sentences. We use the processed dataset made available by HuggingFace Dataset library \footnote{\url{https://github.com/huggingface/datasets}} \cite{wolf-2020HuggingFace-datasets}.

\paragraph{Language Models.} We use pre-trained BERT \cite{devlin2019bert}, DistilBERT \cite{sanh2020distilbert}, and RoBERTa \cite{liu2019roberta} text classification models individually fine-tuned for SST2, IMDB, and RT datasets. The fine-tuned checkpoints used are provided by the HuggingFace library \cite{wolf-etal-2020-transformers}.

\paragraph{Evaluation Metrics.}
Following prior literature, we use the following three automated metrics:
\begin{itemize}
	\item \textbf{Log-odds} (LO) score \cite{shrikumar2017learning} is defined as the average difference of the negative logarithmic probabilities on the predicted class before and after masking the top $k\%$ words with zero padding. Lower scores are better.
	\vspace{-0.1cm}

	\item \textbf{Comprehensiveness} (Comp) score \cite{eraser} is the average difference of the change in predicted class probability before and after removing the top $k\%$ words. Similar to Log-odds, this measures the influence of the top-attributed words on the model's prediction. Higher scores are better.
	\vspace{-0.1cm}

	\item \textbf{Sufficiency} (Suff) score \cite{eraser} is defined as the average difference of the change in predicted class probability before and after keeping only the top $k\%$ words. This measures the adequacy of the top $k\%$ attributions for model's prediction.
	\vspace{-0.1cm}
\end{itemize}

Please refer to Appendix \ref{sec:eval_metrics_details} for more details about the evaluation metrics. We use $k=20\%$ in our experiments. In Appendix \ref{sec:topk} we further analyze the effect of changing top-k\% on the metrics. Additionally, we use our proposed word-approximation error (WAE) metric to compare DIG{} with IG.

\begin{table*}[]
	\centering
	\resizebox{\textwidth}{!}{%
		\begin{tabular}{lcccccccccccc}
			\toprule
			\multirow{2}{*}{\textbf{Method}}		& \multicolumn{4}{c}{\textbf{DistilBERT}}	& \multicolumn{4}{c}{\textbf{RoBERTa}} & \multicolumn{4}{c}{\textbf{BERT}}	\\
			\cmidrule(r){2-5} \cmidrule(r){6-9} \cmidrule(r){10-13}
			& LO $\downarrow$ & Comp $\uparrow$ & Suff $\downarrow$ & WAE $\downarrow$ & LO $\downarrow$ & Comp $\uparrow$ & Suff $\downarrow$ & WAE $\downarrow$ & LO $\downarrow$ & Comp $\uparrow$ & Suff $\downarrow$ & WAE $\downarrow$	\\
			\midrule
			Grad*Inp				& -0.152	& 0.068	& 0.315	& -	& -0.158	& 0.054	& 0.406	& -	& -0.801	& 0.204	& 0.398	& - \\
			DeepLIFT				& -0.077	& 0.017	& 0.372	& -	& -0.150	& 0.050	& 0.413	& -	& -0.388	& 0.096	& 0.438	& - \\
			GradShap				& -0.298	& 0.156	& 0.270	& -	& -0.290	& 0.128	& 0.338	& -	& -0.809	& 0.235	& 0.388	& - \\
			IG 						& -0.424	& 0.208	& 0.189	& 0.348	& -0.368	& \textbf{0.149}	& 0.317	& 0.677	& -0.789	& 0.203	& 0.418	& 0.305		\\
			\midrule
			DIG{}-\textsc{Greedy}{} 	& \textbf{-0.501}	& \textbf{0.257}	& \textbf{0.184}	& 0.329	& \textbf{-0.393}	& 0.148	& \textbf{0.294}	& 0.465	& \textbf{-1.056}	& \textbf{0.267}	& \textbf{0.416}	& 0.251		\\
			DIG{}-\textsc{MaxCount}{} 	& -0.467	& 0.231	& 0.190	& \textbf{0.230}	& -0.361	& 0.133	& 0.332	& \textbf{0.444}	& -0.874	& 0.237	& 0.430	& \textbf{0.178}		\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:results_rotten} Comparison of variants of DIG{} with baselines on three LMs fine-tuned on Rotten Tomatoes dataset. We observe that DIG{} outperforms the baselines on all three LMs. Please refer to Section \ref{sec:results} for more details.}
\end{table*}

\section{Results}

\subsection{Performance Comparison}
\label{sec:results}
\vspace{-0.1cm}

We compare DIG{} with four representative gradient-based explanation methods - Gradient*Input (Grad*Inp) \cite{shrikumar2016not}, DeepLIFT \cite{shrikumar2017learning}, GradShap \cite{lundberg2017unified}, and integrated gradients \cite{sundararajan2017axiomatic}. For the IMDB and RT datasets, we randomly sample a subset of 2,000 reviews from the public test sets to compare the different methods, due to computation costs. For the SST2 dataset, we use the complete set of 1,821 test sentences. The results are shown in Tables \ref{tab:results_sst}, \ref{tab:results_imdb}, and \ref{tab:results_rotten} for SST2, IMDB, and Rotten Tomatoes respectively.

\paragraph{Comparison with baselines.} First, we observe that across the nine different settings we studied (three language models per dataset), DIG{} consistently outperforms the baselines on eight of the settings. This is valid for all the metrics. We also note that the WAE metric is lower for all variants of DIG{} compared to IG. This validates that our proposed interpolation strategies for DIG{} is able to considerably reduce the word-approximation error in the interpolated paths and consistently improving performance on all three explanation evaluation metrics considered.

\paragraph{Comparison between variants of DIG{}.} Second, we observe that on average, DIG{}-\textsc{Greedy}{} performs better than DIG{}-\textsc{MaxCount}{}. Specifically, we find that DIG{}-\textsc{MaxCount}{} doesn't outperform DIG{}-\textsc{Greedy}{} by significantly large margins on any setting (while the opposite is true for one setting - RoBERTa fine-tuned on IMDB dataset). This could be because the DIG{}-\textsc{Greedy}{} strategy ensures that the monotonic point $c$ is always close to the anchor $a$ due to the locally greedy selection at each step which is not explicitly guaranteed by DIG{}-\textsc{MaxCount}{}. But overall, we do not find any specific performance trend between the two proposed variants and plan to study the influence of the embedding distribution in future works.

\paragraph{Analysis.} Finally, though we are able to achieve good reductions in WAE, we note that the WAE for our interpolation algorithms are not close to zero yet. This leaves some scope to design better interpolation algorithms in future. Moreover, we find that the average Pearson correlation between log-odds and WAE is 0.32 and the correlation is 0.45 if we consider the eight settings where we outperform IG. We discuss the correlations of all the settings in Appendix \ref{sec:correlation}. While this suggests a weak correlation between the two metrics, it is hard to comment if there is a causality between the two. This is partially because we believe selection of interpolation points should also take the semantics of the perturbed sentences into consideration, which we don't strongly enforce in our strategies. Hence, we think that constraining interpolations in a semantically meaningful way is a promising direction to explore.

\begin{figure}[t]
\vspace{-0.2cm}
	\centering
	\includegraphics[width=\columnwidth]{figures/human_eval.pdf}
	\vspace{-0.6cm}
	\caption{Result of human evaluation on DistilBERT model fine-tuned on SST2 dataset and BERT model fine-tuned on Rotten Tomatoes dataset. A lower mean rank means higher trustworthy explanation algorithm. For more details, refer to Section \ref{sec:human}}
	\label{fig:human_eval}
	\vspace{-0.2cm}
\end{figure}

\subsection{Human Evaluation}
\label{sec:human}
To further understand the impact of our algorithm on end users, we conduct human evaluations of explanations from our method and the two top baselines - IG and GradShap. We perform the study on the DistilBERT model fine-tuned on SST2 dataset and the BERT model fine-tuned on Rotten Tomatoes dataset. Further, we select the best variant of DIG{} on each dataset for explanation comparisons. First, we pick 50 sample sentences from each dataset with lengths between 5 and 25 words for easier visualizations. Then, we convert the attributions from each method into word highlights, whose intensity is determined by the magnitude of the attributions. Finally, we show the highlighted sentence and the model's predicted label to the annotators and ask them to rank the explanations on a scale of 1-3, ``1'' being the most comprehensive explanation that best justifies the prediction.

Figure \ref{fig:human_eval} shows the mean rank of each explanation algorithm across the two datasets. We find that DIG{} has a significantly lower mean rank compared to IG ($p-$value less than 0.001 and $10^{-4}$ on SST2 and Rotten Tomatoes respectively). Thus, we conclude that explanations generated by DIG{} are also trustworthy according to humans. Please refer to Appendix \ref{sec:visualizations} for visualizations and discussion on explanations generated by our methods.

\begin{table}[]
	\resizebox{\columnwidth}{!}{%
		\begin{tabular}{lccc}
			\toprule
			\multirow{1}{*}{\textbf{Method}} & \multicolumn{1}{c}{\textbf{SST2}} & \multicolumn{1}{c}{\textbf{IMDB}} & \multicolumn{1}{c}{\textbf{RT}}	\\
			\midrule
			IG						& -0.950				& -0.446				& -0.424				\\
			\midrule
			\small{DIG-\textsc{RandomAnchor}{}}			& -1.217 $\pm$ 0.024	& -0.834 $\pm$ 0.021 	& -0.474 $\pm$ 0.003	\\
			\small{DIG-\textsc{RandomNeighbor}{}}			& -1.247 $\pm$ 0.013	& -0.854 $\pm$ 0.015 	& -0.460 $\pm$ 0.010	\\
			\midrule
			DIG{} (best)		& \textbf{-1.259}		& \textbf{-0.878}		& \textbf{-0.501}		\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:ablation} Comparison of DIG{} with two ablation variants - DIG-\textsc{RandomAnchor}{} and DIG-\textsc{RandomNeighbor}{} on the DistilBERT model. We report 5-seed average log-odds score for the randomized methods. Please refer to Section \ref{sec:ablation} for more details.}
\end{table}

\subsection{Performance Analysis}
In this section, we report the ablation of \textsc{AnchorSearch} and the effect of path density on DIG{}. Please refer to Appendix \ref{sec:ablation_extended} for ablations on neighborhood size and discussions on computational complexity.

\label{sec:ablation}
\paragraph{Ablation Study on \textsc{AnchorSearch}.}
We ablate our methods with two random variants - DIG-\textsc{RandomAnchor}{} and DIG-\textsc{RandomNeighbor}{}, in which the \textsc{AnchorSearch} step uses a random anchor selection heuristic. Specifically, in DIG-\textsc{RandomAnchor}{}, the anchor is selected randomly from the complete vocabulary. Thus, this variant just ensures that the selected anchor is close to some word in the vocabulary which is not necessarily in the neighborhood. In contrast, the DIG-\textsc{RandomNeighbor}{} selects the anchor randomly from the neighborhood without using our proposed heuristics \textsc{MaxCount}{} or \textsc{Greedy}{}. The log-odds metrics of IG, the two ablations, and our best variant of DIG{} for DistilBERT fine-tuned individually on all three datasets are reported in Table \ref{tab:ablation}. We report 5-seed average for the randomized baselines. We observe that DIG-\textsc{RandomAnchor}{} improves upon IG on all three datasets. This shows that generating interpolation points close to the words in the vocabulary improve the explanation quality. Further, we observe that DIG-\textsc{RandomNeighbor}{} improves upon DIG-\textsc{RandomAnchor}{} on log-odds metric. One reason could be that the words in a neighborhood are more semantically relevant to the original word, leading to more coherent perturbations for evaluating model gradients. Finally, we observe that, on average, our proposed method is better compared to selecting a random anchor in the neighborhood. This shows that our search strategies \textsc{MaxCount}{} and \textsc{Greedy}{} are indeed helpful.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/interpolation_fig.pdf}
    \caption{Effect of increasing number of interpolation points $m$ on IG and DIG{}.}
    \label{fig:step_analysis}
\end{figure}

\begin{table}[t]
	\centering
	\resizebox{0.8\columnwidth}{!}{%
		\begin{tabular}{cccc}
			\toprule
			Factor $f$ & Log-Odds $\downarrow$ & WAE $\downarrow$ & Delta \% $\downarrow$	\\
			\midrule
			0		& -1.259 & 0.227	& 4.926	\\
			1		& -1.229 & 0.230	& 3.728	\\
			2		& -1.184 & 0.232	& 2.752	\\
			3		& -1.181 & 0.233	& 1.862	\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:density_analysis} Effect of up-sampling a path by a factor $f$ on Delta \% for DIG{} using $m=30$.}
\end{table}

\paragraph{Effect of Increasing Path Density.}
In integrated gradients, the \textit{completeness} axiom (Section \ref{sec:axioms}) is used to estimate if the integral approximation (Equation \ref{eq:ig_sum}) error is low enough. This error is denoted as the Delta \% error. If the error is high, users can increase the number of interpolation points $m$.

While DIG{} also satisfies the \textit{completeness} axiom, error reduction by increasing $m$ is infeasible. This is because increasing $m$ in Equation \ref{eq:dig_sum} implicitly changes the integral path rather than increasing the density. Hence, to achieve an error reduction in DIG{}, we up-sample the interpolation path $P = \{w, w_1, w_2, \dots, w_{m-2}, w^{\prime}\}$ with an up-sampling factor ($f$) of one as follows:
\[
\resizebox{\columnwidth}{!}{
$P_1 = \{w, \frac{w + w_1}{2}, w_1, \frac{w_1 + w_2}{2}, \dots, \frac{w_{m-2} + w^{\prime}}{2}, w^{\prime}\},$
}
\]
i.e., we insert the mean of two consecutive points to the path. This essentially doubles the density of points in the path.
Similarly, $P_2$ can be obtained by up-sampling $P_1$, etc. DIG{}($m,f=0$) refers to the standard DIG{} with no up-sampling.

Given that we have two hyperparameters $m$ and $f$ that determine the overall path density, we analyze the effect of each of these in Figure \ref{fig:step_analysis} and Table \ref{tab:density_analysis} respectively. The results are shown for DIG{}-\textsc{MaxCount}{} applied on DistilBERT model finetuned on SST2 dataset. In Figure \ref{fig:step_analysis}, we observe that as $m$ increases, the Delta \% of IG decreases as expected. But the trend is opposite for DIG{}. As discussed above, for DIG{}, the path length increases with increasing $m$, and hence, we attribute this trend to increasing difficulty in effectively approximating the integral for longer paths. Next, in Table \ref{tab:density_analysis}, we observe that as the up-sampling factor $f$ increases, the Delta \% consistently decreases. We also find that our up-sampling strategy does not increase the WAE by a significant amount with increasing $f$, which is desirable. Thus, this confirms that our up-sampling strategy is a good substitute of increasing $m$ for IG to effectively reduce the integral approximation error Delta \%. Following \citet{sundararajan2017axiomatic}, we choose a threshold of 5\% average Delta to select the hyperparameters. For more discussions, please refer to Appendix \ref{sec:density_details}.

\section{Related Works}
There has been an increasing effort in developing interpretability algorithms that can help understand a neural network model's behavior by explaining their predictions \cite{doshivelez2017rigorous,gilpin2019explaining}. Attributions are a post-hoc explanation class where input features are quantified by scalar scores indicating the magnitude of contribution of the features toward the predicted label. Explanation algorithms that generate attributions can be broadly classified into two categories - model-agnostic algorithms, like LIME \cite{lime}, Input occlusion \cite{li2016understanding}, Integrated gradients \footnote{Note that IG is strictly not a model-agnostic algorithm since it is defined for neural networks, but we still classify it as one since the scope of this work is limited to working on neural networks.}\cite{sundararajan2017axiomatic}, SHAP \cite{lundberg2017unified}, etc. and model-dependent algorithms, like LRP \cite{binder2016layer}, DeepLIFT \cite{shrikumar2017learning}, CD \cite{murdoch2018beyond}, ACD \cite{singh2018hierarchical}, SOC \cite{jin2020Towards}, etc. While the model-agnostic algorithms can be used as black-box explanation tools that can work for any neural network architecture, for the latter, one needs to understand the network's architectural details to implement the explanation algorithm. Typically, model-dependent algorithms require specific layer decomposition rules \cite{ancona2017towards,murdoch2018beyond} which needs to be defined for all the components in the model. Model-agnostic methods usually work directly with the model outputs and gradients which are universally available.

Due to the many desirable explanation axioms and ease of gradient computation, there has been several extensions of integrated gradients. For example, \citet{miglani2020investigating} study the effect of saturation in the saliency maps generated by integrated gradients. \citet{merrill2019generalized} extend integrated gradients to certain classes of discontinuous functions in financial domains. Further, \citet{jha2020enhanced} use KNNs and auto-encoders to learn latent paths for RNAs. Different from prior work, our focus here is to improve integrated gradients specifically for the discrete textual domain. While the idea of learning latent paths for text data is quite interesting, it brings a significant amount of challenge in successfully modeling such a complex latent space and hence, we leave this for future work.

\section{Preliminaries}
\label{sec:background}

\subsection{Attribution-based Explanations}
Attribution-based explanations generate a scalar score for a given input feature that indicates the contribution (or importance) of that feature towards particular label \cite{ancona2017towards}. Formally, let $x=\left[x_{1}, \ldots, x_{N}\right] \in \mathbb{R}^{N}$ be an input to a model which produces an output $y = [y_1, \ldots, y_C]$, where $C$ is the total number of labels. For a given label (usually the label predicted by the model), attribution-based explanation methods compute the contribution $R^{c}=\left[R_{1}^{c}, \ldots, R_{N}^{c}\right] \in \mathbb{R}^{N}$ of each feature.

\subsection{Integrated gradients}
\label{sec:background_ig}
Integrated gradients (IG) \cite{sundararajan2017axiomatic} for an input $x$ along the $i^{th}$ dimension is defined as follows:
\begin{equation}
\label{eq:ig_intergral}
\resizebox{\columnwidth}{!}{
$\textrm{ IG }_{i}(x) =\left(x_{i}-x_{i}^{\prime}\right) \times \int_{\alpha=0}^{1} \frac{\partial F\left(x^{\prime}+\alpha \times\left(x-x^{\prime}\right)\right)}{\partial x_{i}} d \alpha.$
}
\end{equation}
Here, $F$ is the neural network, $x^{\prime}$ is a baseline embedding, and $\alpha$ is the step size. Simply put, integrated gradients algorithm works by sampling points at a uniform spacing along a straight-line between the input and the baseline, and summing the model's gradient at the inputs for each interpolated points. To compute this integral efficiently, the authors propose a Riemann summation approximation defined below:
\begin{equation}
\small
\label{eq:ig_sum}
\resizebox{\columnwidth}{!}{
$\textrm{IG}^{\textrm{approx}}_{i}(x) = \left(x_{i}-x_{i}^{\prime}\right) \times \Sigma_{k=1}^{m} \frac{\left.\partial F\left(x^{\prime}+\frac{k}{m} \times\left(x-x^{\prime}\right)\right)\right)}{\partial x_{i}} \times \frac{1}{m},$
}
\end{equation}
where $m$ is the total number of steps considered for the approximation.

Next, we briefly describe how IG is used to explain a model's prediction which takes a sentence as input (for example, the model can be a text classification network).  Let $S = [w_0..w_n]$ be a sentence of length $n$ and $w_i$ be the $i^{th}$ word embedding of the sentence. Also, let $F$ be a text-classification model, i.e., $y = F(S)$. Then, IG calculates the attribution for each dimension of a word embedding $w_i$. The interpolation points required for Equation \ref{eq:ig_sum} are generated by linearly interpolating the word embedding between $w_i$ and a baseline word embedding (usually chosen as the pad embedding). Then, using Eq. \ref{eq:ig_sum}, the attribution for the $i^{th}$ dimension of $w$ is calculated. The final word attribution is the sum of the attributions for each dimension of the word embedding.

\section{Comparison with Integrated gradients and Path methods}
\label{sec:path}
It is easy to see that the approximation of integrated gradients is a special case of DIG{}. Note that the $k^{th}$ linear interpolation of the $i^{th}$ dimension of input $x$ for IG can be represented as:
\begin{equation}
\label{eq:interpolate}
x^k_i = x^{\prime}_i + \frac{k}{m} \times (x_i - x^{\prime}_i).
\end{equation}
Substituting Eq. \ref{eq:interpolate} in Eq. \ref{eq:dig_sum} gives us Eq. \ref{eq:ig_sum}.

\citet{sundararajan2017axiomatic} define path methods as the general form of integrated gradients that are applicable for all monotonic paths between the input and the baseline. Our DIG{} approach is a reformulation of the path method where the paths are not necessarily parameterized by $\alpha$, making it more applicable for discrete data domain. Hence, DIG{} also satisfies all the theoretical properties applicable for path methods - Implementation Invariance, Sensitivity, Linearity, and Completeness. We refer the readers to Proposition 2 in \citet{sundararajan2017axiomatic} for more technical details.

\section{Evaluation Metrics}
\label{sec:eval_metrics_details}
In this section, we redefine the evaluation metrics and state the formulations for each of them. In this work, we use the following three automated metrics:
\begin{itemize}
	\item \textbf{Log-odds} (LO) score \cite{shrikumar2017learning} is defined as the average difference of the negative logarithmic probabilities on the predicted class before and after masking the top $k\%$ features with zero padding. Given the attribution scores generated by an explanation algorithm, we select the top k\% words based on their attributions replace them with zero padding. More concretely, for a dataset with $N$ sentences, it is defined as:
	$$
	\log -\operatorname{odds}(k)=\frac{1}{N} \sum_{i=1}^{N} \log \frac{p\left(\hat{y} \mid \boldsymbol{x}_{i}^{(k)}\right)}{p\left(\hat{y} \mid \boldsymbol{x}_{i}\right)},
	$$
	where $\hat{y}$ is the predicted class, $\boldsymbol{x}_{i}$ is the $i^{th}$ sentence, and $\boldsymbol{x}_{i}^{(k)}$ is the modified sentence with top k\% words replaced with zero padding. Lower scores are better.

	\item \textbf{Comprehensiveness} (Comp) score \cite{eraser} is the average difference of the change in predicted class probability before and after removing the top $k\%$ features. Similar to Log-odds, this measures the influence of the top-attributed words on the model's prediction. It is defined as:
	$$
	\operatorname{Comp}(k)=\frac{1}{N} \sum_{i=1}^{N} p(\hat{y} \mid \boldsymbol{x}_{i}^{(k)}) - p(\hat{y} \mid \boldsymbol{x}_{i}).
	$$
	Here $\boldsymbol{x}_{i}^{(k)}$ denotes the modified sentence with top $k\%$ words deleted from the sentence. Higher scores are better.

	\item \textbf{Sufficiency} (Suff) score \cite{eraser} is defined as the average difference of the change in predicted class probability before and after keeping only the top $k\%$ features. This measures the adequacy of the top $k\%$ attributions for model's prediction. It is defined in a similar fashion as comprehensiveness, except the $\boldsymbol{x}_{i}^{(k)}$ is defined as the sentence containing only the top $k\%$ words. Lower scores are better.
\end{itemize}

\begin{figure}
	\centering
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/topk_logodds_sst2.pdf}
		\caption{Log-odds $\downarrow$}
	\end{subfigure}%
\\
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/topk_comp_sst2.pdf}
		\caption{Comprehensiveness $\uparrow$}
	\end{subfigure}%
\\
	\begin{subfigure}{.5\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/topk_suff_sst2.pdf}
		\caption{Sufficiency $\downarrow$}
	\end{subfigure}
	\caption{\label{fig:topk}Effect of changing top-k\% in log-odds, comprehensiveness, and sufficiency metric for the DistilBERT model fine-tuned on SST2 dataset.}
\end{figure}

\section{Effect of top-k in evaluation metrics}
\label{sec:topk}
In Figure \ref{fig:topk}, we visualize the effect of changing top-k\% on log-odds, comprehensiveness, and sufficiency metrics for DistilBERT model fine-tuned on the SST2 dataset. We compare the two variants of our method: DIG{}-\textsc{Greedy}{} and DIG{}-\textsc{MaxCount}{} with Integrated Gradients. We observe that our method outperforms IG for all values of k. Specifically, we note that the gap between DIG{} and IG is initially non-existent but then gradually increases with increasing k in Figure \ref{fig:topk} (a) and eventually saturates. This shows that although IG might be equally good as DIG{} at finding the top-5\% important words, the explanations from IG are significantly misaligned from true model behavior for higher top-k values.

\begin{table}[]
	\resizebox{\columnwidth}{!}{%
		\begin{tabular}{lccc}
			\toprule
			\multirow{1}{*}{\textbf{Dataset}} & \multicolumn{1}{c}{\textbf{DistilBERT}} & \multicolumn{1}{c}{\textbf{RoBERTa}} & \multicolumn{1}{c}{\textbf{BERT}}	\\
			\midrule
			SST2				& 1.00 & 0.00 & 0.42	\\
			IMDB				& 0.98 & -0.68 & 0.51	\\
			Rotten Tomatoes		& 0.21 & 0.22 & 0.23	\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:correlation} Pearson correlation between log-odds and WAE metrics for different dataset+LM settings. Please refer to Appendix \ref{sec:correlation} for more details.}
\end{table}

\section{Correlation between Log-odds and WAE}
\label{sec:correlation}
We compute the Pearson correlation between log-odds and WAE for each dataset + LM pair. For this, we consider the metric values for IG, DIG{}-\textsc{Greedy}{}, and DIG{}-\textsc{MaxCount}{} and report the correlations for each setting in Table \ref{tab:correlation}. We observe that, there is a strong correlation on average for DistilBERT. For BERT and RoBERTa we find a weak positive and negative correlation respectively.

\section{Ablation Studies}
\label{sec:ablation_extended}

\subsection{Effect of increasing path density}
\label{sec:density_details}
Here, we report the detailed analysis of the effect of increasing $m$ and $f$ in Tables \ref{tab:step_analysis} and \ref{tab:density_analysis_appendix} respectively. In Table \ref{tab:step_analysis}, we report the Log-odds score along with Delta \%. We do not note any consistent trend in Log-odds with increasing $m$ for both IG and DIG{}. The results of IG suggest that, as long as the Delta \% is \textit{sufficiently} low, decreasing Delta \% any further doesn't impact the explanations very significantly. Further, in Table \ref{tab:density_analysis_appendix}, we report the WAE metrics to emphasize that our up-sampling strategy doesn't increase the WAE by a significant amount, which is desirable. Also, we note a consistent increase (although marginally) in Log-odds with decreasing Delta \%. But per our previous observations on IG, we believe these changes do not imply a causal relation between the two.

\begin{table}[]
	\centering
	\resizebox{\columnwidth}{!}{%
		\begin{tabular}{lcccc}
			\toprule
			\multirow{2}{*}{$m$}		& \multicolumn{2}{c}{IG} & \multicolumn{2}{c}{DIG} \\
			\cmidrule(r){2-3} \cmidrule(r){4-5}
			& Log-Odds $\downarrow$ & Delta \% $\downarrow$ & Log-Odds $\downarrow$ & Delta \% $\downarrow$	\\
			\midrule
			10		& -0.984 & 8.064 & -1.252 & 2.263	 \\
			30		& -0.950 & 3.394 & -1.259 & 4.926	 \\
			100		& -0.933 & 1.235 & -1.258 & 9.849	 \\
			300		& -0.940 & 0.703 & -1.242 & 10.955	 \\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:step_analysis} Effect of increasing number of interpolation points $m$ on Delta \% for IG and DIG{}. Please refer to Appendix \ref{sec:density_details} for more details.}
\end{table}

\begin{table}[]
	\centering
	\resizebox{\columnwidth}{!}{%
		\begin{tabular}{lccc}
			\toprule
			Up-sampling factor $f$ & Log-Odds $\downarrow$ & WAE $\downarrow$ & Delta \% $\downarrow$	\\
			\midrule
			DIG{} ($m=30, f=0$)		& -1.259 & 0.227	& 4.926	\\
			DIG{} ($m=30, f=1$)		& -1.229 & 0.230	& 3.728	\\
			DIG{} ($m=30, f=2$)		& -1.184 & 0.232	& 2.752	\\
			DIG{} ($m=30, f=3$)		& -1.181 & 0.233	& 1.862	\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:density_analysis_appendix} Effect of up-sampling a path by a factor $f$ on Delta \% for DIG{}. For more details, refer to Appendix \ref{sec:density_details}.}
\end{table}

\begin{table}[]
	\centering
	\resizebox{0.75\columnwidth}{!}{%
		\begin{tabular}{lccc}
			\toprule
			\multirow{1}{*}{K} & \multicolumn{1}{c}{Log-odds $\downarrow$} & \multicolumn{1}{c}{WAE $\downarrow$} & \multicolumn{1}{c}{Delta \% $\downarrow$}	\\
			\midrule
			10		& -1.258 & 0.276 & 21.405	\\
			30		& -1.263 & 0.310 & 12.228	\\
			100		& -1.277 & 0.276 & 14.155	\\
			200		& -1.194 & 0.295 & 10.647	\\
			300		& -1.216 & 0.286 & 8.523	\\
			500		& -1.259 & 0.227 & 4.926	\\
			\bottomrule
		\end{tabular}%
	}
	\caption{\label{tab:nbrhood_analysis} Effect of increasing the neighborhood size $K$ of $KNN_V$ for DIG{}. Please refer to Appendix \ref{sec:nbrhood_analysis} for more details.}
\end{table}

\subsection{Effect of increasing neighborhood size}
\label{sec:nbrhood_analysis}
In this section, we study the effect of increasing the neighborhood size in DIG{}. The results are shown in Table \ref{tab:nbrhood_analysis}. We observe a clear decreasing trend in Delta \% with increasing neighborhood size, but there is no clear trend on Log-odds or WAE. Hence, we believe that the neighborhood size has little impact on the explanation quality, but we should still ensure sufficiently low Delta.

\subsection{Discussion on computational complexity}
In this section, we briefly discuss the computational complexity of our proposed interpolation strategies. The algorithms for DIG{}-\textsc{Greedy}{} and DIG{}-\textsc{MaxCount}{} are presented in Algorithms \ref{algo:dig_greedy} and \ref{algo:dig_count} respectively. From there, we observe that both our algorithms have a running time complexity of $\mathcal{O}(nmK)$, where $n$ is the number of words, $m$ is the number of interpolation points, and $K$ is the $KNN_V$ neighborhood size. While it is computationally feasible to parallelize the loops corresponding to $n$ and $K$, the same cannot be said for the loop corresponding to $m$ because we select the interpolation points iteratively. Although we empirically find in Section \ref{sec:density_details} that a small number of interpolation points are sufficient to calculate the explanations, we believe this bottleneck can be further tackled through efficient design of non-iterative search algorithms. We leave this for future works.

\section{Visualizations of explanations}
\label{sec:visualizations}
In this section, we present some interesting sentence visualizations based on explanations from DIG{} and IG for SST2 dataset in Figure \ref{fig:example_visuals}. We show the sentence visualization and the model's predicted sentiment for the sentence for each explanation algorithm. In the visualizations, the red highlighted words denote positive attributions and blue denotes negative attributions. That is, the explanation model suggests that the red highlighted words support the predicted label whereas the blue ones oppose (or undermine) the prediction. We observe that in many cases, DIG{} is able to highlight more plausible explanations. For example, in sentence pairs 1-7, clearly the DIG{} highlights are more inline with the model prediction. But we want to emphasize that it does not mean that our method always produces more plausible highlights. For example, for sentences 8-10, we observe that highlights from IG are more plausible than those of DIG{}. Hence, this shows that, while it could be a good exercise to visualize the attributions as a sanity check, we should rely more on automated metrics and human evaluations to correctly compare explanation algorithms.

\begin{figure*}
	\centering
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/great1.pdf}
	\end{subfigure}%
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/great2.pdf}
	\end{subfigure}%
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/great3.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/good1.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/good2.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/good3.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/good4.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/neg1.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/neg2.pdf}
	\end{subfigure}
\vspace{1em}
\\
	\begin{subfigure}{\textwidth}
		\centering
		\includegraphics[width=\columnwidth]{figures/sents/neg3.pdf}
	\end{subfigure}
\vspace{1em}
	\caption{\label{fig:example_visuals} Some example visualizations of attributions from DIG{} and IG for the DistilBERT model fine-tuned on SST2 dataset. The sentence visualization is followed by model's sentiment prediction for the sentence. Here, the red highlighted words denote positive attributions and blue denotes negative attributions. For more details, please refer to Appendix \ref{sec:visualizations}}
\end{figure*}

\begin{algorithm}
	\SetKwFunction{monotonize}{\textsc{Monotonize}}
	\SetKwFunction{distance}{Distance}
	\SetKwInOut{KwIn}{Input}
	\SetKwInOut{KwOut}{Output}

	\KwIn{Sentence $s = [w_1, w_2, \dots w_n]$, $k$-nearest neighbor graph for the vocabulary $KNN_V$, number of interpolation points $m$}
	\KwOut{Interpolations}

	$points = [\ ]_{n*m}$

	\For{$i \leftarrow 1$ \KwTo $n$}{
		\For{$j \leftarrow 1$ \KwTo $m$}{
		    $dists = \{\ \}$\\
		    \For{$k \leftarrow 1$ \KwTo $K$}{
		        $nbr \leftarrow KNN_V(w_i)[k]$\\
		        $c' \leftarrow \monotonize(nbr, w_i)$\\
		        $dists[nbr] \leftarrow \distance(nbr, c')$\\
			}
			$a \leftarrow \argmin_{a' \in dists} {dists[a']}$\\
			$c \leftarrow \monotonize(a, w_i)$\\
			$points[i, j] \leftarrow c$
		}
	}

	\KwRet{$points$}
	\caption{\label{algo:dig_greedy}DIG{}-\textsc{Greedy}{}}
\end{algorithm}

\begin{algorithm}
	\SetKwFunction{monotonize}{\textsc{Monotonize}}
	\SetKwInOut{KwIn}{Input}
	\SetKwInOut{KwOut}{Output}

	\KwIn{Sentence $s = [w_1, w_2, \dots w_n]$, $k$-nearest neighbor graph for the vocabulary $KNN_V$, number of interpolation points $m$}
	\KwOut{Interpolations}

	$points = [\ ]_{n*m}$

	\For{$i \leftarrow 1$ \KwTo $n$}{
		\For{$j \leftarrow 1$ \KwTo $m$}{
			$a \leftarrow \argmax_{a' \in KNN_V(w_i)} {|M_{a'}|}$\\
			$c \leftarrow \monotonize(a, w_i)$\\
			$points[i, j] \leftarrow c$
		}
	}

	\KwRet{$points$}
	\caption{\label{algo:dig_count} DIG{}-\textsc{MaxCount}{}}
\end{algorithm}