# A Convex Optimization Framework for Computing Robustness Margins of Kalman
Filters

Himanshu Prabhat Raktim Bhattacharya<EMAIL_ADDRESS><EMAIL_ADDRESS>
Aerospace Engineering, Texas A&M University,
College Station, TX, 77843-3141.

###### Abstract

This paper proposes a novel convex optimization framework for designing robust
Kalman filters that guarantee a user-specified steady-state error while
maximizing process and sensor noise. The proposed framework simultaneously
determines the Kalman gain and the robustness margin in terms of the process
and sensor noise. This is the first paper to present such a joint formulation
for Kalman filtering. The proposed methodology is validated through two
distinct examples: the Clohessy-Wiltshire-Hill equations for a chaser
spacecraft in an elliptical orbit and the longitudinal motion model of an F-16
aircraft.

###### Index Terms:

Kalman filtering, robustness margin, convex optimization.

##  I Introduction

The robustness margin in Kalman filtering is essential in quantifying the
filter’s performance under various uncertainties. It measures how well the
filter copes with model uncertainties, deviations in noise characteristics,
and environmental changes. A higher robustness margin reflects a Kalman
filter’s increased reliability and adaptability in dealing with system
variations and uncertainties. Specifically, the upper bounds of process and
sensor noise covariances are vital in assessing a Kalman filter’s robustness.

Consider a discrete-time dynamical system

State Dynamics: $\displaystyle\boldsymbol{x}_{k+1}$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}_{k}+\boldsymbol{B}\boldsymbol{w}_{k},$
(1a) Measurement Model: $\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=\boldsymbol{C}\boldsymbol{x}_{k}+\boldsymbol{n}_{k},$ (1b)

where $\boldsymbol{x}_{k}\in^{n}$ is the state vector,
$\boldsymbol{w}_{k}\in^{m}$ is the process noise, $\boldsymbol{y}_{k}\in^{p}$
is the measurement data, and $\boldsymbol{n}_{k}$ is the sensor noise. We
assume the process and the sensor noise to be zero mean Gaussian, i.e.,
$\mathbb{E}\left[\boldsymbol{w}_{k}\right]=\boldsymbol{0}$,
$\mathbb{E}\left[\boldsymbol{n}_{k}\right]=\boldsymbol{0}$,
$\mathbb{E}\left[\boldsymbol{w}_{k}\boldsymbol{w}^{T}_{k}\right]=\boldsymbol{Q}$,
and
$\mathbb{E}\left[\boldsymbol{n}_{k}\boldsymbol{n}^{T}_{k}\right]=\boldsymbol{R}$.

In conventional Kalman filtering, the system’s model is used to predict the
future uncertainty in the state, resulting in the prior state uncertainty.
When new measurements are available, they are combined with the prior to
obtain a posterior state uncertainty with minimum error variance, where the
error is defined as the difference between the true state and the predicted
state.

The minimum steady-state estimation error’s covariance
$\boldsymbol{\Sigma}_{\infty}$ for a given $\boldsymbol{R}$ and
$\boldsymbol{Q}$ is given by the solution of the following algebraic Riccati
equation

$\boldsymbol{A}\boldsymbol{\Sigma}_{\infty}\boldsymbol{A}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}-\boldsymbol{\Sigma}_{\infty}\\\
-\boldsymbol{A}\boldsymbol{\Sigma}_{\infty}\boldsymbol{C}^{T}(\boldsymbol{C}\boldsymbol{\Sigma}_{\infty}\boldsymbol{C}^{T}+\boldsymbol{R})^{-1}\boldsymbol{C}\boldsymbol{\Sigma}_{\infty}\boldsymbol{A}^{T}=\boldsymbol{0},$

assuming $(\boldsymbol{A},\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T})$ is
controllable and $(\boldsymbol{C},\boldsymbol{A})$ is observable. The steady-
state estimation error worsens with increased process and sensor noise.

The sensor noise is primarily dictated by its internal characteristics,
impacting its performance and cost. High-precision sensors with low noise
(covariance) are typically more expensive, leading to a tradeoff driven by
design economics.

The process noise represents the uncertainty in the model used to predict the
system’s state. This uncertainty could be due to factors like unmodeled
dynamics, external disturbances, or simplifications in the model. A larger
$\boldsymbol{Q}$ indicates higher uncertainty in the model-based predictions.

Significant research efforts have been dedicated to the realm of low-precision
sparse sensing, with a specific focus on sensor scheduling and sensor
selection [1, 2, 3, 4, 5, 6, 7]. Earlier studies addressing optimal sensing
architecture design assumed fixed sensor precision [1, 2, 3]. The challenge of
sparse sensing, coupled with sensor precision minimization, has been tackled
within the Kalman filtering approach in [4] and the
$\mathcal{H}_{2}/\mathcal{H}_{\infty}$ framework in [5, 6, 7]. Notably, these
existing algorithms lack considerations for uncertainty in process noise
statistics.

To address this gap, uncertainties in both process and sensor noise within the
Kalman filtering framework have been approached using robust methods such as
risk-sensitive filtering [8], M-estimation [9], and min-max formulation [10,
11]. However, these approaches do not provide a mechanism for error budgeting,
hindering the ability to assign user-specified tradeoffs between process and
sensor noise.

Overall performance typically mandates specific error budgets for each
subsystem when designing complex systems. In the context of Kalman filters,
the main challenge lies in achieving these target error budgets by selecting
sensors with the appropriate level of precision and identifying the maximum
process noise that the filter can effectively tolerate. In cases where the
system has a particularly stringent error budget, even the optimal error
covariance, derived from a combination of assumed process and sensor noise
characteristics, might not suffice. In general, increasing sensing precision
increases tolerance to process noise. However, this tradeoff is difficult to
determine for a given filter performance. This challenge becomes even more
pronounced in large-scale distributed systems, where determining the most
effective sensing architecture — sensor precision and placement — and ensuring
maximum robustness against process noise is complex and non-trivial. In these
scenarios, designers might inadvertently incorporate sensors with excessively
high precision or position them in locations that do not significantly enhance
the system’s robustness to process noise. Such decisions lead to an increased
system cost without yielding any tangible improvements in robustness. This
paper addresses this issue by proposing a solution as a convex optimization
formulation, offering an efficient approach to sensor selection and placement
for increased filter robustness.

The main contribution of this paper is a novel convex optimization framework
to design Kalman filters where the Kalman gain and sensor and process noise
characteristics are simultaneously determined. This is the first paper that
jointly determines the Kalman gain and the robustness margin in terms of the
process and sensor noise covariances. We consider continuous and discrete-time
systems. The new results are presented in theorems 1, 2 and corollaries 1, 2.
The paper also applies the proposed framework to two realistic aerospace
problems – orbit estimation of a satellite and state estimation of an F16
aircraft model.

##  II Technical Results

###  II-A Discrete-Time Kalman Filtering

######  Theorem 1.

The largest process and sensor noise covariance for which the Kalman filter
can meet a user-specified steady-state estimation error for a discrete-time
dynamical system is given by the solution of the following convex optimization
problem:

$\displaystyle\min_{\boldsymbol{K}\in^{n\times
p},\boldsymbol{\eta}\in^{m}_{+},\boldsymbol{\zeta}\in^{p}_{+}}\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{\lambda}$
subject to
$\displaystyle\begin{bmatrix}\boldsymbol{\Sigma}_{\infty}&(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{A}\sqrt{\boldsymbol{\Sigma}_{\infty}}&(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{B}&\boldsymbol{K}\\\
\ast&\boldsymbol{I}&\boldsymbol{0}&\boldsymbol{0}\\\
\ast&\boldsymbol{0}&\textbf{diag}\left(\boldsymbol{\eta}\right)&\boldsymbol{0}\\\
\ast&\ast&\ast&\textbf{diag}\left(\boldsymbol{\zeta}\right)\end{bmatrix}\geq\boldsymbol{0},$
(2a)

where $\gamma$ is used to weigh the relative importance of sensor and process
noise, and $\lambda\in[1,2]$ defines the suitable norm for optimizing the
sensor precision. The process and noise covariance are recovered as
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$, and
$\boldsymbol{K}$ is the Kalman gain.

###### Proof.

For the system defined in (1), the uncertainty propagation in the covariance
is given by

$\boldsymbol{\Sigma}^{-}_{k+1}=\boldsymbol{A}\boldsymbol{\Sigma}^{+}_{k}\boldsymbol{A}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T},$
(3)

where $\boldsymbol{\Sigma}^{+}_{k}$ is the posterior at time $k$ and
$\boldsymbol{\Sigma}^{-}_{k+1}$ is the prior at time $k+1$. The Joseph form of
the covariance update equation at time $k+1$ is given by

$\boldsymbol{\Sigma}^{+}_{k+1}=(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}^{-}_{k+1}(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})^{T}+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}.$
(4)

Combining, (3) and (4), we get

$\boldsymbol{\Sigma}^{+}_{k+1}=(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\left(\boldsymbol{A}\boldsymbol{\Sigma}^{+}_{k}\boldsymbol{A}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}\right)(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})^{T}\\\
+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}.$

In the steady-state when
$\boldsymbol{\Sigma}^{+}_{k+1}=\boldsymbol{\Sigma}^{+}_{k}=\boldsymbol{\Sigma}_{\infty}$,
we get

$\boldsymbol{\Sigma}_{\infty}=(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\left(\boldsymbol{A}\boldsymbol{\Sigma}_{\infty}\boldsymbol{A}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}\right)(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})^{T}\\\
+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}.$ (5)

We next relax the equality condition in (5) [12] to the following linear
matrix inequality using Schur complement [13],

$\displaystyle\begin{bmatrix}\boldsymbol{\Sigma}_{\infty}&(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{A}\sqrt{\boldsymbol{\Sigma}_{\infty}}&(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{B}&\boldsymbol{K}\\\
\ast&\boldsymbol{I}&\boldsymbol{0}&\boldsymbol{0}\\\
\ast&\boldsymbol{0}&\boldsymbol{Q}^{-1}&\boldsymbol{0}\\\
\ast&\ast&\ast&\boldsymbol{R}^{-1}\end{bmatrix}\geq\boldsymbol{0},$

which is convex in $\boldsymbol{K},\boldsymbol{Q}^{-1}$ and
$\boldsymbol{R}^{-1}$. The $\ast$ denotes symmetric terms. Noting that
$\boldsymbol{Q}$ and $\boldsymbol{R}$ are diagonal, we define new variables
$\boldsymbol{\eta}\in^{m}_{+}$ and $\boldsymbol{\zeta}\in^{p}_{+}$, such that
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$. The
largest process and sensor noise covariance for which the Kalman filter is
able to meet the user-specified steady-state estimation error is achieved by
minimizing the norms of $\boldsymbol{\eta}$ and $\boldsymbol{\zeta}$. The gain
in the Kalman update equations is given by $\boldsymbol{K}$. ∎

######  Corollary 1.

Consider the case where the desired steady-state performance is specified as
$\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq\theta$ instead of
the exact steady-state estimation error. For this performance specification,
Kalman gain with the optimal process and sensor noise covariances is obtained
by solving the following optimization problem:

$\displaystyle\min_{\boldsymbol{W}\in^{n\times
p},\boldsymbol{\eta}\in^{m}_{+},\boldsymbol{\zeta}\in^{p}_{+},\boldsymbol{Z}\in^{n\times
n}_{+},\boldsymbol{X}\in^{n\times
n}_{+}}\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{\lambda}$
subject to
$\displaystyle\begin{bmatrix}\boldsymbol{Z}&(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\boldsymbol{A}&(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\boldsymbol{B}&\boldsymbol{W}\\\
\ast&\boldsymbol{Z}&\boldsymbol{0}&\boldsymbol{0}\\\
\ast&\ast&\textbf{diag}\left(\boldsymbol{\eta}\right)&\boldsymbol{0}\\\
\ast&\ast&\ast&\textbf{diag}\left(\boldsymbol{\zeta}\right)\end{bmatrix}\geq\boldsymbol{0},$
(6a) $\displaystyle\begin{bmatrix}\boldsymbol{X}&\boldsymbol{I}\\\
\ast&\boldsymbol{Z}\\\ \end{bmatrix}\geq\boldsymbol{0},$ (6b)
$\displaystyle\textbf{tr}\left[\boldsymbol{X}\right]\leq\theta$ (6c)

where Kalman gain and steady-state error covariance are obtained as
$\boldsymbol{K}:=\boldsymbol{Z}^{-1}\boldsymbol{W}$ and
$\boldsymbol{\Sigma}_{\infty}:=\boldsymbol{Z}^{-1}$ respectively.The process
and noise covariance are recovered as
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$.

###### Proof.

Starting with the relaxation of (5), we arrive at the inequality condition,

$\boldsymbol{\Sigma}_{\infty}-(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\left(\boldsymbol{A}\boldsymbol{\Sigma}_{\infty}\boldsymbol{A}^{T}\right)(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})^{T}\\\
-(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})\left(\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}\right)(\boldsymbol{I}-\boldsymbol{K}\boldsymbol{C})^{T}-\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}\geq\boldsymbol{0}.$
(7)

Next, we perform congruent transform with
$\boldsymbol{Z}:=\boldsymbol{\Sigma}_{\infty}^{-1}$, and substitute
$\boldsymbol{W}:=\boldsymbol{Z}\boldsymbol{K}$ in (7) to obtain,

$\boldsymbol{Z}-(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\left(\boldsymbol{A}\boldsymbol{Z}^{-1}\boldsymbol{A}^{T}\right)(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})^{T}\\\
-(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\left(\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}\right)(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})^{T}-\boldsymbol{W}\boldsymbol{R}\boldsymbol{W}^{T}\geq\boldsymbol{0}.$
(8)

Applying Schur complement on (8) and substituting
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$, we get

$\displaystyle\begin{bmatrix}\boldsymbol{Z}&(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\boldsymbol{A}&(\boldsymbol{Z}-\boldsymbol{W}\boldsymbol{C})\boldsymbol{B}&\boldsymbol{W}\\\
\ast&\boldsymbol{Z}&\boldsymbol{0}&\boldsymbol{0}\\\
\ast&\ast&\textbf{diag}\left(\boldsymbol{\eta}\right)&\boldsymbol{0}\\\
\ast&\ast&\ast&\textbf{diag}\left(\boldsymbol{\zeta}\right)\end{bmatrix}\geq\boldsymbol{0}.$

The constraint
$\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq\theta$ is equivalent
to
$\boldsymbol{X}-\boldsymbol{Z}^{-1}\geq\boldsymbol{0},\textbf{tr}\left[\boldsymbol{X}\right]\leq\theta$.
By applying Schur complement, this can be written as a linear matrix
inequality,

$\displaystyle\begin{bmatrix}X&\boldsymbol{I}\\\ \ast&\boldsymbol{Z}\\\
\end{bmatrix}\geq\boldsymbol{0}.$

Minimizing the norms of $\boldsymbol{\eta}$ and $\boldsymbol{\zeta}$ subject
to (6a), (6b) and (6c) achieve the largest process and sensor noise covariance
satisfying $\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq\theta$. ∎

######  Remark 1.

Sparse Sensing: With overlapping sensing modalities, minimizing the
$\|\boldsymbol{\zeta}\|_{1}$ norm may lead to a sparse $\boldsymbol{\zeta}$
vector, which in turn results in a Kalman filtering process with a reduced set
of active sensors. If the Kalman gain $\boldsymbol{K}$ can achieve the desired
steady-state accuracy even when some sensors have $\zeta_{i}=0$ (or infinite
noise covariance), it implies that those sensors do not contribute to the
update law. Consequently, the corresponding column of the optimal
$\boldsymbol{K}$ matrix will be zero.

######  Remark 2.

Quantifying Largest Tolerable Sensor Degradation: Minimizing the
$\|\boldsymbol{\zeta}\|_{2}$ norm is recommended to optimize a sensing
architecture without necessarily focusing on sparsity. This approach leads to
the largest sensor noise covariance achievable for a given sensor set.
Practically, this means installing sensors with precisions higher than the
$l_{2}$ minimum $\boldsymbol{\zeta}$. The difference between the minimum
$\boldsymbol{\zeta}$ and the installed sensors’ actual (higher) precision of
the installed sensors can then be used to quantify the permissible degradation
in sensor quality.

###  II-B Continuous-Time Kalman-Bucy Filtering

In this section, consider the continuous-time system

State Dynamics: $\displaystyle\dot{\boldsymbol{x}}(t)$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}(t)+\boldsymbol{B}\boldsymbol{w}(t),$
(9a) Measurement Model: $\displaystyle\boldsymbol{y}(t)$
$\displaystyle=\boldsymbol{C}\boldsymbol{x}(t)+\boldsymbol{n}(t),$ (9b)

where $\boldsymbol{x}(t)\in^{n}$ is the state vector,
$\boldsymbol{w}(t)\in^{m}$ is the process noise, $\boldsymbol{y}(t)\in^{p}$ is
the measurement data, and $\boldsymbol{n}(t)$ is the sensor noise. We assume
the process and the sensor noise to be zero mean Gaussian, i.e.,
$\mathbb{E}\left[\boldsymbol{w}(t)\right]=\boldsymbol{0}$,
$\mathbb{E}\left[\boldsymbol{n}(t)\right]=\boldsymbol{0}$,
$\mathbb{E}\left[\boldsymbol{w}(t)\boldsymbol{w}^{T}(t)\right]=\boldsymbol{Q}$,
and
$\mathbb{E}\left[\boldsymbol{n}(t)\boldsymbol{n}^{T}(t)\right]=\boldsymbol{R}$.

The differential equation for the error covariance is given by

$\dot{\boldsymbol{\Sigma}}=(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}+\boldsymbol{\Sigma}(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})^{T}+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T},$

which in the steady state becomes

$(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}_{\infty}+\boldsymbol{\Sigma}_{\infty}(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})^{T}+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}=\boldsymbol{0}.$
(10)

######  Theorem 2.

The largest process and sensor noise covariance for which the Kalman filter
can meet a user-specified steady-state estimation error for a continuous-time
dynamical system is given by the solution of the following convex optimization
problem:

$\displaystyle\min_{\boldsymbol{K}\in^{n\times
p},\boldsymbol{\eta}\in^{m}_{+},\boldsymbol{\zeta}\in^{p}_{+}}\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{\lambda}$
subject to
$\displaystyle\begin{bmatrix}\textbf{sym}\left((\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}_{\infty}\right)&\boldsymbol{B}&\boldsymbol{K}\\\
\ast&-\textbf{diag}\left(\boldsymbol{\eta}\right)&\boldsymbol{0}\\\
\ast&\ast&-\textbf{diag}\left(\boldsymbol{\zeta}\right)\end{bmatrix}\leq\boldsymbol{0},$
(11a)

where $\gamma$ is used to weigh the relative importance of sensor and process
noise, $\lambda\in[1,2]$ defines the suitable norm for optimizing the sensor
precisions, and $\textbf{sym}\left((\cdot)\right):=(\cdot)+(\cdot)^{T}$. The
process and noise covariance are recovered as
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$, and
$\boldsymbol{K}$ is the Kalman gain.

###### Proof.

Starting with the inequality constraint obtained by the relaxation of (10)

$(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}_{\infty}+\boldsymbol{\Sigma}_{\infty}(\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})^{T}+\boldsymbol{K}\boldsymbol{R}\boldsymbol{K}^{T}+\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^{T}\leq\boldsymbol{0},$
(12)

we apply Schur complement to arrive at the following LMI,

$\begin{bmatrix}\textbf{sym}\left((\boldsymbol{A}-\boldsymbol{K}\boldsymbol{C})\boldsymbol{\Sigma}_{\infty}\right)&\boldsymbol{B}&\boldsymbol{K}\\\
\ast&-\boldsymbol{Q}^{-1}&\boldsymbol{0}\\\
\ast&\ast&-\boldsymbol{R}^{-1}\end{bmatrix}\leq\boldsymbol{0},$ (13)

Substitutions of
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$ and
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$ in the
inequality (13) leads to the constraint in (11).

∎

######  Corollary 2.

In the case where the target steady-state performance is defined as
$\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq\theta$ rather than a
precise steady-state estimation error, the Kalman gain with the optimal
process and sensor noise covariances is determined by solving the subsequent
optimization problem,

$\displaystyle\min_{\boldsymbol{W}\in^{n\times
p},\boldsymbol{\eta}\in^{m}_{+},\boldsymbol{\zeta}\in^{p}_{+},\boldsymbol{Z}\in^{n\times
n}_{+}}\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{\lambda}$
subject to
$\displaystyle\begin{bmatrix}\textbf{sym}\left(\boldsymbol{Z}\boldsymbol{A}-\boldsymbol{W}\boldsymbol{C}\right)&\boldsymbol{Z}\boldsymbol{B}&\boldsymbol{W}\\\
\ast&-\textbf{diag}\left(\boldsymbol{\eta}\right)&\boldsymbol{0}\\\
\ast&\ast&-\textbf{diag}\left(\boldsymbol{\zeta}\right)\end{bmatrix}\leq\boldsymbol{0},$
(14a) $\displaystyle\begin{bmatrix}\boldsymbol{X}&\boldsymbol{I}\\\
\ast&\boldsymbol{Z}\\\ \end{bmatrix}\geq\boldsymbol{0},$ (14b)
$\displaystyle\textbf{tr}\left[\boldsymbol{X}\right]\leq\theta$ (14c)

where the solution is given by
$\boldsymbol{Q}^{-1}:=\textbf{diag}\left(\boldsymbol{\eta}\right)$,
$\boldsymbol{R}^{-1}:=\textbf{diag}\left(\boldsymbol{\zeta}\right)$,
$\boldsymbol{\Sigma}_{\infty}:=\boldsymbol{Z}^{-1}$,
$\boldsymbol{K}:=\boldsymbol{Z}^{-1}\boldsymbol{W}$.

###### Proof.

The proof follows similar steps as corollary 1 starting with the relaxation of
(10). Further steps are omitted due to space restrictions. ∎

##  III Applications

This section applies the proposed theorems and corollaries to sensing problems
in spacecraft proximity operation in low earth orbit (LEO) and longitudinal
flight dynamics for F-16 aircraft.

###  III-A Spacecraft rendezvous maneuver

The Clohessy–Wiltshire-Hill(CWH) equations [14, 15] provide a simplified
representation of orbital relative motion, assuming a circular orbit for the
target and an elliptical or circular orbit for the chaser spacecraft. This
model offers a first-order approximation of the chaser’s motion within a
target-centered coordinate system and is frequently used in planning
rendezvous maneuvers in Low Earth Orbit(LEO) and formation flight dynamics
[16]. The CHW model with full state measurement is described by the following
linear time-invariant (LTI) system:

$\displaystyle\dot{\boldsymbol{x}}(t)$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}(t)+\boldsymbol{B}_{u}\boldsymbol{u}(t)+\boldsymbol{B}_{w}\boldsymbol{w}(t),$
(15a) $\displaystyle\boldsymbol{y}(t)$
$\displaystyle=\boldsymbol{C}\boldsymbol{x}(t)+\boldsymbol{n}(t),$ (15b)

where

$\displaystyle\boldsymbol{A}$ $\displaystyle={\begin{bmatrix}0&0&0&1&0&0\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ 3\omega_{ref}^{2}&0&0&0&2\omega_{ref}&0\\\
0&0&0&-2\omega_{ref}&0&0\\\ 0&0&-\omega_{ref}^{2}&0&0&0\end{bmatrix}},$
$\displaystyle\boldsymbol{B}_{u}$
$\displaystyle={\begin{bmatrix}0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
\end{bmatrix}}^{T},$ $\displaystyle\boldsymbol{B}_{d}$ $\displaystyle=B_{u},$
$\displaystyle\boldsymbol{C}$ $\displaystyle=\boldsymbol{I}^{n\times n},$

$\boldsymbol{x}:=[x\,y\,z\,\dot{x}\,\dot{y}\,\dot{z}]^{\prime}$,
$\boldsymbol{u}:=[F_{x}/m\,F_{y}/m\,F_{z}/m]^{\prime}$ are state and control
vectors respectively. Constant $\omega_{ref}:=\sqrt{\mu/a_{ref}}$ corresponds
to the angular velocity of the target spacecraft in the circular orbit where
$\mu$ is the earth’s gravitational parameter and $a_{ref}$ is orbital radius
of the target spacecraft. For a LEO target, $\omega_{ref}$ is approximately
equal to $0.00113$ $rad/s$. The state of the CHW model comprises the relative
position and velocities along different directions:

  * •

Along the radial direction of the target spacecraft motion: $(x,\dot{x}),$

  * •

Tangential to the target satellite orbit: $(y,\dot{y}),$

  * •

Direction normal to the orbital plane: $(z,\dot{z}).$

Impulsive thrust along respective directions
$\boldsymbol{F}:=[F_{x}\,F_{y}\,F_{z}]^{\prime}$ divided by the mass of the
chaser spacecraft $m$ acts as the control input. We assume small perturbation
in thrusters as the source of process noise $\boldsymbol{w}$. The sensor noise
$\boldsymbol{n}$ corrupts the full-state measurement model. The discrete-time
model is obtained by discretizing the continuous model (15) using the Tustin
method [17] with sampling time $T_{s}=0.01$ sec.

Simulation studies are performed on two cases for discrete and continuous time
systems. The cost functions considered for the optimization problems in
corollaries 1 and 2 for each case are

$J=\begin{cases}\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{2}&\text{Case
1 (c1)}\\\
\|\boldsymbol{\boldsymbol{W}_{q}\eta}\|_{2}+\gamma\|\boldsymbol{W_{r}\zeta}\|_{2}&\text{Case
2 (c2)}\\\ \end{cases}$ (16)

with $\gamma=1$. Case 1 assigns equal weight to process and sensor noise
variances during optimization. Internal weights
$W_{q}=\textbf{diag}\left(\begin{bmatrix}1&100&10\end{bmatrix}\right)$ and
$W_{r}=\textbf{diag}\left(\begin{bmatrix}100&10&1&100&10&1\end{bmatrix}\right)$
in case 2 adds varying priorities to each noise with the expectation of
increment in optimal variance for those linked with larger weight. Both cases
are constrained by steady-state filter performance criteria,
$\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq 0.1$.

We arbitrarily consider the continuous time model with parameters in case 2 to
simulate real-time filter behavior with process and sensor noise
characteristics and Kalman gain obtained from the optimization in corollary 2.
The plots of mean and variance propagation for relative position errors
$(e_{x},e_{y},e_{z})$ and velocity errors
$(e_{\dot{x}},e_{\dot{y}},e_{\dot{z}})$ respectively presented in Fig.1 and
Fig.2, demonstrate convergence of the filter.

Next, we study the effect of internal weights $\boldsymbol{W}_{q}$,
$\boldsymbol{W}_{r}$ on process and sensor noise variances by comparing
optimal outcomes for cases 1 and 2. Fig.3 and Fig.4 show bar plot comparisons
between both cases for continuous and discrete time models respectively. As
expected from weight allocation, variance of noises
$(w_{F_{y}},n_{x},n_{\dot{x}})$ in case 2 are larger than case 1.

Figure 1: Positional error evolution. Solid lines denote mean error in
relative position $(x,y,z)$, whereas shaded regions show
$\mu_{i}\pm\sqrt{\Sigma_{ii}}$ where $i\in(e_{x},e_{y},e_{z})$ Figure 2:
Velocity error evolution. Solid lines denote mean error in relative velocity
$(\dot{x},\dot{y},\dot{z})$, whereas shaded regions show
$\mu_{i}\pm\sqrt{\Sigma_{ii}}$ where
$i\in(e_{\dot{x}},e_{\dot{y}},e_{\dot{z}})$ Figure 3: Bar plot comparing
process(left) and sensor(right) noise variances between case 1(c1) and case
2(c2) for the continuous-time system. The units for
$(w_{F_{x}},w_{F_{y}},w_{F_{z}})$, $(n_{x},n_{y},n_{z})$, and
$(n_{\dot{x}},n_{\dot{y}},n_{\dot{z}})$ are $(m/s^{2})^{2}$, $m^{2}$, and
$(m/s)^{2}$ respectively. Figure 4: Comparison of process(left) and
sensor(right) noise variances between case 1(c1) and case 2(c2) for the
discrete-time system. The units for $(w_{F_{x}},w_{F_{y}},w_{F_{z}})$,
$(n_{x},n_{y},n_{z})$, and $(n_{\dot{x}},n_{\dot{y}},n_{\dot{z}})$ are
$(m/s^{2})^{2}$, $m^{2}$, and $(m/s)^{2}$ respectively.

###  III-B Flight control

Let us consider the longitudinal dynamics of F-16 aircraft linearized about
the equilibrium/trim points

$\displaystyle\boldsymbol{x}_{0}$
$\displaystyle={\begin{bmatrix}1000&&-3.02\times 10^{-3}&&-3.02\times
10^{-3}&&0\end{bmatrix}}^{T},$ $\displaystyle\boldsymbol{u}_{0}$
$\displaystyle={\begin{bmatrix}6041.2&-1.38\end{bmatrix}}^{T},$

for steady flight conditions. The LTI system for linearized longitudinal
dynamics is given by[18]

$\displaystyle\dot{\boldsymbol{x}}(t)$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}(t)+\boldsymbol{B}_{w}\boldsymbol{w}(t),$
(17a) $\displaystyle\boldsymbol{y}(t)$
$\displaystyle=\boldsymbol{C}_{y}\boldsymbol{x}(t)+\boldsymbol{n}(t),$ (17b)

where

$\displaystyle\boldsymbol{A}$ $\displaystyle={\begin{bmatrix}-1.8969\times
10^{-2}&-0.4052&-32.17&0.8915\\\ -6.4397\times 10^{-5}&-1.6176&0&0.9325\\\
0&0&0&1\\\ 0&-2.3683&0&-1.9696\end{bmatrix}},$
$\displaystyle\boldsymbol{B}_{d}$ $\displaystyle={\begin{bmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&0&1\\\ \end{bmatrix}}^{T},$ $\displaystyle\boldsymbol{C}_{y}$
$\displaystyle=\begin{bmatrix}-0.0191&-5.2893&-32.17&3.7071\\\
-0.0643&-1.6176&0.0971&932.5332\\\ 0&1&0&0\\\ 0&0&0&1\\\
1.7578&0&0&0\end{bmatrix}.$

The state vector is denoted as $\boldsymbol{x}=[V\ \alpha\ \theta\
q]^{\prime}$, representing the velocity $V$ in $ft/s$, angle of attack
$\alpha$ in radians, pitch angle $\theta$ in radians, and pitch rate $q$ in
$rad/s$. The measurement vector is $y=[\dot{u}\ \dot{w}\ \alpha\ q\
\bar{q}]^{\prime}$, including body acceleration along the roll axis $\dot{u}$
in $ft/s^{2}$, body acceleration along the yaw axis $\dot{w}$ in $ft/s^{2}$,
angle of attack $\alpha$ in radians, pitch rate $q$ in rad/s, and dynamic
pressure $\bar{q}:=\rho_{atm}V^{2}/2$ in $lb/ft^{2}$, where $\rho_{atm}$
defines atmospheric density at an altitude corresponding to the steady flight.
Similar to the previous example, we consider two cases for the continuous time
system where cost functions are given by (16). In case 1(c1), where
$\gamma=1$, equal weightage is assigned to both process and sensor noise
variances during optimization. Case 2(c2) assigns internal weights
$W_{q}=\textbf{diag}\left(\begin{bmatrix}1&10&1\end{bmatrix}\right)$ to
process noise and
$W_{r}=\textbf{diag}\left(\begin{bmatrix}1&1&0.1&1&1\end{bmatrix}\right)$ to
sensor noise. The steady-state filter performance constraint
$\textbf{tr}\left[\boldsymbol{\Sigma}_{\infty}\right]\leq 0.1$ is imposed on
both cases.

We investigate the influence of internal weights, denoted as
$\boldsymbol{W}_{q}$ and $\boldsymbol{W}_{r}$, on the process and sensor noise
variances by comparing the optimal results between cases 1 and 2. Fig.5 shows
bar plot comparison between both cases. As anticipated from the allocation of
weights, there is an observed increase in the variance of the process noise
$w_{\alpha}$ and a decrease in the variance of the sensor noise $n_{\alpha}$
in case 2. The internal weights $\boldsymbol{W}_{q}$ and $\boldsymbol{W}_{r}$
shape the optimization outcome, influencing the tradeoff between process and
sensor noise variances. This optimization framework facilitates adding
constraints that explicitly specify lower (upper) bounds on the variance
(precision) of the noises, i.e.

$\eta_{i}\leq\eta_{max},\zeta_{i}\leq\zeta_{max}.$

Adding these constraints may lead to infeasibility, which can be addressed by
tweaking the performance constraints.

Next, we consider the case of sparse sensing, which is simulated by using the
cost function $J:=\|\boldsymbol{\eta}\|_{2}+\gamma\|\boldsymbol{\zeta}\|_{1}$
with $\gamma=1$. Stem plot of sensor precisions $\zeta_{i}$ in Fig.6 suggests
a sparse sensor configuration by eliminating sensor requirement for $\alpha$
and $q$. In some cases, we may obtain significantly small non-zero sensor
precision for specific sensors. Iteratively reweighted minimization proposed
in [5] can be implemented alongside the proposed algorithms for such cases to
obtain sparse sensing architectures while maximizing process noise variance.

Figure 5: Bar plot comparing process(left) and sensor(right) noise variances
between case 1(c1) and case 2(c2) for the longitudinal flight dynamics. The
units for $w_{V}$, $(w_{\alpha},n_{\alpha})$, $(w_{q},n_{q})$,
$(n_{\dot{u}},n_{\dot{w}})$ and $n_{\bar{q}}$ are $(ft/s)^{2}$, $rad^{2}$,
$(rad/s)^{2}$, $(ft/s^{2})^{2}$, and $(lb/ft^{2})^{2}$ respectively. Figure 6:
Stem plot of process noise variances(left) and sensor precisions(right). The
units for $w_{V}$, $w_{\alpha}$, $w_{q}$, $\zeta_{\alpha}$, $\zeta_{q}$,
$(\zeta_{\dot{u}},\zeta_{\dot{w}})$ and $\zeta_{\bar{q}}$ are $(ft/s)^{2}$,
$rad^{2}$, $(rad/s)^{2}$,$rad^{-2}$, $(rad/s)^{-2}$, $(ft/s^{2})^{-2}$, and
$(lb/ft^{2})^{-2}$ respectively.

##  IV Summary & Conclusion

This paper introduced a novel convex optimization framework maximizing the
robustness margin for the Kalman filter while ensuring adherence to a
predefined steady-state error budget. The formulation is detailed for both
discrete and continuous linear time-varying dynamics. The proposed algorithms
can yield an optimal sensing architecture with or without sparseness while
maximizing the process noise variance. This is demonstrated using the examples
of spacecraft rendezvous maneuvers and longitudinal dynamics of F-16 aircraft.

##  V Acknowledgements

This research is sponsored by the AFOSR Dynamic Data and Information
Processing Program, grant FA9550-22-1-0539.

## References

  * [1] Siddharth Joshi and Stephen Boyd.  Sensor selection via convex optimization.  IEEE Transactions on Signal Processing, 57(2):451–462, 2009.
  * [2] Haotian Zhang, Raid Ayoub, and Shreyas Sundaram.  Sensor selection for kalman filtering of linear dynamical systems: Complexity, limitations and greedy algorithms.  Automatica, 78:202–210, 2017.
  * [3] R. Cumbo, L. Mazzanti, T. Tamarozzi, P. Jiranek, W. Desmet, and F. Naets.  Advanced optimal sensor placement for kalman-based multiple-input estimation.  Mechanical Systems and Signal Processing, 160:107830, 2021.
  * [4] Niladri Das and R. Bhattacharya.  Optimal sensor precision for multirate sensing for bounded estimation error.  IEEE Transactions on Aerospace and Electronic Systems, 58:844–854, 2021.
  * [5] Vedang M. Deshpande and Raktim Bhattacharya.  Sparse sensing and optimal precision: Robust $\mathcal{H}_{\infty}$ optimal observer design with model uncertainty.  In 2021 American Control Conference (ACC), pages 4105–4110, 2021.
  * [6] Vedang M. Deshpande, Raktim Bhattacharya, and Kamesh Subbarao.  Sensor placement with optimal precision for temperature estimation of battery systems.  IEEE Control Systems Letters, 6:1082–1087, 2022.
  * [7] Vedang M. Deshpande and Raktim Bhattacharya.  Guaranteed robust performance of $\mathcal{H}_{\infty}$ filters with sparse and low precision sensing.  IEEE Transactions on Automatic Control, 69(2):1029–1036, 2024.
  * [8] J. Speyer, J. Deyst, and D. Jacobson.  Optimization of stochastic linear systems with additive measurement and process noise using exponential performance criteria.  IEEE Transactions on Automatic Control, 19(4):358–366, 1974.
  * [9] Z.M. Durovic and B.D. Kovacevic.  Robust estimation with unknown noise statistics.  IEEE Transactions on Automatic Control, 44(6):1292–1296, 1999.
  * [10] Lars Peter Hansen.  Uncertainty outside and inside economic models.  Working Paper 20394, National Bureau of Economic Research, August 2014.
  * [11] Bingyan Han.  Distributionally robust kalman filtering with volatility uncertainty.  2023\.
  * [12] Jan Willems.  Least squares stationary optimal control and the algebraic riccati equation.  IEEE Transactions on automatic control, 16(6):621–634, 1971.
  * [13] Fuzhen Zhang.  The Schur complement and its applications, volume 4.  Springer Science & Business Media, 2006.
  * [14] W. H. CLOHESSY and R. S. WILTSHIRE.  Terminal guidance system for satellite rendezvous.  Journal of the Aerospace Sciences, 27(9):653–658, 1960.
  * [15] G. W. Hill.  Researches in the lunar theory.  American Journal of Mathematics, 1(1):5–26, 1878.
  * [16] Joseph A. Starek, Edward Schmerling, Gabriel D. Maher, Brent W. Barbee, and Marco Pavone.  Fast, safe, propellant-efficient spacecraft motion planning under clohessy–wiltshire–hill dynamics.  Journal of Guidance, Control, and Dynamics, 40(2):418–438, 2017.
  * [17] Gene Franklin, J.D. Powell, and M.L. Workman.  Digital Control of Dynamic Systems-Third Edition.  11 2022.
  * [18] Brian L. Stevens, Frank L. Lewis, and Eric N. Johnson.  Aircraft control and simulation: Dynamics, controls design, and autonomous systems: Third edition.  wiley, November 2015.  Publisher Copyright: © 2016 by John Wiley & Sons, Inc. All rights reserved.


# Revisiting Attention Weights as Interpretations
of Message-Passing Neural Networks

Yong-Min Shin  Siqing Li  Xin Cao  Won-Yong Shin

###### Abstract

The self-attention mechanism has been adopted in several widely-used message-
passing neural networks (MPNNs) (e.g., GATs), which adaptively controls the
amount of information that flows along the edges of the underlying graph. This
usage of attention has made such models a baseline for studies on explainable
AI (XAI) since interpretations via attention have been popularized in various
domains (e.g., natural language processing and computer vision). However,
existing studies often use naïve calculations to derive attribution scores
from attention, and do not take the precise and careful calculation of edge
attribution into consideration. In our study, we aim to fill the gap between
the widespread usage of attention-enabled MPNNs and their potential in largely
under-explored explainability, a topic that has been actively investigated in
other areas. To this end, as the first attempt, we formalize the problem of
edge attribution from attention weights in GNNs. Then, we propose GAtt, an
edge attribution calculation method built upon the computation tree. Through
comprehensive experiments, we demonstrate the effectiveness of our proposed
method when evaluating attributions from GATs. Conversely, we empirically
validate that simply averaging attention weights over graph attention layers
is insufficient to interpret the GAT model’s behavior. Code is publicly
available at https://github.com/jordan7186/GAtt/tree/main.

##  1 Introduction

###  1.1 Background

In graph learning, graph neural networks (GNNs) (Wu et al., 2021) have been
used as the de facto architecture, since they can effectively encode the graph
structure along with the node (or edge) features. Among various GNNs (Kipf &
Welling, 2017; Hamilton et al., 2017; Xu et al., 2019), several models have
successfully incorporated the self-attention mechanism (Bahdanau et al., 2015;
Vaswani et al., 2017) into the message passing neural networks (MPNNs) (Gilmer
et al., 2017; Bronstein et al., 2021), e.g., GAT (Velickovic et al., 2018).
Such architectures have effectively allowed to adaptively weigh the incoming
information during message passing.

Figure 1: A visualization of two different edge attribution methods given the
attention weights in a GAT model, a representative GNN architecture with self-
attention, trained on the Infection dataset. Contrast to simply averaging over
graph attention layers, our proposed method, GAtt, is built upon the
computation tree, and reveals that the GAT model focuses mostly on the edges
included in the ground truth explanation (blue edges in the input graph).

Despite their widespread usage, GNNs, similarly as in many other neural
networks, are considered as black-box models that lack interpretability. To
shed light on the inner workings of such GNNs, a large number of studies have
focused on developing new post-hoc explainability methods, i.e., methods that
are applied after the underlying GNN model completed its training (Li et al.,
2022; Yuan et al., 2023). However, such approaches typically necessitate an
optimization framework, making the explanation performance dependent on
various hyperparameters (e.g., number of iterations and random seeds). Such
post-hoc explanations are also known to be sub-optimal since they essentially
perform a single-step projection to an information-controlled space (Miao et
al., 2022). This motivates us to explore an exciting opportunity for GNNs that
produces attention weights (dubbed Att-GNNs), as it would be beneficial to
introduce a new method that directly utilizes the attention weights for
explanation without any additional optimization/learning frameworks and
hyperparameters. Furthermore, attention itself is also viewed as a direct way
of providing model interpretations without a separate post-hoc explanation
method (Lee et al., 2017; Ghaeini et al., 2018; Hao et al., 2021). This
viewpoint has already been extensively performed on transformers, the
representative model architecture employing attention (Bahdanau et al., 2015;
Xu et al., 2015; Vig, 2019; Dosovitskiy et al., 2021; Caron et al., 2021).
Also, there has been a significant body of research debating the effectiveness
of attention as explanation, particularly within the domain of natural
language processing (NLP) (Jain & Wallace, 2019; Wiegreffe & Pinter, 2019;
Bibal et al., 2022).

###  1.2 Motivation

In contrast to other domains and architectures, the exploration of the
attention weight’s potential as explanation has been largely overlooked
outside transformer architectures. While a few studies incorporated GATs as a
baseline in their evaluations (Ying et al., 2019; Luo et al., 2020; Sánchez-
Lengeling et al., 2020), they used attention as edge attribution (i.e.,
highlighting relevant edges of the input graph) by simply averaging the
attention weights over different layers, the precise computation of edge
attribution from attention warrants a more intricate methodology.

Our study seeks to formulate post-processing for attention weights in Att-
GNNs, enabling to extract higher-quality edge attributions that can better
capture the essence of Att-GNNs.111As mentioned, although graph transformers
(Ying et al., 2021; Kreuzer et al., 2021; Chen et al., 2023) also heavily
utilize attention, we focus on Att-GNNs in this study due to their widespread
usage as well as their message passing-based architecture. While analogous
investigations have been conducted in the context of transformers (Abnar &
Zuidema, 2020; Chefer et al., 2021b, a), to the best of our knowledge, we are
the first to address this issue within the domain of Att-GNNs (i.e., MPNNs
using attention).

###  1.3 Main Contributions

In this study, we tackle the intricate challenge of attributing edges using
attention weights in Att-GNNs by dissecting the Att-GNNs feed-forward process,
while taking into account the model’s underlying computation tree. Drawing
insights from this analysis, we assert that the edge attribution function to
be designed should encompass two crucial dimensions: 1) proximity to the
target node and 2) its position in the computation tree, thus aligning with
the intrinsic feed-forward process.

To this end, we introduce GAtt, which adds attention weights across the
computation tree while adjusting their influence by employing targeted
multiplication factors for attention weights guiding towards the target node.
When using a GAT model as a representative GNN architecture, we run extensive
experiments by answering pivotal facets of interpretation—faithfulness and
explanation accuracy—of the underlying GAT across diverse real-world and
synthetic datasets. Empirical results demonstrate that the integration of GAtt
to process attention weights yields substantively enhanced explanations,
excelling in both faithfulness and explanation accuracy. As an example, Figure
1 visualizes attribution scores from different edge attribution methods using
the same GAT model. The attribution scores from GAtt (see red box) imply that
the GAT model places high emphasis on the edges included in the ground truth
explanation. Such conclusion could not have been reached if we were to use
simple averaging as the tool for interpretation, validating the superiority of
our method. Finally, we also perform an ablation study where we introduce two
variants of GAtt, namely GAtt${}_{\textsc{sim}}$ and GAtt${}_{\textsc{avg}}$,
where we each remove one critical design element of our method. Our analysis
shows that although the variants still outperforms naïve layer-wise averaging,
it deteriorates the quality of the attribution in all measures. In summary, we
conclude that Att-GNNs such as GAT are indeed highly explainable when adopting
the proper interpretation, i.e., proper adjustment of attention weights by
taking the computation tree viewpoint. Note that, although the principles of
GAtt is generally applicable to all Att-GNNs, we focus on the original GAT
model in our study due to its widespread usage and significance in the
literature. Our contributions are summarized as follows:

  * •

Key observations: We make key observations and design principles that need to
be considered in the edge attribution calculation by integrating the
computation tree of the target node during its calculation.

  * •

Novel methodology: We propose GAtt, a new method to calculate edge attribution
from attention weights in a Att-GNNs by integrating the computation tree of
the given GNN model.

  * •

Extensive evaluations: We extensively demonstrate that the GAT model is shown
to be more faithful and accurate when using our proposed method compared to
the simple alternative.

##  2 Related Work

In this section, we discuss relevant prior studies in two major themes.
Specifically, we first address previous attempts to use attention weights as
interpretations in other domains, and then we provide an overview of
explainable AI (XAI).

###  2.1 Interpreting Models with Attention

Attention weights have been used as a useful tool to visualize the inner
workings of the underlying model (Bahdanau et al., 2015; Xu et al., 2015; Vig,
2019; Dosovitskiy et al., 2021; Caron et al., 2021), and there have been
multiple studies that systematically questioned the utility of attention as
interpretation in NLP (Bibal et al., 2022). In (Jain & Wallace, 2019),
attention was found to have low correlations to other importance measures.
However, (Wiegreffe & Pinter, 2019) pointed out an unfair setting in the prior
work, and argued that attention can still be an effective explanation.
Furthermore, there have been studies that focused on post-processing attention
in transformers for token attribution. (Abnar & Zuidema, 2020) proposed
attention rollout and attention flow, and follow-up attempts (Chefer et al.,
2021a, b) were built upon attention rollout by incorporating other attribution
methods or attention gradients. Although our study lies in a similar
objective, our focus is on graph data with the Att-GNN architecture, which has
been largly under-explored in the literature.

###  2.2 Explainability in GNNs

The primary goal of XAI is to provide a comprehensible understanding of neural
network models’ decisions. In recent years, various studies have developed
methods to explain GNN models. As one of the pioneering work, GNNExplainer
(Ying et al., 2019) identified a subset of edges and node features around the
target node that affect the underlying model’s decision. PGExplainer (Luo et
al., 2020) trained a separate parameterized mask predictor to generate edge
masks that identify important edges. Although explanations of GNN models are
still an active research area (Li et al., 2022; Yuan et al., 2023), most
studies overlooked Att-GNN as an inherently explainable model. Several
attempts (Ying et al., 2019; Luo et al., 2020; Sánchez-Lengeling et al., 2020)
introduced GATs, a representitative model in Att-GNNs, as a baseline by
averaging attention over layers. In light of this, our study aims to
scrutinize attention as a paramount candidate for explaining Att-GNNs.

##  3 Edge Attribution Calculation in Att-GNNs

In this section, we first describe our notations used in the paper. Then, we
formalize the problem of calculating the edge attribution in Att-GNNs, and
then propose GAtt, an approach for incorporating the computation tree to
compute edge attributions.

###  3.1 Notations

Let us denote a given (undirected) graph as a tuple of two sets
$G=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the set of nodes and
$\mathcal{E}$ is the set of edges. We denote the edge connecting two nodes
$v_{i},v_{j}\in\mathcal{V}$ as $e_{ij}\in\mathcal{E}$. For the undirected $G$,
$e_{j,i}\in\mathcal{E}$ if $e_{i,j}\in\mathcal{E}$ and vice versa. We also
denote the set of neighbors of node $v_{i}$ as $\mathcal{N}_{i}$.

###  3.2 Problem Statement

We are given a graph $G=(\mathcal{V},\mathcal{E})$, the Att-GNN model $f$ with
$L$ layers and a target node $v_{i}\in\mathcal{V}$ of interest. The attention
weights calculated from $f$ are denoted as $\mathcal{A}=\\{{\bf
A}(l)\\}_{l=1}^{L}$, where ${\bf
A}(l)\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$ and $[{\bf
A}(l)]_{j,i}=\alpha_{i,j}^{l}$ is the attention weight of edge $e_{i,j}$ in
the $l$-th layer ($l=1$ being the input layer). The problem of edge
attribution calculation is characterized by an edge attribution function
$\Phi(v,\mathcal{A},e_{i,j})\triangleq\phi_{i,j}^{v}$ such that the edge
attribution score $\phi_{i,j}^{v}$ accounts for the importance of edge
$e_{i,j}$ to the model’s calculation for node $v$ (i.e., faithfulness of $f$).

In our study, we design $\Phi$ using the computation tree in Att-GNNs
alongside several observations and key design principles, which will be
specified later. To design such a function $\Phi$, we argue that the
computation tree of Att-GNNs should be considered, and incorporating its
several key properties. Note that most post-hoc instance-level explanation
methods for GNNs (Ying et al., 2019; Luo et al., 2020) also have a similar
objective in terms of calculating $\phi_{i,j}^{v}$, but they do not leverage
the attention weights $\mathcal{A}$. Additionally, we mainly consider GATs as
the representative Att-GNN model in our study, but our methodology can be
generally applied to other Att-GNN models. We will focus on the original GAT
model (Velickovic et al., 2018) as a default from now on.

Figure 2: A visualization for a 2-layer GAT on node 27 on the infection
dataset. LABEL:subfig:infectiondatasetexample shows the local 2-hop subgraph
with the edge $e_{40,27}$ marked as red.
LABEL:subfig:infectiondatasetcomputationgraph shows the computation tree in
the GAT, where the information flows from leaf nodes to node 27 at the root.
The edges are colored by the attention weights from the model, while
highlighting the two occurrences of edge $e_{40,27}$.

###  3.3 From Attention to Attribution

We first visualize the computation tree in a GAT, which will lead to several
observations to guide GAtt, an edge attribution calculation method given the
attention weights in the GAT model.

####  3.3.1 Visualizing the computation tree

To provide an illustrative example, we train a 2-layer GAT model (Velickovic
et al., 2018) with a single attention head on the synthetic infection
benchmark dataset (Faber et al., 2021). LABEL:subfig:infectiondatasetexample
shows the 2-hop subgraph from target node 27, which contains all nodes and
edges that the GAT model utilizes from node 27’s point of view. The
computation tree in the GAT is commonly expressed as a rooted subtree (Sato et
al., 2021), which is shown in LABEL:subfig:infectiondatasetcomputationgraph
for node 27. In the figure, the information flows from leaf nodes at depth 2
to the root node 27 at depth 0, which exhibits a different structure from that
of the subgraph in LABEL:subfig:infectiondatasetexample. Since each graph
attention layer calculates attention weights for each edge in $\mathcal{E}$,
we color-code those weights in the computation tree (note that self-attention
always includes self-loops regardless of the original graph structure).

####  3.3.2 Observations and Design Principles

We make several observations from the computation tree:

  1. (O1)

Identical edges can appear multiple times in the computation tree. For
example, edge $e_{40,27}$ in LABEL:subfig:infectiondatasetexample appears
twice in LABEL:subfig:infectiondatasetcomputationgraph.

  2. (O2)

Nodes do not appear uniformly in the computation tree. Specifically, nodes
that are $k$-hops away from the target node do not exist in depth $k^{\prime}$
for $0<k^{\prime}<k$ (e.g., node 70 appears only at depth 2 while node 40
appears three times).

  3. (O3)

The graph attention layer always includes self-loops during its feed-forward
process.

Based on the above observations, we would like to state two design principles
that must be reflected when building the edge attribution function $\Phi$.

  1. (P1)

Proximity effect: Edges within closer proximity to the target node tend to
highly impact the model’s prediction compared with faraway edges, since they
tend to appear more frequently in the computation tree.

  2. (P2)

Contribution adjustment: The contribution of an edge in the computation tree
should be adjusted by its position. (i.e., depth in the tree)

By these standards, we look again at the attention weights in
LABEL:subfig:infectiondatasetcomputationgraph. We first see that edges close
to the target node such as $e_{40,27}$ appear twice, whereas distant edges
such as $e_{70,40}$ appear only once (P1). Additionally, for the attention
weights from the last graph attention layer (i.e., edges connecting nodes at
depth 1 to the root node), each edge tends to have roughly the value of 0.25
for the attention weights. In consequence, the information flowing from the
first graph attention layer (i.e., edges connecting leaf nodes to nodes at
depth 1) will be diminished by $0.25$ as it reaches the root node (P2).

####  3.3.3 Proposed method

To design the edge attribution function $\Phi$, we start by formally defining
the computation tree alongside flow and attention flow.

######  Definition 3.1 (Computation tree).

The computation tree for an $L$-layer GAT in our study is defined to as a
rooted subtree of height $L$ with the root (target) node. For each node in the
tree at depth $d$, the neighboring nodes are at depth $d+1$ with edges
directed towards node $v$.

According to Definition 3.1, we define the concept of flows.

######  Definition 3.2 (Flow in a computation tree).

Given a computation tree as a rooted subtree of height $L$ with the root
(target) node $v$, we define a flow $\lambda_{i,j,v}^{l}$ as the list of edges
that sequentially appear in a path of length $l$ starting from a given edge
$e_{i,j}$ within the computation tree and ending with some edge
$e_{*,v}$.222As stated in (O1), nodes/edges are not unique in the computation
tree. Nonetheless, we will use the node indices from the original graph and
avoid differentiating them in the computation tree as long as it does not
cause any confusion. We indicate the $k$-th position within the flow as
$\lambda_{i,j,v}^{l}(k)$ for $k\in[1,L]$. We denote the set of all flows in
the computation tree with node $v$ at its root that starts from edge $e_{i,j}$
with length $m\in[1,L]$ as $\Lambda_{v}^{m}(e_{i,j})$.

From Definition 3.2, it follows that $\lambda_{i,j,v}^{l}(1)=e_{i,j}$ and
$\lambda_{i,j,v}^{l}(l)=e_{*,v}$ for all flows in $\Lambda_{v}(e_{i,j})$.

######  Definition 3.3 (Attention flow in a computation tree).

Given a flow $\lambda_{v_{0},v_{1},w}^{m}=[e_{v_{0},v_{1}},\cdots,e_{*,w}]$ of
length $m\leq L$ for an $L$-layer GAT model, we define an attention flow
$\alpha[\lambda_{v_{0},v_{1},w}^{m}]$ as the corresponding attention weights
assigned to each edge by the associated graph attention layers:

$\alpha[\lambda_{v_{0},v_{1},w}^{m}]=[\alpha_{v_{0},v_{1}}^{L-m+1},\cdots,\alpha_{*,w}^{L}].$
(1)

Then, it follows that
$\alpha[\lambda_{v_{0},v_{1},w}^{m}](i)=\alpha_{v_{i-1},v_{i}}^{L-m+i}$.

######  Example 1.

In Figure LABEL:subfig:infectiondatasetcomputationgraph,
$\Lambda_{27}(e_{40,27})$ includes two flows, i.e.,
$\lambda_{40,27,27}^{1}=[e_{40,27}]$ and
$\lambda_{40,27,27}^{2}=[e_{40,27},e_{27,27}]$, along with the corresponding
attention flows $\alpha[\lambda_{40,27,27}^{1}]=[0.25]$ and
$\alpha[\lambda_{40,27,27}^{2}]=[0.9,0.25]$, respectively.

Finally, we are ready to present GAtt.

######  Definition 3.4 (GAtt).

Given a target node $v$, an edge $e_{i,j}$ of interest, the set of flows,
$\Lambda_{v}(e_{i,j})$, and the attention flows for all flows in
$\Lambda_{v}(e_{i,j})$, we define the edge attribution of $e_{i,j}$ in
$L$-layer GATs as

$\phi^{v}_{i,j}=\sum_{m^{\prime}=1}^{L}\sum_{\lambda_{i,j,v}^{m^{\prime}}\in\Lambda_{v}^{m^{\prime}}(e_{i,j})}C(\alpha[\lambda_{i,j,v}^{m^{\prime}}])\alpha[\lambda_{i,j,v}^{m^{\prime}}](1),$
(2)

where $C(\alpha[\lambda_{i,j,v}^{m}])=\prod_{2\leq k\leq
m}\alpha[\lambda_{i,j,v}^{m}](k)$ (or $1$ if $m=1$). Equation (2) can be
interpreted as follows. We first find all occurrences of the target edge
$e_{i,j}$ in the computation tree, and then re-weight its attention score
(i.e., $\alpha[\lambda_{i,j,v}^{m}](1)$) by the product of all attention
weights that appear after $e_{i,j}$ (i.e., $\alpha[\lambda_{i,j,v}^{m}](k)$
for $k\geq 2$) in the flow before the summation over all relevant flows. Next,
let us turn to addressing how our design principles (P1) and (P2) are met.
First, (P1) holds as we add the contributions from each flow rather than
taking the average, as this will offset the total number of occurrences of
$e_{i,j}$ in the computation tree. Next, (P2) is fulfilled by the adjustment
factor from $C(\alpha[\lambda_{i,j,v}^{m}])$, since its value is dependent on
the position of $\lambda_{i,j,v}^{m}(1)$. Essentially,
$C(\alpha[\lambda_{i,j,v}^{m}])$ takes the chain of calculation from an edge
to the target node into account. We provide an insightful example below.

######  Example 2.

Let us recall $\lambda^{2}_{40,27,27}=[e_{40,27},e_{27,27}]$ and its attention
flow $\alpha[\lambda^{2}_{40,27,27}]=[0.9,0.25]$ on node 27 from Example 1. At
face value, the contribution of edge $e_{40,27}$ within the flow
$\lambda^{2}_{40,27,27}$ should be $0.9$. However, this is inappropriate since
the information will eventually get muted significantly by
$\alpha_{27,27}^{2}=0.25$; thus, we need to consider
$C(\alpha[\lambda^{2}_{40,27,27}])$ before calculating the final edge
attribution. From Definition 3.4, the edge attribution $\phi^{27}_{40,27}$
from the attention weights is calculated as:

$\phi^{27}_{40,27}=1\times 0.25+0.25\times 0.9=0.475.$ (3)

###  3.4 Practical Calculation of GAtt

Although GAtt is defined by Equation 2, directly using this to compute the
edge attribution $\phi_{i,j}^{v}$ is not desirable since it involves
constructing the computation tree in the form of a rooted subtree for each
node $v$, as well as computing over all relevant attention flows, resulting in
high redundancy during computation and not being proper for batch computation.
To overcome these challenges, we introduce a matrix-based computation method
that is much preferred in practice.

We first define

$\mathbf{C}_{L}(k)=\begin{cases}{\bf I},&\text{if }k=0,\\\ {\bf A}(L){\bf
A}(L-1)\cdots{\bf A}(L-k+1),&\text{otherwise.}\end{cases}$

Then, we would like to establish the following proposition.

######  Proposition 3.5.

For a given set of attention weights $\mathcal{A}=\\{{\bf A}(l)\\}_{l=1}^{L}$
for an $L$-layer GAT with $L\geq 1$, GAtt in Definition 3.4 is equivalent to

$\phi^{v}_{i,j}=\sum_{m=1}^{L}[\mathbf{C}_{L}(L-m)]_{v,j}[{\bf A}(m)]_{j,i}.$
(4)

Proposition 3.5 signifies that GAtt sums the attention scores, weighted by the
sum of the products of attention weights $[{\bf A}(m)]_{j,i}$ along the paths
from node $j$ to node $v$, over all graph attention layers.

###  3.5 Complexity Analysis

We first analyse the computational complexity of GAtt. According to Equation
4, the bottleneck for calculating $\phi_{i,j}^{v}$ is to acquire
$\prod_{k=m+1}^{L}{\bf A}(k)$. However, this matrix can be pre-computed and
does not require re-calculation after its initial acquirement. Since we only
count the number of multiplications in the summation, the computation
complexity is $O(L)$, which is extremely efficient. Next, according to
Equation 4, the memory complexity requires analyzing $\mathbf{C}_{L}(L-m)$ and
$\mathbf{A}(m)$. For an $L$-layer GAT, while storing all attention weights in
$\mathbf{A}(m)$ requires $O(L|\mathcal{E}|)$, $\mathbf{C}_{L}(L-m)$ requires
at most $O(L||T^{L-1}||_{0})$, where $T$ denotes the adjacency matrix, and
$||\cdot||_{0}$ is the 0-norm. In conclusion, the total memory complexity is
$O(L||T^{L-1}||_{0}+L|\mathcal{E}|)$.

##  4 Can Attention Interpret Att-GNNs?

In this section, we carry out empirical studies to validate the effectiveness
of GAtt, the proposed method interpreting a given GAT model.

###  4.1 Is Attention Faithful to the Model?

Table 1: Experimental results with respect to the faithfulness for GAtt, AvgAtt, and random attribution on four real-world datasets. Results for 2-layer and 3-layer GATs are shown for each case. The best performer is highlighted as bold. Dataset | Measure | Metric |  | 2-layer GAT |  |  | 3-layer GAT |
---|---|---|---|---|---|---|---|---
GAtt | AvgAtt | Random | GAtt | AvgAtt | Random
Cora | $\Delta_{\text{PC}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.8468 | 0.1764 | -0.0056 | 0.8642 | 0.0967 | 0.0045
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.7051 | -0.1826 | 0.0082 | 0.6512 | -0.0537 | -0.0025
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.6516 | -0.1240 | 0.0061 | 0.5679 | -0.0379 | -0.0018
$\Delta_{\text{NE}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.7112 | 0.1526 | -0.0076 | 0.7690 | 0.0859 | 0.0040
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.7948 | -0.2463 | 0.0060 | 0.7616 | -0.0820 | 0.0007
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.7371 | -0.1736 | 0.0044 | 0.6737 | -0.0580 | 0.0005
$\Delta_{\text{P}}$ |  AUROC ($\uparrow$) | 0.9755 | 0.7251 | 0.4389 | 0.9875 | 0.7075 | 0.5235
Citeseer | $\Delta_{\text{PC}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.8516 | 0.3096 | 0.0012 | 0.8711 | 0.2110 | -0.0073
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.7584 | -0.0106 | 0.0041 | 0.6456 | -0.0130 | 0.0021
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.8321 | -0.0187 | 0.0057 | 0.7318 | -0.0191 | 0.0031
$\Delta_{\text{NE}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.7653 | 0.2780 | 0.0021 | 0.8291 | 0.2006 | -0.0058
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.8469 | -0.0312 | -0.0003 | 0.7235 | -0.0263 | 0.0023
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.9206 | -0.0517 | -0.0004 | 0.8204 | -0.0376 | 0.0032
$\Delta_{\text{P}}$ |  AUROC ($\uparrow$) | 0.9846 | 0.9213 | 0.3695 | 0.9920 | 0.8979 | 0.4039
Pubmed | $\Delta_{\text{PC}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.8812 | 0.1648 | -0.0064 | 0.8489 | 0.0592 | 0.0009
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.6268 | -0.0797 | 0.0002 | 0.5349 | -0.0964 | -0.0003
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.6746 | -0.1097 | -0.0003 | 0.5946 | -0.1348 | -0.0004
$\Delta_{\text{NE}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.8201 | 0.1477 | -0.0068 | 0.8612 | 0.0600 | 0.0015
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.7031 | -0.0823 | 0.0025 | 0.5378 | -0.1138 | -0.0004
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.7568 | -0.1133 | 0.0033 | 0.6187 | -0.1628 | -0.0006
$\Delta_{\text{P}}$ |  AUROC ($\uparrow$) | 0.9915 | 0.8834 | 0.3974 | 0.9993 | 0.8932 | 0.5172
OGB-Arxiv | $\Delta_{\text{PC}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.7790 | 0.0794 | 0.0007 | 0.7721 | 0.0465 | -0.0004
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.2047 | 0.0128 | 0.0009 | 0.1327 | -0.0158 | -0.0041
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.2590 | 0.0187 | 0.0013 | 0.1778 | -0.0246 | -0.0061
$\Delta_{\text{NE}}$ |  $\rho_{\text{Pearson}}$ ($\uparrow$) | 0.8287 | 0.0804 | 0.0016 | 0.8282 | 0.0478 | -0.0017
$\tau_{\text{Kendall}}$ ($\uparrow$) | 0.2619 | 0.0053 | -0.0010 | 0.1557 | -0.0086 | -0.0038
$\rho_{\text{Spearman}}$ ($\uparrow$) | 0.3275 | -0.0066 | -0.0015 | 0.2106 | -0.0142 | -0.0056
$\Delta_{\text{P}}$ |  AUROC ($\uparrow$) | 0.9908 | 0.8470 | 0.4962 | 0.9985 | 0.8331 | 0.5004

We focus primarily on one of the most important properties in evaluating an
explanation method: faithfulness, which indicates how closely it truly
reflects the underlying model’s inner workings (Jacovi & Goldberg, 2020;
Chrysostomou & Aletras, 2021; Liu et al., 2022; Li et al., 2022). Measuring
the faithfulness involves 1) manipulating the input (such as masking a part of
the input) according to the attribution suggested by the explanation method
and 2) observing the change of the model’s response. We specify our experiment
settings below.

####  4.1.1 Attention erasure

To quantitatively analyze the importance of edge $e_{i,j}$, we mask the
attention coefficients (i.e.,, attention weights before softmax) corresponding
to edge $e_{i,j}$ with zeros in the computation tree, which is compared with
the original response of the GAT model (Tomsett et al., 2020). Due to the fact
that removing $e_{i,j}$ instead of erasing attention weights causes the whole
attention distribution to change, we refrain from edge masking in our study.
Moreover, we do not mask the attention weights after softmax, which never
happens in a normal feed-forward procedure of GATs and makes the attention
distribution not properly normalized. In other words, we only remove attention
coefficients from one edge at a time to prevent the model from out-of-
distribution computation, a common pitfall for perturbation-based approaches
(Hooker et al., 2019; Hase et al., 2021).

####  4.1.2 Datasets

In this experiment, we use four citation datasets, Cora, Citeseer, Pubmed
(Yang et al., 2016), and a large-scale dataset, OGB-Arxiv (Hu et al., 2020),
for node classification. They are real-world benchmark datasets where nodes
correspond to papers and edges represent citations between papers. The nodes
are labeled according to their topics.

####  4.1.3 Baseline methods

We compare the proposed GAtt against another attention-based explanation
method (Ying et al., 2019; Luo et al., 2020; Sánchez-Lengeling et al., 2020),
named as AvgAtt, which attributes each edge as the average of the attention
weights over different layers and attention heads. As another baseline, we
additionally include random attribution, called ‘Random’, by randomly
assigning scores between $[0,1]$ to each edge.

####  4.1.4 Measurement of GAT

Denoting the output probability vector of the GAT for node $v$ as ${\bf
p}_{v}$ and the output probability vector after the attention erasure for
$e_{i,j}$ as ${\bf p}_{v\backslash e_{i,j}}$, we measure the model’s behavior
from three points of view: 1) decline in prediction confidence
$\Delta_{\text{PC}}$ (Guo et al., 2017) defined as the decrease of the
probability for the predicted label (i.e., $\Delta_{\text{PC}}={\bf
p}_{v}[i]-{\bf p}_{v\backslash e_{i,j}}[i]$, where
$i=\operatorname*{arg\,max}_{i}{\bf p}_{v}[i]$), 2) change in negative entropy
$\Delta_{\text{NE}}$ (Moon et al., 2020) defined as the increase of
‘smoothness’ of the probability vector (i.e., $\Delta_{\text{NE}}=-\sum{\bf
p}_{v\backslash e_{i,j}}\log{\bf p}_{v\backslash e_{i,j}}+\sum{\bf
p}_{v}\log{\bf p}_{v}$), which also reflects the model’s confidence, and 3)
change in prediction $\Delta_{\text{P}}$, which measures the difference
between $\operatorname*{arg\,max}_{i}{\bf p}_{v}[i]$ and
$\operatorname*{arg\,max}_{i}{\bf p}_{v\backslash e_{i,j}}[i]$, where ${\bf
p}_{v}[i]$ is the $i$-th entry of ${\bf p}_{v}$.

####  4.1.5 Quantitative analysis with erasure

We investigate the relationship between the attribution scores acquired from
attention erasure and edge attribution calculation methods. In each dataset,
we randomly select 100 nodes as target nodes $v$ and calculate GAtt values for
all edges $(i,j)$ that affect the target node (i.e., $\phi^{v}_{i,j}$). We
also perform attention erasure for the same edges $(i,j)$ and measure
$\Delta_{\text{PC}}$, $\Delta_{\text{NE}}$, and $\Delta_{\text{P}}$ to observe
the correlation between GAtt values. Specifically, for $\Delta_{\text{PC}}$
and $\Delta_{\text{NE}}$, we adopt three metrics: Pearson’s rho
$\rho_{\text{Pearson}}$, Kendall’s tau $\tau_{\text{Kendall}}$, and Spearman’s
rho $\rho_{\text{Spearman}}$. For $\Delta_{\text{P}}$, we use the area under
receiver operating characteristic (AUROC), basically assessing the quality of
an attribution score as a predictor of whether the prediction of the target
node will change after attention erasure.

Table 1 summarizes the experimental results with respect to the faithfulness
on the four real-world citation datasets, using 2-layer and 3-layer trained
GATs with a single attention head for each dataset. The results strongly
indicate using GAtt reflects the behavior of the GAT model much better than
AvgAtt, resulting in a drastic increase in explainability performance compared
to AvgAtt and random attribution. Although AvgAtt shows modest performance in
$\Delta_{\text{P}}$, it performs poorly in terms of changes in confidence
(i.e., $\Delta_{\text{PC}}$ and $\Delta_{\text{NE}}$), sometimes performing
worse than random attribution. This is because AvgAtt does not account for the
proximity effect and contribution adjustment, and rather naïvely averages the
attention weights over different layers with no context of the computation
tree.

Table 2: Experimental results with respect to the faithfulness by increasing the number of attention heads from 1 to 8 for 2 and 3-layer GAT models. The table shows the case for the Cora dataset, measuring $\Delta_{\text{PC}}$ with $\rho_{\text{Spearman}}$.  Model | Method | Number of heads
---|---|---
1 | 2 | 4 | 8
2-layer | GAtt | 0.8477 | 0.8625 | 0.8468 | 0.8496
AvgAtt | 0.1768 | 0.1807 | 0.1697 | 0.1728
Random | -0.0079 | 0.0031 | 0.0021 | 0.0073
3-layer | GAtt | 0.8624 | 0.8857 | 0.8674 | 0.7048
AvgAtt | 0.0966 | 0.0965 | 0.0994 | 0.0857
Random | 0.0092 | -0.0011 | 0.0001 | -0.0033

Additionally, we consider the effect of multi-head attention, which allows
each attention head to potentially capture different patterns in GATs. To this
end, we carry out an experiment while increasing the number of attention heads
from 1 to 8. Table 2 summarizes the performance of edge attribution
calculation methods according to different attention head configurations. We
find that the trend in performance is largely consistent with our findings in
Table 1 even with an increased number of attention heads. However, we find no
distinguishable pattern with respect to the number of heads in both 2-layer
and 3-layer GATs for the Cora dataset; nonetheless, this is not necessarily
the case for other two datasets.

###  4.2 Does Attention Reveal Accurate Explanations?

We evaluate the edge attributions of the GAT model in comparison with ground
truth explanations. Since only the synthetic datasets are equipped with proper
ground truth explanations, we only use these datasets during evaluation.

####  4.2.1 Datasets

We use the BA-shapes and Infection synthetic benchmark datasets, which also
includes ground truth explanations. BA-shapes (Ying et al., 2019) attaches 80
house-shaped motifs to a base graph made from the Barabási-Albert model with
300 nodes, where the edges included in the motif is set as the ground-truth
explanation. Infection benchmark (Faber et al., 2021) generates a backbone
graph from the Erdös-Rényi model. Then, a small portion of the nodes are
assigned as ‘infected’, and the ground-truth explanation is the path from an
infected node to the target node. We expect that edge attributions should
highlight such ground truth explanations for the underlying GAT with high
performance.

####  4.2.2 Baseline methods

In this experiment, we mainly compare the performance among attention-based
edge attribution calculation methods (i.e., GAtt and AvgAtt) including random
attribution. Additionally, we consider seven popular post-hoc explanation
methods: Saliency (Simonyan et al., 2014), Guided Backpropagation
(Springenberg et al., 2015), Integrated Gradient (Sundararajan et al., 2017),
GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), GraphMask
(Schlichtkrull et al., 2021), and FastDnX (Pereira et al., 2023). We emphasize
that post-hoc explanation methods are treated as a complementary category of
inherent explanations, thus belonging to a different category (Du et al.,
2020). However, we include them for a more comprehensive comparison.

####  4.2.3 Measurement

As in prior studies (Ying et al., 2019; Luo et al., 2020), we treat evaluation
as binary classification of edges, aiming to predict whether each edge belongs
to ground truth explanations by using the attribution scores as probability
values. In this context, we adopt the AUROC as our metric.

Table 3: Experimental results on the explanation accuracy for the synthetic datasets using 3-layer GATs, measured in terms of the AUROC. The results for directly using attention weights as explanation are colored in gray. The best and runner up performances are marked as bold and underline, respectively. Expl. type | Method | BA-Shapes | Infection
---|---|---|---
| GAtt | 0.9591 | 0.9976
Attention | AvgAtt | 0.7977 | 0.8786
- | Random | 0.4975 | 0.4811
Post-hoc | Saliency | 0.9563 | 0.8237
Guided Backprop | 0.6231 | 0.8949
Integrated Gradient | 0.6231 | 0.9472
GNNExplainer | 0.8912 | 0.9272
PGExplainer | 0.8289 | 0.7173
GraphMask | 0.5316 | 0.6859
FastDnX | 0.9917 | 0.6574

####  4.2.4 Experimental results

Table 3 summarizes the results on the explanation accuracy for two synthetic
datasets with ground truth explanations. For both datasets, we observe that
GAtt is superior to AvgAtt. Even compared to the representative post-hoc
explanation methods, GAtt performs surprisingly better. This indicates that
the attention weights can inherently capture the GAT model’s behavior without
external post-hoc explanation methods.

###  4.3 Ablation Study

The definition of GAtt by Definition 3.4 satisfies the two design principles
(i.e., proximity effect (P1) and contribution adjustment (P2)). We now perform
an ablation study to validate the effectiveness of the contribution adjustment
principle. To this end, we devise two variants GAtt${}_{\textsc{sim}}$ and
GAtt${}_{\textsc{avg}}$, simply adding all attention weights uniformly
corresponding to the target edge in the computation tree and replacing the
weighted summation in Equation 2 with averaging to remove the effects of the
proximity effect, respectively. Specifically, GAtt${}_{\textsc{sim}}$ and
GAtt${}_{\textsc{avg}}$ are defined as,

$\displaystyle\phi^{v}_{i,j}$
$\displaystyle=\sum_{m^{\prime}=1}^{L}\sum_{\lambda_{i,j,v}^{m^{\prime}}\in\Lambda^{m^{\prime}}_{v}(e_{i,j})}\alpha[\lambda_{i,j,v}^{m^{\prime}}](1),$
(5) $\displaystyle\phi^{v}_{i,j}$
$\displaystyle=\dfrac{1}{|\Lambda_{v}(e_{i,j})|}\sum_{m^{\prime}=1}^{L}\sum_{\lambda_{i,j,v}^{m^{\prime}}\in\Lambda^{m^{\prime}}_{v}(e_{i,j})}C(\alpha[\lambda_{i,j,v}^{m^{\prime}}])\alpha[\lambda_{i,j,v}^{m^{\prime}}](1),$
(6)

respectively. The properties of different edge attribution calculation methods
are summarized in Table 4.

Table 4: Properties of different edge attribution calculation methods.

Method | GAtt | GAtt${}_{\textsc{sim}}$ | GAtt${}_{\textsc{avg}}$ | AvgAtt
---|---|---|---|---
(P1) Proximity effect | ✓ | ✓ | ✗ | ✗
(P2) Contribution adjustment | ✓ | ✗ | ✓ | ✗

Table 5: Properties of different edge attribution calculation methods.

Dataset | Model | GAtt | GAtt${}_{\textsc{sim}}$ | GAtt${}_{\textsc{avg}}$ | AvgAtt
---|---|---|---|---|---
Cora | 2-layer | 0.7793 | 0.7565 | 0.7650 | 0.0461
3-layer | 0.7598 | 0.6603 | 0.6940 | 0.0954
Citeseer | 2-layer | 0.8514 | 0.8295 | 0.7317 | 0.1995
3-layer | 0.8019 | 0.7210 | 0.7890 | 0.1734
Pubmed | 2-layer | 0.7792 | 0.7413 | 0.6796 | 0.1158
3-layer | 0.7136 | 0.5727 | 0.7051 | 0.0721

Figure 3: Case study on the BA-Shapes and Infection datasets for a 2-layer
GAT.

Evaluation of faithfulness. We compare the performance among GAtt,
GAtt${}_{\textsc{sim}}$, GAtt${}_{\textsc{avg}}$, and AvgAtt by running
experiments with respect to the faithfulness on the Cora, Citeseer, and Pubmed
datasets. Table 5 summarizes the performance of each method by averaging over
all 7 measures and metrics (i.e., 3 correlation metrics for
$\Delta_{\text{PC}}$ and $\Delta_{\text{NE}}$ each and AUROC for
$\Delta_{\text{p}}$). Comparing GAtt and its variants, we observe that GAtt
consistently outperforms both GAtt${}_{\textsc{sim}}$ and
GAtt${}_{\textsc{avg}}$ for all cases. In particular, the performance
degradation of the variants is generally more severe for 3-layer GATs. This is
because the effect of the contribution adjustment (P2) and the cardinality of
$\Lambda_{v}(e_{i,j})$ are more significant in 3-layer GATs since the length
of each attention flow are longer as well as the number of flows to consider
is much higher compared to 2-layer GATs.

Case study. We conduct case studies on the BA-Shapes and Infection datasets
for 2-layer GATs. In Figure 3, each of two cases shows the target node (orange
node) and the ground truth explanations (blue edges). We aim to observe how
much the attribution scores from different methods focus on the ground truth
explanation edges. Indeed, for both datasets, GAtt focuses primarily on the
ground truth explanations, while the attribution scores from
GAtt${}_{\textsc{sim}}$ and AvgAtt tend to be more spread throughout the whole
graph. This indicates that the attention weights in GATs indeed recognize the
ground truth explanations. For GAtt${}_{\textsc{avg}}$, it is likely that the
attributions are not much different from GAtt in BA-Shapes. However,
GAtt${}_{\textsc{avg}}$ in the Infection dataset attributes its attribution
scores to a single self-loop edge that does not belong to the ground truth
explanations, failing to provide adequate explanations. Interestingly, on BA-
Shapes, GAtt tends to strongly emphasize local edges even within the house-
shaped motifs, which coincides with the pitfall addressed in (Faber et al.,
2021).

##  5 Conclusion and Future Work

In this study, we have investigated the largely under-explored problem of
interpreting Att-GNNs by utilizing the attention weights. Although Att-GNNs
were not regarded as a candidate as an inherently explainable model, our
empirical evaluations have demonstrated the affirmative results when our
proposed method, GAtt, is used to calculate edge attribution scores. Even
though GAtt is generally applicable, this work does not include a systematic
analysis on how different designs of attention weights will interact with
GAtt, which we leave as future work.

## References

  * Abnar & Zuidema (2020) Abnar, S. and Zuidema, W. H.  Quantifying attention flow in transformers.  In _ACL_ , pp. 4190–4197, Online, Jul. 2020.
  * Bahdanau et al. (2015) Bahdanau, D., Cho, K., and Bengio, Y.  Neural machine translation by jointly learning to align and translate.  In _ICLR_ , San Diego, CA, May 2015.
  * Bibal et al. (2022) Bibal, A., Cardon, R., Alfter, D., Wilkens, R., Wang, X., François, T., and Watrin, P.  Is attention explanation? an introduction to the debate.  In _ACL_ , pp. 3889–3900, Dublin, Ireland, May 2022.
  * Bronstein et al. (2021) Bronstein, M. M., Bruna, J., Cohen, T., and Velickovic, P.  Geometric deep learning: Grids, groups, graphs, geodesics, and gauges, 2021.
  * Caron et al. (2021) Caron, M., Touvron, H., Misra, I., Jégou, H., Mairal, J., Bojanowski, P., and Joulin, A.  Emerging properties in self-supervised vision transformers.  In _ICCV_ , pp. 9630–9640, Montreal, Canada, Oct. 2021.
  * Chefer et al. (2021a) Chefer, H., Gur, S., and Wolf, L.  Transformer interpretability beyond attention visualization.  In _CVPR_ , pp. 782–791, virtual, Jun. 2021a.
  * Chefer et al. (2021b) Chefer, H., Gur, S., and Wolf, L.  Generic attention-model explainability for interpreting bi-modal and encoder-decoder transformers.  In _ICCV_ , pp. 387–396, Montreal, Canada, Oct. 2021b.
  * Chen et al. (2023) Chen, J., Gao, K., Li, G., and He, K.  NAGphormer: A tokenized graph transformer for node classification in large graphs.  In _ICLR_ , May 2023.
  * Chrysostomou & Aletras (2021) Chrysostomou, G. and Aletras, N.  Improving the faithfulness of attention-based explanations with task-specific information for text classification.  In _ACL/IJCNLP_ , pp. 477–488, Virtual Event, Aug. 2021.
  * Dosovitskiy et al. (2021) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _ICLR_ , Virtual event, May 2021.
  * Du et al. (2020) Du, M., Liu, N., and Hu, X.  Techniques for interpretable machine learning.  _Commun. ACM_ , 63(1):68–77, 2020.
  * Faber et al. (2021) Faber, L., Moghaddam, A. K., and Wattenhofer, R.  When comparing to ground truth is wrong: On evaluating GNN explanation methods.  In _KDD_ , pp. 332–341, Virtual Event, Aug. 2021.
  * Ghaeini et al. (2018) Ghaeini, R., Fern, X. Z., and Tadepalli, P.  Interpreting recurrent and attention-based neural models: a case study on natural language inference.  In _EMNLP_ , pp. 4952–4957, Brussels, Belgium, Nov. 2018.
  * Gilmer et al. (2017) Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E.  Neural message passing for quantum chemistry.  In _ICML_ , pp. 1263–1272, Sydney, NSW, Aug. 2017.
  * Guo et al. (2017) Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q.  On calibration of modern neural networks.  In _ICML_ , pp. 1321–1330, Sydney, Australia, Aug. 2017.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  In _NeurIPS_ , pp. 1024–1034, Long Beach, CA, Dec. 2017.
  * Hao et al. (2021) Hao, Y., Dong, L., Wei, F., and Xu, K.  Self-attention attribution: Interpreting information interactions inside transformer.  In _AAAI_ , pp. 12963–12971, Virtual event, Feb. 2021.
  * Hase et al. (2021) Hase, P., Xie, H., and Bansal, M.  The out-of-distribution problem in explainability and search methods for feature importance explanations.  In _NeurIPS_ , pp. 3650–3666, Virtual event, Dec. 2021.
  * Hooker et al. (2019) Hooker, S., Erhan, D., Kindermans, P., and Kim, B.  A benchmark for interpretability methods in deep neural networks.  In _NeurIPS_ , pp. 9734–9745, Vancouver, Canada, Dec. 2019.
  * Hu et al. (2020) Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., and Leskovec, J.  Open graph benchmark: Datasets for machine learning on graphs.  In _NeurIPS_ , Virtual event, Dec. 2020.
  * Jacovi & Goldberg (2020) Jacovi, A. and Goldberg, Y.  Towards faithfully interpretable NLP systems: How should we define and evaluate faithfulness?  In _ACL_ , pp. 4198–4205, Online, Jul. 2020.
  * Jain & Wallace (2019) Jain, S. and Wallace, B. C.  Attention is not explanation.  In _NAACL-HLT_ , pp. 3543–3556, Minneapolis, MN, Jun. 2019.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _ICLR_ , Toulon, France, Apr. 2017.
  * Kreuzer et al. (2021) Kreuzer, D., Beaini, D., Hamilton, W. L., Létourneau, V., and Tossou, P.  Rethinking graph transformers with spectral attention.  In _NeurIPS_ , pp. 21618–21629, Virtual event, Dec. 2021.
  * Lee et al. (2017) Lee, J., Shin, J., and Kim, J.  Interactive visualization and manipulation of attention-based neural machine translation.  In _EMNLP_ , pp. 121–126, Copenhagen, Denmark, Sept. 2017.
  * Li et al. (2022) Li, P., Yang, Y., Pagnucco, M., and Song, Y.  Explainability in graph neural networks: An experimental survey, 2022.
  * Liu et al. (2022) Liu, Y., Li, H., Guo, Y., Kong, C., Li, J., and Wang, S.  Rethinking attention-model explainability through faithfulness violation test.  In _ICML_ , pp. 13807–13824, Baltimore, MD, Jul. 2022.
  * Luo et al. (2020) Luo, D., Cheng, W., Xu, D., Yu, W., Zong, B., Chen, H., and Zhang, X.  Parameterized explainer for graph neural network.  In _NeurIPS_ , Virtual event, Dec. 2020.
  * Miao et al. (2022) Miao, S., Liu, M., and Li, P.  Interpretable and generalizable graph learning via stochastic attention mechanism.  In _ICML_ , Baltimore, MD, Jul. 2022.
  * Moon et al. (2020) Moon, J., Kim, J., Shin, Y., and Hwang, S.  Confidence-aware learning for deep neural networks.  In _ICML_ , pp. 7034–7044, Vienna, Austria, Jul. 2020.
  * Pereira et al. (2023) Pereira, T. A., Nascimento, E., Resck, L. E., Mesquita, D., and Souza, A. H.  Distill n’ explain: explaining graph neural networks using simple surrogates.  In _AISTATS_ , Palau de Congressos, Spain, Apr. 2023.
  * Sánchez-Lengeling et al. (2020) Sánchez-Lengeling, B., Wei, J. N., Lee, B. K., Reif, E., Wang, P., Qian, W. W., McCloskey, K., Colwell, L. J., and Wiltschko, A. B.  Evaluating attribution for graph neural networks.  In _NeurIPS_ , virtual, Dec. 2020.
  * Sato et al. (2021) Sato, R., Yamada, M., and Kashima, H.  Random features strengthen graph neural networks.  In _SDM_ , pp. 333–341, Virtual event, Apr.–May 2021.
  * Schlichtkrull et al. (2021) Schlichtkrull, M. S., Cao, N. D., and Titov, I.  Interpreting graph neural networks for NLP with differentiable edge masking.  In _ICLR_ , Virtual Event, May 2021.
  * Simonyan et al. (2014) Simonyan, K., Vedaldi, A., and Zisserman, A.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  In _ICLR_ , Banff, Canada, Apr. 2014.
  * Springenberg et al. (2015) Springenberg, J. T., Dosovitskiy, A., Brox, T., and Riedmiller, M. A.  Striving for simplicity: The all convolutional net.  In _ICLR_ , San Diego, CA, May 2015.
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., and Yan, Q.  Axiomatic attribution for deep networks.  In _ICML_ , Sydney, Australia, Aug. 2017.
  * Tomsett et al. (2020) Tomsett, R., Harborne, D., Chakraborty, S., Gurram, P., and Preece2, A.  Sanity checks for saliency metrics.  In _AAAI_ , New York, NY, Feb. 2020.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I.  Attention is all you need.  In _NIPS_ , pp. 5998–6008, Long Beach, CA, Dec. 2017.
  * Velickovic et al. (2018) Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y.  Graph attention networks.  In _ICLR_ , Vancouver, Canada, Apr. 2018.
  * Vig (2019) Vig, J.  A multiscale visualization of attention in the transformer model.  In _ACL_ , pp. 37–42, Florence, Italy, Jul.–Aug. 2019.
  * Wiegreffe & Pinter (2019) Wiegreffe, S. and Pinter, Y.  Attention is not not explanation.  In _EMNLP-IJCNLP_ , pp. 11–20, Hong Kong, China, Nov. 2019.
  * Wu et al. (2021) Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., and Yu, P. S.  A comprehensive survey on graph neural networks.  _IEEE Trans. Neural Networks Learn. Syst._ , 32(1):4–24, 2021.
  * Xu et al. (2015) Xu, K., Ba, J., Kiros, R., Cho, K., Courville, A. C., Salakhutdinov, R., Zemel, R. S., and Bengio, Y.  Show, attend and tell: Neural image caption generation with visual attention.  In _ICML_ , pp. 2048–2057, Lille, France, Jul. 2015.
  * Xu et al. (2019) Xu, K., Hu, W., Leskovec, J., and Jegelka, S.  How powerful are graph neural networks?  In _ICLR_ , New Orleans, LA, May 2019.
  * Yang et al. (2016) Yang, Z., Cohen, W. W., and Salakhutdinov, R.  Revisiting semi-supervised learning with graph embeddings.  In _ICML_ , New York City, NY, Jun. 2016.
  * Ying et al. (2021) Ying, C., Cai, T., Luo, S., Zheng, S., Ke, G., He, D., Shen, Y., and Liu, T.  Do transformers really perform badly for graph representation?  In _NeurIPS_ , pp. 28877–28888, virtual, Dec. 2021.
  * Ying et al. (2019) Ying, Z., Bourgeois, D., You, J., Zitnik, M., and Leskovec, J.  GNNExplainer: Generating explanations for graph neural networks.  In _NeurIPS_ , pp. 9240–9251, Vancouver, Canada, Dec. 2019.
  * Yuan et al. (2023) Yuan, H., Yu, H., Gui, S., and Ji, S.  Explainability in graph neural networks: A taxonomic survey.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 45(5):5782–5799, 2023.