# Global Counterfactual Explainer for Graph Neural Networks

Zexi Huang  University of CaliforniaSanta BarbaraCAUSA<EMAIL_ADDRESS>,  Mert Kosan  University of CaliforniaSanta BarbaraCAUSA
<EMAIL_ADDRESS>,  Sourav Medya  University of IllinoisChicagoILUSA
<EMAIL_ADDRESS>,  Sayan Ranu  Indian Institute of TechnologyDelhiIndia
<EMAIL_ADDRESS>and  Ambuj Singh  University of CaliforniaSanta
BarbaraCAUSA<EMAIL_ADDRESS>
(2023)

###### Abstract.

Graph neural networks (GNNs) find applications in various domains such as
computational biology, natural language processing, and computer security.
Owing to their popularity, there is an increasing need to explain GNN
predictions since GNNs are black-box machine learning models. One way to
address this is counterfactual reasoning where the objective is to change the
GNN prediction by minimal changes in the input graph. Existing methods for
counterfactual explanation of GNNs are limited to instance-specific _local_
reasoning. This approach has two major limitations of not being able to offer
global recourse policies and overloading human cognitive ability with too much
information. In this work, we study the global explainability of GNNs through
global counterfactual reasoning. Specifically, we want to find a small set of
representative counterfactual graphs that explains all input graphs. Towards
this goal, we propose GCFExplainer, a novel algorithm powered by vertex-
reinforced random walks on an edit map of graphs with a _greedy summary_.
Extensive experiments on real graph datasets show that the global explanation
from GCFExplainer provides important high-level insights of the model behavior
and achieves a 46.9% gain in recourse coverage and a 9.5% reduction in
recourse cost compared to the state-of-the-art local counterfactual
explainers.

Counterfactual explanation; Graph neural networks

††journalyear: 2023††copyright: rightsretained††conference: Proceedings of the
Sixteenth ACM International Conference on Web Search and Data Mining; February
27–March 3, 2023; Singapore††booktitle: Proceedings of the Sixteenth ACM
International Conference on Web Search and Data Mining (WSDM ’23), February
27–March 3, 2023, Singapore††doi: 10.1145/nnnnnnn.nnnnnnn††isbn:
978-1-4503-9132-0/22/02

##  1\. Introduction

Graph Neural Networks (GNNs) (Kipf and Welling, 2017; Hamilton et al., 2017;
Veličković et al., 2018; Klicpera et al., 2018; Wang et al., 2021, 2022) are
being used in many domains such as drug discovery (Jiang et al., 2020), chip
design (Mirhoseini et al., 2020), recommendation engines (Ying et al., 2018),
and rare event prediction (Kosan et al., 2021). Taking the graph(s) as input,
GNNs are trained to perform various downstream tasks that form the core of
various real-world applications. For example, graph classification has been
applied to predict whether a drug would exhibit the desired chemical activity
(Jiang et al., 2020). Similarly, node prediction is used to predict the
functionality of proteins in protein-protein interaction networks (Borgwardt
et al., 2005) and categorize users into roles on social networks (Yang et al.,
2011).

Despite the impressive success of GNNs on predictive tasks, GNNs are black-box
machine learning models. It is non-trivial to explain or reason why a
particular prediction is made by a GNN. Explainability of a prediction model
is important to understand its shortcomings and identify areas for
improvement. In addition, the ability to explain a model is critical towards
making it trustworthy. Owing to this limitation of GNNs, there has been
significant efforts in recent times towards explanation approaches.

(a) Formaldehyde

(b) Formic acid

Figure 1. Formaldehyde (a) is classified by a GNN to be an undesired mutagenic
molecule with its important subgraph found by factual reasoning highlighted in
red. Formic acid (b) is its non-mutagenic counterfactual example obtained by
removing one edge and adding one node and two edges.

Existing work on explaining GNN predictions can be categorized mainly in two
directions: 1) factual reasoning (Ying et al., 2019; Luo et al., 2020; Vu and
Thai, 2020; Yuan et al., 2020), and 2) counterfactual reasoning (Lucic et al.,
2022; Bajaj et al., 2021; Abrate and Bonchi, 2021; Tan et al., 2022).
Generally speaking, the methods in the first category aim to find an important
subgraph that correlates most with the underlying GNN prediction. In contrast,
the methods with counterfactual reasoning attempt to identify the smallest
amount of perturbation on the input graph that changes the GNN’s prediction,
for example, removal/addition of edges or nodes.

Compared to factual reasoning, counterfactual explainers have the additional
advantage of providing the means for recourse (Voigt and Von dem Bussche,
2017). For example, in the applications of drug discovery (Jiang et al., 2020;
Xiong et al., 2021), mutagenicity is an adverse property of a molecule that
hampers its potential to become a marketable drug (Kazius et al., 2005). In
Figure 1, formaldehyde is classified by a GNN to be mutagenic. Factual
explainers can attribute the subgraph containing the carbon-hydrogen bond to
the cause of mutagenicity, while counterfactual explainers provide an
effective way (i.e., a recourse) to turn formaldehyde into formic acid, which
is non-mutagenic, by replacing a hydrogen atom with a hydroxyl.

In this work, we focus on counterfactual explanations. Our work is based on
the observation that existing counterfactual explainers (Ying et al., 2019;
Luo et al., 2020; Vu and Thai, 2020; Yuan et al., 2020) for graphs take a
local perspective, generating counterfactual examples for individual input
graphs. However, this approach has two key limitations:

  * •

Lack of global insights:  It is desirable to offer insights that generalize
across a multitude of data graphs. For example, instead of providing formic
acid as a counterfactual example to formaldehyde, we can summarize global
recourse rules such as “ _Given any molecule with a carbonyl group (carbon-
oxygen double bond), it needs a hydroxy to be non-mutagenic_ ”. This focus on
global counterfactual explanation promises to provide higher-level insights
that are complementary to those obtained from local counterfactual
explanations.

  * •

Information overload: The primary motivation behind counterfactual analysis is
to provide human-intelligible explanations. With this objective, consider
real-world graph datasets that routinely contain thousands to millions of
graphs. Owing to instance-specific counterfactual explanations, the number of
counterfactual graphs grows linearly with the graph dataset size.
Consequently, the sheer volume of counterfactual graphs overloads human
cognitive ability to process this information. Hence, the initial motivation
of providing human-intelligible insights is lost if one does not obtain a
holistic view of the counterfactual graphs.

Contributions: In this paper, we study the problem of model-agnostic, _global_
counterfactual explanations of GNNs for graph classification. More
specifically, given a graph dataset, our goal is to counterfactually explain
the largest number of input graphs with a small number of counterfactuals. As
we will demonstrate later in our experiments, this formulation naturally
forces us to remove redundancy from instance-specific counterfactual
explanations and hence has higher information density. Algorithmically, the
proposed problem introduces new challenges. We theoretically establish that
the proposed problem is NP-hard. Furthermore, the space of all possible
counterfactual graphs itself is exponential. Our work overcomes these
challenges and makes the following contributions:

  * •

Novel formulation: We formulate the novel problem of global counterfactual
reasoning/explanation of GNNs for graph classification. In contrast to
existing works on counterfactual reasoning that only generate instance-
specific examples, we provide an explanation on the global behavior of the
model.

  * •

Algorithm design: While the problem is NP-hard, we propose GCFExplainer, which
organizes the exponential search space as an edit map. We then perform vertex-
reinforced random walks on it to generate diverse, representative
counterfactual candidates, which are greedily summarized as the global
explanation.

  * •

Experiments: We conduct extensive experiments on real-world datasets to
validate the effectiveness of the proposed method. Results show that
GCFExplainer not only provides important high-level insights on the model
behavior but also outperforms state-of-the-art baselines related to
counterfactual reasoning in various recourse quality metrics.

##  2\. Global Counterfactual Explanations

This section introduces the global counterfactual explanation (GCE) problem
for graph classification. We start with the background on local counterfactual
reasoning. Then, we propose a representation of the global recourse rule that
provides a high-level counterfactual understanding of the classifier behavior.
Finally, we introduce quality measures for recourse rules and formally define
the GCE problem.

###  2.1. Local Counterfactual

Consider a graph $G=(V,E)$, where $V$ and $E$ are the sets of (labelled) nodes
and edges respectively. A (binary) graph classifier (e.g., a GNN) $\phi$
classifies $G$ into either the undesired class ($\phi(G)=0$) or the desired
one ($\phi(G)=1$). An explanation of $\phi$ seeks to answer how these
predictions are made. Those based on factual reasoning analyze what properties
$G$ possesses to be classified in the current class while those based on
counterfactual reasoning find what properties $G$ needs to be assigned to the
opposite class.

Existing counterfactual explanation methods take a local perspective.
Specifically, for each input graph $G$, they find a counterfactual (graph) $C$
that is somewhat similar to $G$ but is assigned to a different class. Without
loss of generality, let $G$ belong to the undesired class, i.e., $\phi(G)=0$,
then the counterfactual $C$ satisfies $\phi(C)=1$. The similarity between $C$
and $G$ is quantified by a predefined distance metric $d$, for example, the
number of added/removed edges (Lucic et al., 2022; Bajaj et al., 2021).

In our work, we consider the graph edit distance (GED) (Sanfeliu and Fu,
1983), a more general distance measure, as the distance function to account
for other types of changes. Specifically, $\operatorname{GED}(G_{1},G_{2})$
counts the minimum number of “edits” to convert $G_{1}$ to $G_{2}$. An “edit”
can be the addition or removal of edges and nodes, or change of node labels
(see Figure 2). Moreover, to account for graphs of different sizes, we
normalize the GED by the sizes of graphs:
$\widehat{\operatorname{GED}}(G_{1},G_{2})=\operatorname{GED}(G_{1},G_{2})/(|V_{1}|+|V_{2}|+|E_{1}|+|E_{2}|)$.
Nonetheless, our method can be applied with other graph distance metrics, such
as those based on graph kernels (e.g., RW (Borgwardt et al., 2006), NSPDG
(Costa and De Grave, 2010), WL (Shervashidze et al., 2011)).

Figure 2. Edits between graphs.

The distance function measures the quality of the counterfactual found by the
explanation model. Ideally, the counterfactual $C$ should be very close to the
input graph $G$ while belonging to a different class. Formally, we define the
counterfactuals that are within a certain distance $\theta$ from the input
graph as close counterfactuals.

######  Definition 0 (Close Counterfactual).

Given the GNN classifier $\phi$, distance parameter $\theta$, and an input
graph $G$ with undesired outcome, i.e., $\phi(G)=0$; a counterfactual graph,
$C$, is a close counterfactual of $G$ when $\phi(C)=1$ and $d(G,C)\leq\theta$.

While the (close) counterfactual $C$ found by existing methods explains the
classifier behavior for the corresponding input graph $G$, it is hard to
generalize to understand the global pattern. Next, we introduce the global
recourse rule that provides a high-level summary of the classifier behavior
across different input graphs.

###  2.2. Global Recourse Representation

The global counterfactual explanation requires a global recourse rule $r$.
Specifically, for any (undesired) input graph $G$ with $\phi(G)=0$, $r$
provides a (close) counterfactual (i.e., a recourse) for $G$: $\phi(r(G))=1$.
While both a recourse rule and a local counterfactual explainer find a
counterfactual given an input graph, their goals are different. The goal of
the local counterfactual explainer is to find the best (closest)
counterfactual possible for each input graph, and therefore, $r$ can be very
complicated, e.g., in the form of an optimization algorithm (Bajaj et al.,
2021; Tan et al., 2022). On the other hand, a recourse rule aims to provide an
explanation of the classifier’s global behavior, which requires a simpler form
that is understandable for domain experts without prior knowledge of deep
learning on graphs.

Existing global recourse rules for classifiers with feature vectors as input
take the form of short decision trees (Rawal and Lakkaraju, 2020). However,
this is hard to be generalized to graph data with rich structure information.
Instead, we propose the representation of a global recourse rule for a graph
classifier to be a collection of counterfactual graphs $\mathbb{C}$ in the
desired class that are diverse and representative enough to capture its global
behavior. This representation does not require any additional knowledge for
domain experts to understand and draw insights from, similar to the local
counterfactual examples. It is also easy to find the local counterfactual for
a given input graph $G$ based on $\mathbb{C}$ by nominating the closest graph
in $\mathbb{C}$: $r(G)=\operatorname*{arg\,min}_{C\in\mathbb{C}}d(G,C)$.

###  2.3. Quantifying Recourse Quality

Given a graph classifier $\phi$ and a set of $n$ input graphs $\mathbb{G}$ in
the undesired class, we want to compare the quality of different recourse
representations $\mathbb{C}$. Similar to the quality metrics introduced for
vector data (Rawal and Lakkaraju, 2020), we aim to account for the following
factors:

  1. (1)

Coverage: Like local counterfactual explainers, we want to ensure that
counterfactuals found for individual input graphs are of high quality.
Specifically, we introduce recourse coverage—the proportion of input graphs
that have close counterfactuals from $\mathbb{C}$ under a given distance
threshold $\theta$:

$\operatorname{\textbf{coverage}}(\mathbb{C})=|\\{G\in\mathbb{G}\mid\min_{C\in\mathbb{C}}\left\\{d(G,C)\right\\}\leq\theta\\}|/|\mathbb{G}|$

  2. (2)

Cost: Another quality metric based on local counterfactual quality is the
recourse cost (i.e., the distance between the input graph and its
counterfactual) across the input graphs:

$\operatorname{\textbf{cost}}(\mathbb{C})=\operatorname*{agg}_{G\in\mathbb{G}}\\{\min_{C\in\mathbb{C}}\\{d(G,C)\\}\\}$

where $\operatorname*{agg}$ is an aggregation function, e.g., mean or median.

  3. (3)

Interpretability: Finally, the recourse rule should be easy (small) enough for
human cognition. We quantify the interpretability as the size of recourse
representation:

$\operatorname{\textbf{size}}(\mathbb{C})=|\mathbb{C}|$

###  2.4. Problem Formulation and Characterization

An ideal recourse representation should maximize the coverage while minimizing
the cost and the size. Formally, we define the global counterfactual
explanation problem as follows:

######  Problem 1 (Global Counterfactual Explanation for Graph Classification
(GCE)).

Given a GNN graph classifier $\phi$ that classifies $n$ input graphs
$\mathbb{G}$ to the undesired class 0 and a budget $k\ll n$, our goal is to
find the best recourse representation $\mathbb{C}$ that maximizes the recourse
coverage with size $k$:

$\max_{\mathbb{C}}\operatorname{\textbf{coverage}}(\mathbb{C})\text{ s.t.
}\operatorname{\textbf{size}}(\mathbb{C})=k$

We note that in our problem formulation only coverage and size are explicitly
accounted for, whereas cost is absent. We make this design choice since cost
and coverage are intrinsically opposing forces. Specifically, if we are
willing to allow a high cost, coverage increases since we allow for higher
individual distances between an input graph and its counterfactual. Therefore,
we take the approach of binding the cost to the distance threshold $\theta$ in
the coverage definition. Nonetheless, an explicit analysis of all these
metrics including cost is performed to quantify recourse quality during our
empirical evaluation in Section 4. Below we discuss the hardness of GCE.

######  Theorem 1 (NP-hardness).

The GCE problem is NP-hard.

###### Proof.

To establish NP-hardness of the proposed problem we reduce it from the
classical Maximum Coverage problem.

######  Definition 0 (Maximum Coverage).

Given a budget $k$ and a collection of subsets
$\mathcal{S}=\\{S_{1},\cdots,S_{m}\\}$ from a universe of items
$U=\\{u_{1},\cdots,u_{n}\\}$, find a subset $S^{\prime}\subseteq\mathcal{S}$
of sets such that $|S^{\prime}|\leq k$ and the number of covered elements is
maximized, i.e., maximize $|\bigcup_{\forall S_{i}\in S^{\prime}}S_{i}|$.

We show that given any instance of a maximum coverage problem
$\langle\mathcal{S},U\rangle$, it can be mapped to the GCE problem.
Specifically, we assume that each element $u_{i}$ in $U$ corresponds to an
input undesired graph $G_{i}$. For $u_{i}$, assume a (star) graph with a
center node with empty label and $n$ nodes with empty labels except the node
at $i$-th position that has label $u_{i}$. For $S_{i}$, assume similar graph
with a center node with blank label and nodes with labels with the elements in
$S_{i}$ in their corresponding positions. The edit operations that are allowed
are either making a node label blank or add a new label in their position and
but not both. So, each $S_{i}$ corresponds to a $C_{i}\in\mathbb{C}$ and
$d(G_{j},C_{i})\leq\theta\\}|$ iff $u_{j}\in S_{i}$ where $\theta$ is assumed
large enough. Note that with this construction, it is easy to see that an
optimal solution for this instance of the GCE problem is the optimal solution
for the corresponding instance of the maximum coverage problem. ∎

Owing to NP-hardness, it is not feasible to identify the optimal solution for
the GCE problem in polynomial time. In the next section, we will introduce
GCFExplainer, an effective and efficient heuristic that solves the GCE
problem.

##  3\. Proposed Method: GCFExplainer

In this section, we propose GCFExplainer, the first global counterfactual
explainer for graph classification. The GCE problem requires us to find a
collection of $k$ counterfactual graphs that maximize the coverage of the
input graphs. Intuitively, we want each individual counterfactual graph to be
a close counterfactual to (i.e., “cover”) as many input graphs as possible.
Additionally, different counterfactual graphs should cover different sets of
input graphs to maximize the overall coverage. These intuitions motivate the
design of our algorithm GCFExplainer, which has three major components:

  1. (1)

Structuring the search space: The search space of counterfactual graphs
consists of all graphs that are in the same domain as the input graphs and
within a distance of $\theta$. In other words, any graph within a distance of
$\theta$ from an input graph may be a potential counterfactual candidate and
therefore needs to be analyzed. The number of potential graphs within $\theta$
increases exponentially with $\theta$ since the space of graph edits is
combinatorial (Liang and Zhao, 2017). GCFExplainer uses an edit map to
organize these graphs as a meta-graph $\mathcal{G}$, where individual nodes
are graphs that are created via a different number of edits from the input
graphs and each edge represents a single edit.

  2. (2)

Vertex-reinforced random walk: To search for good counterfactual candidates,
GCFExplainer leverages vertex-reinforced random walks (VRRW) (Pemantle, 1992)
on the edit map $\mathcal{G}$. VRRW has the nice property of converging to a
set of nodes that are both important (i.e., cover many input graphs) and
diverse (i.e., non-overlapping coverage), which will form a small set of
counterfactual candidates for further processing.

  3. (3)

Iterative computation of the summary: After obtaining good counterfactual
candidates from VRRW, GCFExplainer creates the final set of the counterfactual
graphs (i.e., the summary) as the recourse representation by iteratively
adding the best candidate based on the maximal gain of the coverage given the
already added candidates.

###  3.1. Structuring the Search Space

The search space for counterfactual graphs in GCFExplainer is organized via an
edit map $\mathcal{G}$. The edit map is a meta-graph whose nodes are graphs in
the same domain as the input graphs and edges connect graphs that differ by a
single graph edit. As an example, each graph in Figure 2 represents a node in
the edit map, and the arrows denote edges between graphs (nodes) that are one
edit away. In the edit map, we only include connected graphs since real graphs
of interest are often connected (e.g., molecules, proteins, etc.).

While all potential counterfactual candidates are included as its nodes, the
edit map has an exponential size and it is computationally prohibitive to
fully explore it. However, a key observation is that a counterfactual
candidate can only be a few hops away from some input graph. Otherwise, the
graph distance between the counterfactual and the input graph would be too
large for the counterfactual to cover it. This observation motivates our
exploration of the edit map to be focused on the union of close neighborhoods
of the input graphs (see Section 3.2.3). Additionally, while we cannot compute
the entire edit map, it is easy to chart the close neighborhoods by
iteratively performing all possible edits from the input graphs. Next, we
introduce the vertex-reinforced random walk to efficiently explore the edit
map to find counterfactual candidates.

###  3.2. Vertex-Reinforced Random Walk

Vertex-reinforced random walk (VRRW) (Pemantle, 1992) is a time-variant random
walk. Different from other more widely applied random walk processes such as
the simple random walk and the PageRank (Perozzi et al., 2014; Klicpera et
al., 2018; Huang et al., 2021, 2022), the transition probability $p(u,v)$ of
VRRW from node $u$ to node $v$ depends not only on the edge weight $w(u,v)$
but also the number of previous visits in the walk to the target node $v$,
which we denote using $N(v)$. Specifically,

(1) $p(u,v)\propto w(u,v)N(v)$

GCFExplainer applies VRRW on the edit map and produces $n$ most frequently
visited nodes in the walk as the set of counterfactual candidates
$\mathbb{S}$. Next, we formalize VRRW in our setting and explain how it
surfaces good counterfactual candidates for GCE.

####  3.2.1. Vertex-reinforcement

Our main motivation for using VRRW to explore the edit map instead of other
random walk processes is that VRRW converges to a diverse and representative
set of nodes (Mei et al., 2010; Natarajan and Ranu, 2016) in different regions
of the edit map. In this way, the frequently visited nodes in instances of
VRRW have the potential to be good counterfactual candidates as they would
cover a diverse set of input graphs in the edit map. The reason behind the
diversity of the highly visited nodes is the previous visit count $N(v)$ in
the transition probability. Specifically, nodes with larger visit counts tend
to be visited more often later (“richer gets richer”), and thereby dominating
all other nodes in their neighborhood. This leads to a bunch of highly visited
nodes to “represent” each region of the edit map. We refer the readers to (Mei
et al., 2010) for details on the mathematical basis and the theoretical
correctness of this property. Moreover, as our goal is to find counterfactual
candidates, we only reinforce (i.e., increase the visit counts of) graphs in
the counterfactual class.

####  3.2.2. Importance function

While the vertex-reinforcement mechanism ensures diversity of the highly
visited nodes, we still need to guide the walker to visit graphs that are good
counterfactual candidates. We achieve this by assigning large edge weight
$w(u,v)$ to good counterfactual candidates via an importance function $I(v)$:

(2) $w(u,v)=I(v)$

The importance function $I(v)$ should capture the quality of a graph $v$ as a
counterfactual candidate. It has the following components:

  1. (1)

Counterfactual probability $p_{\phi}(v)$. The graph classifier $\phi$ predicts
a probability for $v$ to be in the counterfactual class ($\phi(v)=1$). By
using it as part of the importance function, the walker is encouraged to visit
regions with rich counterfactual graphs.

  2. (2)

Individual coverage $\operatorname{\textbf{coverage}}(\\{v\\})$. The
individual coverage of a graph $v$ computes the proportion of input graphs
that are close to $v$. This encourages the walker to visit graphs that cover a
large number of input graphs.

  3. (3)

Gain of coverage $\operatorname{\textbf{gain}}(v;\mathbb{S})$. Given a graph
$v$ and the current set of counterfactual candidates $\mathbb{S}$ (i.e., the
$n$ most frequently visited nodes), we can compute the gain between the
current coverage and the coverage after adding $v$ to $\mathbb{S}$:

$\operatorname{\textbf{gain}}(v;\mathbb{S})=\operatorname{\textbf{coverage}}(\mathbb{S}\cup\\{v\\})-\operatorname{\textbf{coverage}}(\mathbb{S})$

This guides the walker to find graphs that complement the current
counterfactual candidates to cover additional input graphs.

The importance function is a combination of these components:

(3)
$I(v)=p_{\phi}(v)(\alpha\operatorname{\textbf{coverage}}(\\{v\\})+(1-\alpha)\operatorname{\textbf{gain}}(v;\mathbb{S}))$

where $\alpha$ is a hyperparameter between 0 and 1. With the above importance
function, the VRRW in GCFExplainer converges to a set of diverse nodes that
have high counterfactual probability and collectively cover a large number of
input graphs.

####  3.2.3. Dynamic teleportation

The last component of VRRW, teleportation, is to help us manage the
exponential search space of the edit map. Since our goal is to find close
counterfactuals to the input graphs, the walker only needs to explore the
nearby regions of the input graphs. Therefore, we start the walk from the
input graphs, and also at each step, let the walker teleport back (i.e.
transit) to a random input graph with probability $\tau$.

To decide which input graph to teleport to, we adopt a dynamic probability
distribution based on the current counterfactual candidate set $\mathbb{S}$.
Specifically, let $g(G)=|\\{v\in\mathbb{S}\mid d(v,G)\leq\theta\text{ and
}\phi(v)=1\\}|$ be the number of close counterfactuals in $\mathbb{S}$
covering an input graph $G$. Then the probability to teleport to $G$ is

(4)
$p_{\tau}(G)=\frac{\exp(-g(G))}{\sum_{G^{\prime}\in\mathbb{G}}\exp(-g(G^{\prime}))}$

This dynamic teleportation favors input graphs that are not well covered by
the current solution set and encourages the walker to explore nearby
counterfactuals to cover them after teleportation.

###  3.3. Iterative Computation of the Summary

We have applied VRRW to generate a good set of $n$ counterfactual candidates
$\mathbb{S}$. In the last step of GCFExplainer, we aim to further refine the
candidate set and create the final recourse representation (i.e., the summary)
with $k$ counterfactual graphs. This summarization problem is also NP-hard and
we propose to build $\mathbb{C}$ in an iterative and greedy manner from
$\mathbb{S}$.

Specifically, we start with an empty solution set $\mathbb{C}_{0}$. Then, for
each iteration $t$, we add the graph $v$ to $\mathbb{C}_{t}$ with the maximal
gain of coverage $\operatorname{\textbf{gain}}(v;\mathbb{C}_{t})$. This is
repeated $k$ times to get the final recourse representation $\mathbb{C}$ with
$k$ graphs. It is easy to show that the summarization problem is submodular
and therefore, our greedy algorithm provides $(1-1/e)$-approximation.

Notice that the greedy algorithm can also be applied to the local
counterfactuals found by existing methods to generate a GCE solution. Here, we
highlight three advantages of GCFExplainer:

  1. (1)

Existing local counterfactual explainers (Lucic et al., 2022; Bajaj et al.,
2021; Abrate and Bonchi, 2021; Tan et al., 2022) are only able to generate
counterfactuals based on one type of graph edits—edge removal, while
GCFExplainer incorporates all types of edits to include a richer set of
counterfactual candidates.

  2. (2)

The set of counterfactual candidates from GCFExplainer is generated with the
GCE objective in mind, while the local counterfactuals from existing methods
are optimized for individual input graphs. Therefore, they may not be good
candidates to capture the global behavior of the classifier.

  3. (3)

It is easy to incorporate domain constraints (e.g., the valence of chemical
bonds) into GCFExplainer by pruning the neighborhood of the edit map, while
existing methods based on optimization require non-trivial efforts to
customize.

We will empirically demonstrate the superiority of GCFExplainer to this two-
stage approach with state-of-the-art local counterfactual explanation methods
in our experiments in Section 4.2.

Pseudocode and complexity: The pseudocode of GCFExplainer is presented in
Algorithm 1. Line 1-16 summarizes the VRRW component of GCFExplainer.
Specifically, Line 3-10 determines the next graph to visit based on VRRW
transition probabilities and dynamic teleportation, and Line 11-16 update the
visit counts and the set of counterfactual candidates. The iterative
computation of the counterfactual summary is described in Line 17-21. The
overall complexity of GCFExplainer is $O(Mhn+kn)$, where $M$ is the number of
iterations for the VRRW, $h$ is the average node degree in the meta-graph, $n$
is the number of input graphs, and $k$ is the size of the global
counterfactual representation. In practice, we store the computed transition
probabilities with a space-saving algorithm (Metwally et al., 2005) to improve
the running time of GCFExplainer.

Algorithm 1 GCFExplainer($\phi$, $\mathbb{G}$)

1: $G\leftarrow$ random input graph from $\mathbb{G}$, $N(G)\leftarrow 1$,
$\mathbb{S}=\\{G\\}$

2: for $i\in 1:M$ do

3: Let $\epsilon\sim Bernoulli(\tau)$

4: if $\epsilon=0$ then

5: for $v\in Neighbors(G)$ do

6: Compute $I(v)$ based on Equation 3

7: Compute $p(G,v)$ based on Equation 1

8: $v\leftarrow$ random neighbor of $G$ based on $p(G,v)$

9: else

10: $v\leftarrow$ random input graph from $\mathbb{G}$ based on Equation 4

11: if $\phi(v)=1$ then

12: if $v\in\mathbb{S}$ then

13: $N(v)\leftarrow N(v)+1$

14: else

15: $\mathbb{S}\leftarrow\mathbb{S}+\\{v\\}$, $N(v)\leftarrow 1$

16: $G\leftarrow v$

17: $\mathbb{S}\leftarrow$ top $n$ frequently visited counterfactuals in
$\mathbb{S}$

18: $\mathbb{C}\leftarrow\emptyset$

19: for $t\in 1:k$ do

20:
$v\leftarrow\operatorname*{arg\,max}_{v\in\mathbb{S}}\operatorname{\textbf{gain}}(v;\mathbb{C})$

21: $\mathbb{C}\leftarrow\mathbb{C}+\\{v\\}$

22: return $\mathbb{C}$

##  4\. Experiments

We provide empirical results for the proposed GCFExplanier along with
baselines on commonly used graph classification datasets. Our code is
available at https://github.com/mertkosan/GCFExplainer.

Table 1. The statistics of the datasets.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
#Graphs | 3978 | 4308 | 1837 | 1113
#Nodes | 118714 | 130719 | 28905 | 43471
#Edges | 128663 | 132707 | 29985 | 81044
#Node Labels | 10 | 10 | 9 | 3
Table 2. Recourse coverage ($\theta=0.1$) and median recourse cost comparison between GCFExplainer and baselines for a 10-graph global explanation. GCFExplainer consistently and significantly outperforms all baselines across different datasets.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
| Coverage | Cost | Coverage | Cost | Coverage | Cost | Coverage | Cost
Ground-Truth | 16.54% | 0.1326 | 28.96% | 0.1275 | 0.41% | 0.2012 | 8.47% | 0.2155
RCExplainer | 15.22% | 0.1370 | 31.99% | 0.1290 | 8.96% | 0.1531 | 8.74% | 0.2283
CFF | 17.61% | 0.1331 | 30.43% | 0.1327 | 3.39% | 0.1669 | 3.83% | 0.2557
GCFExplainer | 27.85% | 0.1281 | 37.08% | 0.1135 | 14.66% | 0.1516 | 10.93% | 0.1856
Figure 3. Coverage and cost performance comparison between GCFExplainer and
baselines based on different counterfactual summary sizes. GCFExplainer
consistently outperforms the baselines across different sizes.

###  4.1. Experimental Settings

####  4.1.1. Datasets

We use four different real-world datasets for graph classification benchmark
with their statistics in Table 1. Specifically, NCI1 (Wale and Karypis, 2006),
Mutagenicity (Riesen and Bunke, 2008; Kazius et al., 2005), and AIDS (Riesen
and Bunke, 2008) are collections of molecules with nodes representing
different atoms and edges representing chemical bonds between them. The
molecules are classified by whether they are anticancer, mutagenic, and active
against HIV, respectively. Proteins (Borgwardt et al., 2005; Dobson and Doig,
2003) is a collection of proteins classified into enzymes and non-enzymes,
with nodes representing secondary structure elements and edges representing
structural proximity. For all datasets, we filter out graphs containing rare
nodes with label frequencies smaller than 50.

####  4.1.2. Graph classifier

We follow (Vu and Thai, 2020) and train a GNN with 3 convolution layers (Kipf
and Welling, 2017) of embedding dimension 20, a max pooling layer, and a fully
connected layer for classification. The model is trained with the Adam
optimizer (Kingma and Ba, 2014) and a learning rate of 0.001 for 1000 epochs.
The datasets are split into 80%/10%/10% for training/validation/testing with
the model accuracy shown in Table 3.

Table 3. Accuracy of the GNN graph classifier.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
Training | 0.8439 | 0.8825 | 0.9980 | 0.7800
Validation | 0.8161 | 0.8302 | 0.9727 | 0.8198
Testing | 0.7809 | 0.8000 | 0.9781 | 0.7297

####  4.1.3. Baselines

To the best of our knowledge, GCFExplainer is the first global counterfactual
explainer. To validate its effectiveness, we compare it against state-of-the-
art local counterfactual explainers combined with the greedy summarization
algorithm described in Section 3.3. The following local counterfactual
generation methods are included in our experiments.

  * •

Ground-Truth: Using graphs belonging to the desired class from the original
dataset as local counterfactuals.

  * •

RCExplainer (Bajaj et al., 2021): Local counterfactual explainer based on the
modeling of implicit decision regions of GNNs.

  * •

CFF (Tan et al., 2022): Local counterfactual explainer based on joint modeling
of factual and counterfactual reasoning.

####  4.1.4. Explainer settings

We use a distance threshold $\theta$ of 0.05 for training all explainers. For
GCFExplainer, we set the teleportation probability $\tau=0.1$ and tune
$\alpha$, the weight between individual coverage and gain of coverage, from
$\\{0,0.5,1\\}$. A sensitivity analysis is presented in Section 4.6. The
number of VRRW iterations $M$ is set to 50000, which is enough for convergence
as shown in Section 4.5 For baselines, we tune their hyperparameters to
achieve the best local counterfactual rates while maintaining an average
distance to input graphs that is smaller than the distance threshold $\theta$.

Figure 4. Recourse coverage comparison between GCFExplainer and baselines
based on different distance threshold values ($\theta$). GCFExplainer
consistently outperforms the baselines across different $\theta$.

###  4.2. Recourse Quality

We start by comparing the recourse quality between GCFExplainer and baselines.
Table 2 shows the recourse coverage with $\theta=0.1$ and median recourse cost
of the top 10 counterfactual graphs (i.e., $k=10$). We first notice that the
two state-of-the-art local counterfactual explainers have similar performance
as Ground-Truth, consistent with our claim that local counterfactual examples
from existing methods are not good candidates for a global explanation. The
proposed GCFExplainer, on the other hand, achieves significantly better
performance for global recourse quality. Compared to the best baseline,
RCExplainer, GCFExplainer realizes a 46.9% gain in recourse coverage and a
9.5% reduction in recourse cost.

Next, we show the recourse coverage and cost for different sizes of
counterfactual summary in Figure 3. As expected, adding more graphs to the
recourse representation increases recourse coverage while decreasing recourse
cost, at the cost of interpretability. And GCFExplainer maintains a constant
edge over the baselines.

We also compare the recourse coverage based on different distance thresholds
$\theta$, with results shown in Figure 4. While coverage increases for all
methods as the threshold increases, GCFExplainer consistently outperforms the
baselines across different sizes.

###  4.3. Global Counterfactual Insight

We have demonstrated the superiority of GCFExplainer based on various quality
metrics for global recourse. Here, we show how GCFExplainer provides global
insights compared to local counterfactual examples. Figure 6 illustrates (a)
four input undesired graphs with a similar structure from the AIDS dataset,
(b) corresponding local counterfactual examples (based on RCExplainer and
CFF), and (c) the representative global counterfactual graph from GCFExplainer
covering the input graphs. Our goal is to understand why the input graphs are
inactive against AIDS (undesired) and how to obtain the desired property with
minimal changes.

The local counterfactuals in (b) attribute the classification results to
different edges in individual graphs (shown as red dotted lines) and recommend
their removal to make input graphs active against HIV. Note that while only
two edits are proposed for each individual graph, they appear at different
locations, which are hard to generalize for a global view of the model
behavior. In contrast, the global counterfactual graph from GCFExplainer
presents a high-level recourse rule. Specifically, the carbon atom with the
carbon-oxygen bond is connected to _two_ other carbon atoms in the input
graphs, making them ketones (with a C=O bond) or ethers (with a C-O bond). On
the other hand, the global counterfactual graph highlights a different
functional group, aldehyde (shown in blue), to be the key for combating AIDS.
In aldehydes, the carbon atom with a carbon-oxygen bond is only connected to
_one_ other carbon atom, leading to different chemical properties compared to
ketones and ethers. Indeed, aldehydes have been shown to be effective HIV
protease inhibitors (Sarubbi et al., 1993).

Finally, this case study also demonstrates that counterfactual candidates
found by GCFExplainer are better for global explanation than local
counterfactuals. We note that while the graph edit distance between the local
counterfactuals and their corresponding input graphs is only 2, they do not
cover other similarly structured input graphs (with distance $>5$). Meanwhile,
our global counterfactual graph covers all input graphs (with distance $\leq
4$).

###  4.4. Ablation Study

We then conduct an ablation study to investigate the effectiveness of
GCFExplainer components. We consider three alternatives:

  * •

GCFExplainer-NVR: no vertex-reinforcement ($N(v)=1$)

  * •

GCFExplainer-NIF: no importance function ($I(v)=1$)

  * •

GCFExplainer-NDT: no dynamic teleportation ($p_{\tau}(G)=1/|\mathbb{G}|$)

The coverage results are shown in Table 4. We observe decreased performance
when any of GCFExplainer components is absent.

Table 4. Ablation study results based on recourse coverage.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
GCFExplainer-NVR | 24.56% | 35.44% | 11.33% | 8.56%
GCFExplainer-NIF | 13.29% | 29.16% | 4.54% | 6.83%
GCFExplainer-NDT | 27.34% | 36.35% | 14.05% | 9.28%
GCFExplainer | 27.85% | 37.08% | 14.66% | 10.93%

###  4.5. Convergence Analysis

In this subsection, we show the empirical convergence of VRRW based on the
mutagenicity dataset in Figure 5. We observe that the coverage performance for
different summary sizes starts to converge after 15000 iterations and fully
converges after 50000 iterations, which is the number we applied in our
experiments.

Figure 5. Convergence of VRRW for the mutagenicity dataset based on recourse
coverage with different summary sizes. VRRW fully converges after $M=50000$
iterations.

C(-[:-18,1.5,,,,,,])([:180]-C(-[:50,2.2,,,,,,])=O)*5(-C([:90,1.0]-C?)-[,,,,thick,dashdotted,,brown,draw=none]C([:45]-C)([:-45]-C)-C?-C([:90]-C)-[,,,,thick,dashdotted,,brown,draw=none])
C([:180]-C(-[:30,1.6,,,,,,])=O)*5(-@bC([:90,1.0,,,thick,dashdotted,,brown,draw=none]-C?)-C([:45]-C)([:315]-C)-C?-C([:90]-@aC)-[,,,,thick,dashdotted,,brown,draw=none])
[-](a)..controls +(-30:3em) and +(60:3em)..(b);
C([:180]-C=O)*5(-[,,,,thick,dashdotted,,brown,draw=none]C(-[:155,1.6,,,,,,])([:90,1.0]-C?)-C([:45]-C)([:315]-C)-[,,,,thick,dashdotted,,brown,draw=none]C?-C([:90]-C)-)
C([:180]-C(-[:26,1.6,,,,,,])-O(-[:18,2.5,,,,,,]))*5(-C([:90,1.0]-C?)-C([:45]-C)([:315]-C)-C?-C([:90]-C)-[,,,,thick,dashdotted,,brown,draw=none])
(a) Input graphs

C(-[:-18,1.4,,,,,,])([:180]-C(-[:50,2.2,,,,,,])=O)*5(-C([:90,1.0]-C?)-[,,,,draw=none]C([:45]-C)(-[:315,1.1,,,thick,dotted,,red]C)-C?-C(-[:90,,,,thick,dotted,,red]C))
C([:180]-C(-[:30,1.6,,,,,,])=[,,,,thick,dotted,,red]O)*5(-@bC([:90,1.0,,,thick,dashdotted,,brown,draw=none]-C?)-C([:45]-C)([:315]-C)-C?-C([:90]-@aC)-[,,,,thick,dashdotted,,brown,draw=none])
[-,thick,dotted,red](a)..controls +(-30:3em) and +(60:3em)..(b);
C([:180]-C=O)*5(-[,,,,thick,dashdotted,,brown,draw=none]C(-[:155,1.6,,,,,,])(-[:90,1.0,,,thick,dotted,,red]C?)-C([:45]-C)([:315]-C)-[,,,,thick,dashdotted,,brown,draw=none]C?-C(-[:90,,,,thick,dotted,,red]C)-)
C([:180]-C(-[:26,1.6,,,,,,])-O(-[:18,2.5,,,,,,]))*5(-C(-[:90,1.0,,,thick,dotted,,red]C?)-C([:45]-C)(-[:315,1.1,,,thick,dotted,,red]C)-C?-C([:90]-C)-[,,,,thick,dashdotted,,brown,draw=none])
(b) Local counterfactuals

C([:180]-[,,,,thick,,,blue]C=[,,,,thick,,,blue]O)*5(-C(-[:90,1.0]C?)-C([:45]-C)([:315]-C)-C?-C([:90]-C)-)
(c) Global counterfactual

Figure 6. Illustration of global and local counterfactual explanations for the
AIDS dataset. The global counterfactual graph (c) presents a high-level
recourse rule—changing ketones and ethers into aldehydes (shown in blue)—to
combat HIV, while the edge removals (shown in red) recommended by local
counterfactual examples (b) are hard to generalize.

###  4.6. Sensitivity Analysis

The only hyperparameter of GCFExplainer we tune is $\alpha$ in Equation 3 that
weights the individual coverage and gain of coverage for the importance
function. Table 6 shows the results based on different $\alpha$. While
GCFExplainer outperforms baselines with all different $\alpha$, we observe
that individual coverage works better for NCI1 and gain of cumulative coverage
works better for other datasets.

###  4.7. Running Time

Table 5 summarizes the running times of generating counterfactual candidates
based on different methods. GCFExplainer has a competitive running time albeit
exploring more counterfactual graphs in the process. We also include results
for GCFExplainer-S which samples a maximum of 10000 neighbors for computing
the importance at each step. It achieves better running time at a negligible
cost of 3.3% performance loss on average. Finally, summarizing the
counterfactual candidates takes less than a second for all methods.

Table 5. Counterfactual candidates generation time comparison. GCFExplainer (-S) has competitive running time albeit exploring more counterfactual graphs.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
RCExplainer | 30454 | 52549 | 29047 | 8444
CFF | 22794 | 31749 | 21296 | 6412
GCFExplainer | 19817 | 24006 | 2615 | 19246
GCFExplainer-S | 19365 | 18798 | 2539 | 7429
Table 6. Sensitivity analysis on $\alpha$, the weight between individual coverage and gain of coverage in the importance function.  | NCI1 | Mutagenicity | AIDS | Proteins
---|---|---|---|---
$\alpha=0.0$ | 27.85% | 36.87% | 12.83% | 10.11%
$\alpha=0.5$ | 27.50% | 36.59% | 14.66% | 10.38%
$\alpha=1.0$ | 22.27% | 37.08% | 13.99% | 10.93%

##  5\. Related Work

Explanations for Graph Neural Networks. There is much research (Ying et al.,
2019; Luo et al., 2020; Vu and Thai, 2020; Yuan et al., 2020) on explaining
graph neural networks (GNNs). The first proposed method, GNNExplainer (Ying et
al., 2019), finds the explanatory subgraph and sub-features by maximizing the
mutual information between the original prediction and the prediction based on
the subgraph and sub-features. Later, PGExplainer (Luo et al., 2020) provides
an inductive framework that extracts GNN node embeddings and learns to map
embedding pairs to the probability of edge existence in the explanatory
weighted subgraph. PGMExplainer (Vu and Thai, 2020) builds a probabilistic
explanation model that learns new predictions from perturbed node features,
performs variable selection using Markov blanket of variables, and then
produces a Bayesian network via structure learning. In XGNN (Yuan et al.,
2020), the authors find model-level explanations by a graph generation module
that outputs a sequence of edges using reinforcement learning. These
explanation methods focus on _factual_ reasoning while the goal of our work is
to provide a global _counterfactual_ explanation for GNNs.

Counterfactual Explanations. Recently, there are several attempts to have
explanations of graph neural networks (GNNs) via counterfactual reasoning
(Lucic et al., 2022; Bajaj et al., 2021; Abrate and Bonchi, 2021; Tan et al.,
2022). One of the earlier methods, CF-GNNExplainer (Lucic et al., 2022),
provides counterfactual explanations in terms of a learnable perturbed
adjacency matrix that leads to the flipping of classifier prediction for a
node. On the other hand, RCExplainer (Bajaj et al., 2021) aims to find a
robust subset of edges whose removal changes the prediction of the remaining
graph by modeling the implicit decision regions based on GNN graph embeddings.
In (Abrate and Bonchi, 2021), the authors investigate counterfactual
explanations for a more specific class of graphs—the brain networks—that share
the same set of nodes by greedily adding or removing edges using a heuristic.
More recently, the authors of CFF (Tan et al., 2022) argue that a good
explanation for GNNs should consider both factual and counterfactual reasoning
and they explicitly incorporate those objective functions when searching for
the best explanatory subgraphs and sub-features. Counterfactual reasoning has
also been applied for link prediction (Zhao et al., 2022). All the above
methods produce local counterfactual examples while our work aims to provide a
_global_ explanation in terms of a summary of representative counterfactual
graphs.

##  6\. Conclusion

We have proposed GCFExplainer, the first global counterfactual explainer for
graph classification. Compared to local explainers, GCFExplainer provides a
high-level picture of the model behavior and effective global recourse rules.
We hope that our work will not only deepen our understanding of graph neural
networks but also build a bridge for experts from other domains to leverage
deep learning models for high-stakes decision-making.

###### Acknowledgements.

This work is partially funded by NSF via grant IIS 1817046. The authors thank
Sevgi Kosan and Ilhan Kosan for their helpful comments on the chemical
properties in the global counterfactual analysis.

## References

  * (1)
  * Abrate and Bonchi (2021) Carlo Abrate and Francesco Bonchi. 2021.  Counterfactual graphs for explainable classification of brain networks. In _SIGKDD_.
  * Bajaj et al. (2021) Mohit Bajaj, Lingyang Chu, Zi Yu Xue, Jian Pei, Lanjun Wang, Peter Cho-Ho Lam, and Yong Zhang. 2021.  Robust Counterfactual Explanations on Graph Neural Networks. In _NeurIPS_.
  * Borgwardt et al. (2006) Karsten Borgwardt, Nicol Schraudolph, and SVN Vishwanathan. 2006\.  Fast computation of graph kernels. In _NeurIPS_.
  * Borgwardt et al. (2005) Karsten M Borgwardt, Cheng Soon Ong, Stefan Schönauer, SVN Vishwanathan, Alex J Smola, and Hans-Peter Kriegel. 2005.  Protein function prediction via graph kernels.  _Bioinformatics_ 21, suppl_1 (2005), i47–i56.
  * Costa and De Grave (2010) Fabrizio Costa and Kurt De Grave. 2010.  Fast neighborhood subgraph pairwise distance kernel. In _ICML_.
  * Dobson and Doig (2003) Paul D Dobson and Andrew J Doig. 2003.  Distinguishing enzyme structures from non-enzymes without alignments.  _Journal of molecular biology_ 330, 4 (2003), 771–783.
  * Hamilton et al. (2017) Will Hamilton, Zhitao Ying, and Jure Leskovec. 2017\.  Inductive representation learning on large graphs.  _NeurIPS_ (2017).
  * Huang et al. (2021) Zexi Huang, Arlei Silva, and Ambuj Singh. 2021.  A broader picture of random-walk based graph embedding. In _SIGKDD_.
  * Huang et al. (2022) Zexi Huang, Arlei Silva, and Ambuj Singh. 2022.  POLE: Polarized Embedding for Signed Networks. In _WSDM_.
  * Jiang et al. (2020) Mingjian Jiang, Zhen Li, Shugang Zhang, Shuang Wang, Xiaofeng Wang, Qing Yuan, and Zhiqiang Wei. 2020.  Drug–target affinity prediction using graph neural network and contact maps.  _RSC advances_ 10, 35 (2020), 20701–20712.
  * Kazius et al. (2005) Jeroen Kazius, Ross McGuire, and Roberta Bursi. 2005\.  Derivation and validation of toxicophores for mutagenicity prediction.  _Journal of medicinal chemistry_ 48, 1 (2005), 312–320.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv:1412.6980_ (2014).
  * Kipf and Welling (2017) Thomas N Kipf and Max Welling. 2017.  Semi-supervised classification with graph convolutional networks. In _ICLR_.
  * Klicpera et al. (2018) Johannes Klicpera, Aleksandar Bojchevski, and Stephan Günnemann. 2018.  Predict then propagate: Graph neural networks meet personalized pagerank. In _ICLR_.
  * Kosan et al. (2021) Mert Kosan, Arlei Silva, Sourav Medya, Brian Uzzi, and Ambuj Singh. 2021.  Event detection on dynamic graphs.  _arXiv preprint arXiv:2110.12148_ (2021).
  * Liang and Zhao (2017) Yongjiang Liang and Peixiang Zhao. 2017.  Similarity Search in Graph Databases: A Multi-Layered Indexing Approach. In _ICDE_. 783–794.
  * Lucic et al. (2022) Ana Lucic, Maartje A Ter Hoeve, Gabriele Tolomei, Maarten De Rijke, and Fabrizio Silvestri. 2022\.  Cf-gnnexplainer: Counterfactual explanations for graph neural networks. In _AISTATS_.
  * Luo et al. (2020) Dongsheng Luo, Wei Cheng, Dongkuan Xu, Wenchao Yu, Bo Zong, Haifeng Chen, and Xiang Zhang. 2020.  Parameterized explainer for graph neural network. In _NeurIPS_.
  * Mei et al. (2010) Qiaozhu Mei, Jian Guo, and Dragomir Radev. 2010.  Divrank: the interplay of prestige and diversity in information networks. In _SIGKDD_.
  * Metwally et al. (2005) Ahmed Metwally, Divyakant Agrawal, and Amr El Abbadi. 2005\.  Efficient computation of frequent and top-k elements in data streams. In _ICDT_.
  * Mirhoseini et al. (2020) Azalia Mirhoseini, Anna Goldie, Mustafa Yazgan, Joe W. J. Jiang, Ebrahim M. Songhori, Shen Wang, Young-Joon Lee, Eric Johnson, Omkar Pathak, Sungmin Bae, Azade Nazi, Jiwoo Pak, Andy Tong, Kavya Srinivasa, William Hang, Emre Tuncer, Anand Babu, Quoc V. Le, James Laudon, Richard Ho, Roger Carpenter, and Jeff Dean. 2020\.  Chip Placement with Deep Reinforcement Learning.  _CoRR_ abs/2004.10746 (2020).
  * Natarajan and Ranu (2016) Dheepikaa Natarajan and Sayan Ranu. 2016.  A scalable and generic framework to mine top-k representative subgraph patterns. In _ICDM_.
  * Pemantle (1992) Robin Pemantle. 1992\.  Vertex-reinforced random walk.  _Probability Theory and Related Fields_ 92, 1 (1992), 117–136.
  * Perozzi et al. (2014) Bryan Perozzi, Rami Al-Rfou, and Steven Skiena. 2014\.  Deepwalk: Online learning of social representations. In _SIGKDD_.
  * Rawal and Lakkaraju (2020) Kaivalya Rawal and Himabindu Lakkaraju. 2020.  Beyond individualized recourse: Interpretable and interactive summaries of actionable recourses. In _NeurIPS_.
  * Riesen and Bunke (2008) Kaspar Riesen and Horst Bunke. 2008.  IAM graph database repository for graph based pattern recognition and machine learning. In _Joint IAPR International Workshops on Statistical Techniques in Pattern Recognition (SPR) and Structural and Syntactic Pattern Recognition (SSPR)_. Springer, 287–297.
  * Sanfeliu and Fu (1983) Alberto Sanfeliu and King-Sun Fu. 1983.  A distance measure between attributed relational graphs for pattern recognition.  _IEEE transactions on systems, man, and cybernetics_ 3 (1983), 353–362.
  * Sarubbi et al. (1993) Edoardo Sarubbi, Pier Fausto Seneci, Michael R Angelastro, Norton P Peet, Maurizio Denaro, and Khalid Islam. 1993.  Peptide aldehydes as inhibitors of HIV protease.  _FEBS letters_ 319, 3 (1993), 253–256.
  * Shervashidze et al. (2011) Nino Shervashidze, Pascal Schweitzer, Erik Jan Van Leeuwen, Kurt Mehlhorn, and Karsten M Borgwardt. 2011\.  Weisfeiler-lehman graph kernels.  _JMLR_ 12, 9 (2011).
  * Tan et al. (2022) Juntao Tan, Shijie Geng, Zuohui Fu, Yingqiang Ge, Shuyuan Xu, Yunqi Li, and Yongfeng Zhang. 2022.  Learning and evaluating graph neural network explanations based on counterfactual and factual reasoning. In _WebConf_.
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. 2018.  Graph Attention Networks. In _ICLR_.
  * Voigt and Von dem Bussche (2017) Paul Voigt and Axel Von dem Bussche. 2017.  The eu general data protection regulation (gdpr).  _A Practical Guide, 1st Ed., Cham: Springer International Publishing_ 10, 3152676 (2017), 10–5555.
  * Vu and Thai (2020) Minh Vu and My T Thai. 2020\.  Pgm-explainer: Probabilistic graphical model explanations for graph neural networks. In _NeurIPS_.
  * Wale and Karypis (2006) Nikil Wale and George Karypis. 2006.  Comparison of Descriptor Spaces for Chemical Compound Retrieval and Classification. In _ICDM_.
  * Wang et al. (2022) Yuke Wang, Boyuan Feng, and Yufei Ding. 2022.  QGTC: accelerating quantized graph neural networks via GPU tensor core. In _PPoPP_.
  * Wang et al. (2021) Yuke Wang, Boyuan Feng, Gushu Li, Shuangchen Li, Lei Deng, Yuan Xie, and Yufei Ding. 2021.  GNNAdvisor: An Efficient Runtime System for GNN Acceleration on GPUs. In _OSDI_.
  * Xiong et al. (2021) Jiacheng Xiong, Zhaoping Xiong, Kaixian Chen, Hualiang Jiang, and Mingyue Zheng. 2021.  Graph neural networks for automated de novo drug design.  _Drug Discovery Today_ 26, 6 (2021), 1382–1393.
  * Yang et al. (2011) Shuang-Hong Yang, Bo Long, Alex Smola, Narayanan Sadagopan, Zhaohui Zheng, and Hongyuan Zha. 2011\.  Like like alike: joint friendship and interest propagation in social networks. In _WebConf_.
  * Ying et al. (2019) Rex Ying, Dylan Bourgeois, Jiaxuan You, Marinka Zitnik, and Jure Leskovec. 2019.  Gnnexplainer: Generating explanations for graph neural networks. In _NeurIPS_.
  * Ying et al. (2018) Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L. Hamilton, and Jure Leskovec. 2018.  Graph Convolutional Neural Networks for Web-Scale Recommender Systems. In _SIGKDD_.
  * Yuan et al. (2020) Hao Yuan, Jiliang Tang, Xia Hu, and Shuiwang Ji. 2020\.  XGNN: Towards Model-Level Explanations of Graph Neural Networks. In _SIGKDD_.
  * Zhao et al. (2022) Tong Zhao, Gang Liu, Daheng Wang, Wenhao Yu, and Meng Jiang. 2022.  Learning from Counterfactual Links for Link Prediction. In _ICML_.


# Accurate Spectral Collocation Computation of High Order Eigenvalues for
Singular Schrödinger Equations

Călin-Ioan Gheorghiu Tiberiu Popoviciu Institute of Numerical Analysis,
Romanian Academy<EMAIL_ADDRESS>
###### Abstract

We are concerned with the study of some classical spectral collocation methods
as well as with the new software system Chebfun in computing high order
eigenpairs of singular and regular Schrödinger eigenproblems. We want to
highlight both the qualities as well as the shortcomings of these methods and
evaluate them in conjunction with the usual ones. In order to resolve a
boundary singularity we use Chebfun with domain truncation. Although it is
applicable with spectral collocation, a special technique to introduce
boundary conditions as well as a coordinate transform, which maps an unbounded
domain to a finite one, are the special ingredients. A challenging set of
“hard”benchmark problems, for which usual numerical methods (f. d., f. e. m.,
shooting etc.) fail, are analyzed. In order to separate “good”and
“bad”eigenvalues we estimate the drift of the set of eigenvalues of interest
with respect to the order of approximation and/or scaling of domain parameter.
It automatically provides us with a measure of the error within which the
eigenvalues are computed and a hint on numerical stability. We pay a
particular attention to problems with almost multiple eigenvalues as well as
to problems with a mixed spectrum.

Keywords: spectral collocation; Chebfun; singular Schrödinger; high index
eigenpairs; multiple eigenpairs; accuracy; numerical stability.

MSC2010: 34L40; 34L16; 65L15; 65L20; 65L60

##  1 Introduction

There is clearly an increasing interest to develop accurate and efficient
methods of solution to singular Schrödinger eigenproblems. Our main interest
here is to compare the capabilities of the new Chebfun package with those of
classical spectral methods in solving such problems, showing their
capabilities and weaknesses. The latter employ basis functions and/or grid
points based on Chebyshev, Laguerre or Hermite polynomials as well as on sinc
or Fourier functions. The effort expended by both classes of methods is also
of real interest. It can be assessed in terms of the ease of implementation of
the methods as well as in terms of computer resources required to achieve a
specified accuracy.

Spectral methods have been shown to provide exponential convergence for a
large variety of problems, generally with smooth solutions, and are often
preferred. For details on Chebfun we refer to [1], [2], [3], [4], [5] and [6].
For Chebyshev collocation (ChC), Laguerre-Gauss-Radau collocation (LGRC),
Hermite (HC) and sinc collocation (SiC) we refer among other sources to our
contributions [7], [8] as well as to the seminal paper [9].

For problems on the entire real axis SiC proved to be particularly well
suited. Moreover, this method has given excellent results recorded in our
contribution [8] and in the works cited there. The so-called generalized
pseudospectral (GPS) method, actually the Legendre collocation, is employed in
[10] to calculate the bound states of the Hulthén and the Yukawa potentials in
quantum mechanics, with special emphasis on higher excited states and stronger
couplings. The author uses a two parameter dependent nonlinear transformation
in order to map the half-line into the canonical interval $\left[-1,1\right].$
In contrast, we will use an analogous transformation but which depends on only
one parameter. Also, very recently spectral methods based on non-classical
orthogonal polynomials have been used in [11] in order to solve some
Schrödinger problems connected with Fokker-Planck operator. A particular
attention will be paid in this paper to the challenging issue of continuous
spectra vs. discrete (numerical) eigenvalues. It is well known that some
eigenvalue problems (see for instance the well known text [12]) for
differential operators which are naturally posed on the whole real line or the
half-line, often lead to some discrete eigenvalues plus a continuous spectrum.
Actually, the usual numerical approximation typically involves three
processes:

  1. 1.

reduction to a finite interval;

  2. 2.

discretization;

  3. 3.

application of a numerical eigenvalue solver.

Reduction to a finite interval and discretization typically eliminate the
continuous spectrum. Even if we do not truncate a priory the domain on which
the problem is formulated such an inherent reduction can not be avoided.

It can be argued that, generally speaking, in solving various differential
problems, the Chebfun software provides a greater flexibility than the
classical spectral methods. This fact is fully true for regular problems.

Unfortunately, in the presence of various singularities, the maximum order of
approximation $N,$ of the unknowns can be reached ($N\geq 4000$) and then
Chebfun issues a message that warns about the possible inaccuracy of the
results provided.

We came out of this tangle using modified classical spectral methods. In this
way, when we had serious doubts about the accuracy of the solutions given by
Chebfun, we managed to establish the correctness of the numerical results.

As a matter of fact, in order to resolve a singularity on the ends of an
unbounded integration interval, Chebfun uses only the arbitrary truncation of
the domain. Classical spectral methods can also use this method, but it is not
recommended. For singular points at finite distances (mainly origin) we will
use the so-called removing technique of independent boundary conditions. The
boundary conditions at infinity can be enforced using basis functions that
satisfy these conditions (Laguerre, Hermite, sinc). An alternative method,
which proved to be very accurate, is the Chebyshev collocation (ChC) method in
combination with a change of variables (coordinates) which transform the half
line into the canonical Chebyshev interval $\left[-1,1\right].$ Then the
removing technique of independent boundary conditions is essential in order to
remove the singularities at the end points of integration interval.

A Chebfun code and two MATLAB codes, one for ChC and another for SiC method,
are provided in order to exemplify. With minor modifications they could be
fairly useful for various numerical experiments.

The structure of this work is as follows. In Section 2 we recall some specific
issues for the regular as well as singular Schrödinger eigenproblems. The main
comment refers to the notion of mixed spectrum. In Section 3 we review on the
Chebfun structure and the classical spectral methods (differentiation
matrices, enforcing boundary conditions, etc.). The Section 4 is the central
part of the paper. Here we analyze some benchmark problems. In order to
separate the “good”from the “bad”eigenvalues we estimate their relative drift
with respect to some parameters. The accuracy in computing eigenfunctions is
estimated by their departure from orthogonality. We end up with Section 5
where we underline some conclusions and suggest some open problems.

##  2 Regular and singular Schrödinger eigenproblems

The Schrödinger equation reads

$u^{\prime\prime}+\left[\lambda-q\left(x\right)\right]u=0,\ -\infty\leq
a<x<b\leq+\infty.$ (1)

It is a Liouville normal form of a general Sturm-Liouville (SL) equation where
$\lambda$ is proportional with the energy levels of the physical system,
$q\left(x\right)$ is directly proportional with the potential energy and the
“wave function”$u$ may be real or complex such that $u\
u^{\ast}dx=\left|u\right|^{2}dx$ is the probability that the particle under
consideration will be “observed”in the interval $\left(x,\ x+dx\right).$ In
problems involving Schrödinger equations, it is customary among chemists and
physicists to define the spectrum of this Sturm-Liouville problem as the all
eigenvalues $\lambda$ for which eigenfunctions $u$ exist. The set of isolated
points (if any) in this spectrum is called the discrete spectrum; the part (if
any) that consists of entire interval is called continuous spectrum. We shall
adopt this suggestive terminology here.

The equation (1) can be given on a finite, semi-infinite, or infinite
interval. Only on a closed and finite interval $a\leq x\leq b$ can Schrödinger
equation be associated with a regular Sturm-Liouville problem. If the interval
of definition is semi-infinite or infinite, or is finite and $q\left(x\right)$
vanishes at one or both endpoints, or if $q$ is discontinuous, we can not
obtain from (1) a regular Sturm-Liouville problem. In any such case, the
Schrödinger equation (1) is called singular. We obtain a singular eigenproblem
from a singular Schrödinger equation by imposing suitable homogeneous boundary
conditions. They can not always be described by formulae like $\alpha
u\left(e\right)+\beta u^{\prime}\left(e\right)=0,$ where $e$ can be the end
point $a$ or $b.$ For instance, the condition that $u$ be bounded near a
singular end point, which can be finite or $\pm\infty,$ is a common boundary
condition defining a singular eigenproblem. For regular SL problems, it is
proved (see for instance [12]) that the spectrum is always discrete, and the
eigenfunctions are (trivially) square-integrable. For singular problems the
situation is completely different. For instance in their textbook [12]
Birkhoff and Rota consider the eigenproblem attached to the free particle
equation and show that the spectrum of the free particle is continuous. Some
software packages have been designed over time to solve various singular SL
problems. The most important would be SLEIGN and SLEIGN2, SLEDGE, SL02F and
MATSLISE. The SLDRIVER interactive package supports exploration of a set of SL
problems with the four previously mentioned packages. In [13] (see also [14])
the authors designed the software package SLEDGE. They observed that for a
class of singular problems their method either fails or converges very slowly.
Essentially, the numerical method used in this software package replaces the
coefficient function $q(x)$ by step function approximation. Similar behavior
has been observed on the NAG code SL02F introduced in [15] and [16] as well as
on the packages SLEIGN and SLEIGN2 introduced in [17] and [18]. The MATSLISE
code introduced in [19] can solve some Schrödinger eigenvalue problem by a
constant perturbation method of a higher order.

The main purpose of this paper is to argue that Chebfun, along with the
spectral collocation methods, can be a very feasible alternative to these
software packages regarding accuracy, robustness as well as simplicity of
implementation. In addition, these methods can compute the “whole”set of
eigenvectors and provide some details on the accuracy and numerical stability
of the results provided.

##  3 Chebfun vs. spectral collocation (ChC, LGRC, SiC)

###  3.1 Chebfun

The Chebfun system, in object-oriented MATLAB, contains algorithms which
amount to spectral collocation methods on Chebyshev grids of automatically
determined resolution. Its properties are briefly summarized in [2]. In [1]
the authors explain that chebops are the fundamental Chebfun tools for solving
ordinary differential (or integral) equations. One may then use them as tools
for more complicated computations that may be nonlinear and may involve
partial differential equations. This is analogous to the situation in MATLAB
itself. The implementation of chebops combines the numerical analysis idea of
spectral collocation with the computer science idea of lazy or delayed
evaluation of the associated spectral discretization matrices. The grammar of
chebops along with a lot of illustrative examples is displayed in the above
quoted paper as well as in the text [5]. Thus one can get a suggestive image
of what they can do.

In [1] p.12 the authors explain clearly how the Chebfun works, i.e., it solves
the eigenproblem for two different orders of approximation, automatically
chooses a reference eigenvalue and checks the convergence of the process. At
the same time, it warns about the possible failures due to the high non-
normality of the analyzed operator (matrix).

Actually, we want to show in this paper that Chebfun along with chebops can do
much more, i.e., can accurately solve highly (double) singular Schrödinger
eigenproblems.

###  3.2 ChC, LGRC and SiC

In the spectral collocation method the unknown solution to a differential
equation is expanded as a global interpolant, such as a trigonometric or
polynomial interpolant. In other methods, such as f. e. and f. d., the
underlying expansion involves local interpolants such as piecewise
polynomials. This means that the accuracy of spectral collocation is superior.
For problems with smooth solutions convergence rates are typically of order
$e^{-cN}$ or $e^{-c\sqrt{N}}$ where $N$ is the order of approximation or
resolution, i.e., the number of degree of freedom in expansion. In contrast,
f. e. or f. d. yield convergence rates that are only algebraic in $N$,
typically of orders $N^{-2}$ or $N^{-4}.$ The net superiority of global
spectral methods on local methods is discussed in detail in [20]. In all
spectral collocation methods designed so far we have used the collocation
differentiation matrices from the seminal paper [9]. We preferred this MATLAB
differentiation suite for the accuracy, efficiency as well as for the
ingenious way of introducing various boundary conditions.

In order to impose (enforce) the boundary conditions we have used two methods
that are conceptually different, namely the boundary bordering as well as the
basis recombination. A very efficient way to accomplish the boundary bordering
is available in [21] and is called removing technique of independent boundary
conditions. We have used this technique in the large majority of our papers
except [22] where the latter technique has been employed. In the last quoted
paper a modified Chebyshev tau method based on basis recombination has been
used in order to solve an Orr-Sommerfeld problem with an eigenparameter
dependent boundary condition. Even in eigenproblems that contain the spectral
parameter in the definition of the boundary conditions we have used the
boundary bordering technique (see [23]).

In [24] (see also [25]) we have solved some multiparameter eigenproblems (MEP)
which come from separation of variables, in several orthogonal coordinate
systems, applied to the Helmholtz, Laplace, or Schrödinger equation. Important
cases include Mathieu’s system, Lame’s system, and a system of spheroidal wave
functions. We show that by combining spectral collocation methods, ChC and
LGRC, and new efficient numerical methods for solving algebraic MEPs, it is
possible to solve such problems both very efficiently and accurately. We
improve on several previous results available in the literature, and also
present a MATLAB toolbox for solving a wide range of problems.

###  3.3 The drift of eigenvalues

Two techniques are used in order to eliminate the “bad”eigenvalues as well as
to estimate the stability (accuracy) of computations. The first one is the
drift, with respect to the order of approximation or the scaling factor, of a
set of eigenvalues of interest. In a simplified form this concept has been
introduced by J. P. Boyd in [26]. The second one is based on the check of the
eigenvectors’ orthogonality.

In other words, we want to separate the “good”eigenvalues from the “bad”ones,
i.e., inaccurate eigenvalues. An obvious way to achieve this goal is to
compare the eigenvalues computed for different orders of some parameters such
as the approximation order (cut-off parameter) $N$ or the scaling factor. Only
those whose difference or “resolution-dependent drift”is “small”can be
believed. Actually, in [26] the so called absolute (ordinal) drift with
respect to the order of approximation has been introduced.

We extend this definition to the following one. The absolute (ordinal) drift
of the $jth$ eigenvalue with respect to the parameter $\alpha$ is defined as

$\delta_{j,absolute,\alpha}:=\left|\lambda_{j}^{\left(\alpha_{1}\right)}-\lambda_{j}^{\left(\alpha_{2}\right)}\right|,\quad\alpha_{1}\neq\alpha_{2},$
(2)

where $\lambda_{j}^{\left(\alpha\right)}$ is the $jth$ eigenvalue, after the
eigenvalues have been sorted, as computed using a specific value of the
parameter. In the most common cases this parameter can be $N$ or $c.$

The dependence of $\delta_{j,absolute,\alpha},$ $j=1,2,\ldots,Ne,$ where $Ne$
is the number of analyzed eigenvalues, on the index (mode) $j$ will be
displayed in a log-linear plot. If we divide the right hand side of (2) by
$\left|\lambda_{j}^{\left(\alpha_{1}\right)}\right|$ we get the so called
relative drift denote by $\delta_{j,relative,\alpha}.$

##  4 Numerical benchmark problems and discussions

###  4.1 A regular Schrödinger eigenproblem

The bounded Coffey-Evans potential reads

$q\left(x\right):=-2\beta\cos\left(2x\right)+\beta^{2}\sin^{2}\left(2x\right),\
\beta\in{\mathbb{R}}.$ (3)

We attach to equation (1) the homogeneous Dirichlet boundary conditions
$u\left(\pm\pi/2\right)=0$ and use $\beta:=30.$ In spite of being regular, the
Coffey-Evans problem is one of the most difficult test problems in the
literature because there are very close eigenvalue triplets as $\beta$
increases.

We have used a short Chebfun code in order to solve this regular eigenproblem.
It is available in the next lines. With appropriate changes, this code can be
used to analyze any other Schrödinger eigenproblem.

    dom=[-pi/2,pi/2];
    x=chebfun(’x’,dom);beta=30; sigma=-1;
    L=chebop(dom);
    L.op =@(x,y) -diff(y,2)+(-2*beta*cos(2*x)+(beta*sin(2*x))^2)*y;
    L.rbc=0; L.lbc=0; N=201;
    [V,D]=eigs(L,N,sigma); D=diag(D)

The eigenvalues obtained by ChC and Chebfun are extremely close, practically
indistinguishable. They also compare very well with those computed in [27] by
some coefficient approximation methods of orders $2$, $4$ and $8$ and reported
in Table 3 of this paper.

The absolute drift reported in Fig. 1 means that we can compute the first
hundred eigenvalues with better accuracy than $10^{-10}$. Unfortunately no
accuracy analysis is reported in [27] (see also [19] and [28]).

Figure 1: The drift of the first $100$ eigenvalues of Coffey-Evans problem.
ChC used the orders of approximations $N_{1}:=256$ and $N_{2}:=512$-red dotted
line and respectively $N_{1}:=400$ and $N_{2}:=512$-green circled line.

The eigenvectors are either symmetric (the even ones) or anti-symmetric (the
odd ones). The first four of them are depicted in Fig. 2. They look fairly
smooth and satisfy the boundary conditions.

Figure 2: The first four eigenvectors of Coffey-Evans problem computed by
Chebfun.

The Chebyshev coefficients of the first four eigenvectors of Coffey-Evans
problem computed by Chebfun are displayed in Fig. 3. These coefficients
decrease sharply and smoothly to a rounding-off plateau below $10^{-15}$.
Roughly speaking this means they are computed with the machine precision.

Figure 3: The Chebyshev coefficients of the first four eigenvectors of Coffey-Evans problem computed by Chebfun. Table 1: High index eigenvalues of Schrödinger eigenproblem equipped with Coffey-Evans potential (3) computed by three different methods. $\;j$ |  $\;\lambda_{j}$ by Chebfun |  $\;\lambda_{j}$ computed in [27] |  $\;\lambda_{j}$ by ChC
---|---|---|---
$\;20$ | $\;\numprint{951.8788067958783}$ | $\;951.878806796591$ | $\;951.8788067965993$
$\;30$ | $\;1438.295244640637$ | $\;1438.295244640802$ | $\;1438.295244640797$
$\;40$ | $\;2146.405360539156$ | $\;2146.405360539854$ | $\;2146.405360539845$
$\;50$ | $\;3060.923491511540$ | $\;3060.923491511421$ | $\;3060.923491511401$
$\;100$ | $\;10653.52543568510$ | $\;10653.525435875921$ | $\;10653.52543587600$
$\;200$ | $\;40851.63764596094$ | $\;40851.637646050455$ | $\;40851.63764605047$

For this potential the first eigenvalue $\lambda_{0}$ is close to zero
(actually we have got $\lambda_{0}=\numprint{-6.254959429708980e-12}$) and
there are very close eigenvalue triplets
$\left(\lambda_{2},\;\lambda_{3},\;\lambda_{4}\right)$,
$\left(\lambda_{6},\;\lambda_{7},\;\lambda_{8}\right)$, $\ldots$ as $\beta$
increases. The common numeric part of the eigenvalues in the first triplet is
$2.31664929$ and that of the eigenvalues in the second is a little shorter,
i.e. $4.45283.$

###  4.2 Two singular Schrödinger eigenproblem on the half line

In this section we study the mixed spectrum of some Schrödinger eigenproblems
having a “potential well”dying out at infinity.

####  4.2.1 Hydrogen atom equation

The first example consists in the equation (1) equipped with the potential

$q\left(x\right):=-\frac{1}{x}+\frac{l\left(l+1\right)}{x^{2}},\
l\in{\mathbb{R}},$ (4)

along with the boundary conditions

$u(0)=0,\ u\rightarrow 0\ as\ x\rightarrow\infty.$ (5)

The problem (1)-(4)-(5) is clearly singular. We must mention from the
beginning that our numerical experiments performed with LGRC and Chebfun
together with the truncation of the domain did not produce satisfactory
results. In these conditions we have resorted to the mapped ChC method.

In order to implement this method we use the algebraic map

$x:=c\frac{1+s}{1-s},\ s\in\left[-1,1\right],\
x\in[0,+\infty),c\in{\mathbb{R}},c>0,$ (6)

which, for each $c$, transforms the interval $\left[-1,1\right]$ into the half
line, and its inverse. The parameter $c$ is free to be tuned for optimum
accuracy.

The mapping (6) has been introduced in [29] where its practical effects have
been discussed. The author observed that the convergence of the Chebyshev
expansion is governed by the closeness of the singularities of the function
being expanded to the expansion region. The major effect of such mapping is to
allow us to move the singularities further away from the expansion region. The
mapping may also weaken the effect of the singularities by modifying the
strength of the singularity as well as moving it. The value of the scaling
parameter $c$ used in our computation was chosen essentially by trial and
error along with the drift with respect to this parameter.

In order to write down any second order differential equation in independent
variable $s$, we need the following derivatives:

$\begin{array}[]{c}u^{\prime}\left(x\right)=u^{\prime}\left(s\right)\frac{1}{x^{\prime}_{s}},\
x^{\prime}_{s}=\frac{2c}{\left(1-s\right)^{2}},\\\
u^{\prime\prime}\left(x\right)=u^{\prime\prime}\left(s\right)\frac{1}{\left(x^{\prime}_{s}\right)^{2}}-u^{\prime}\left(s\right)\frac{x_{s}^{\prime\prime}}{\left(x^{\prime}_{x}\right)^{3}}.\end{array}$
(7)

Now it is easy to write the differential equation for $u\left(s\right)$ and to
attach to this new equation the homogeneous Dirichlet boundary conditions
$u\left(\pm 1\right)=0.$ As usual we implement these conditions by deleting
the first and the last rows and columns of the collocation matrix attached to
the left hand side of the equation (1)-(4). This is the most simplified
version of removing technique of the independent boundary conditions
introduced in [21]. With (7) the following MATLAB code solve this problem.

    % approximation order and parameter l
    N=1600; l=1;
    % number of displayed eigenvalues
    Ne=50;
    c=2;         % scaling factor
    % Chebyshev differentiation matrices (Weideman & Reddy)
    [st,D]=chebdif(N,2);
    % N-2 kept nodes; 1 and N are removed nodes
    k=2:N-1;  s=st(k);
    % enforced boundary conditions in differentiation matrices
    D2=D(k,k,2); D1=D(k,k,1);
    % collocation matrix of the system
    A=-diag((1-x).^4)*D2/(4*(c^2))+diag((1-x).^3)*D1/(2*(c^2))+...
        diag(l*(l+1)*((1-x).^2)./(((1+x).^2)*(c*2))-(1-x)./(c*(1+x)));
    % computed and sorted eigenpairs
    [U,S]=eig(A); S=diag(S); [t,o]=sort(S); S=S(o); U=U(:,o);
    disp(S(1:Ne))
    % Chebyshev coefficients of the first four eigenvectors by FCT
    Ucoeff=fcgltran(U(:,1:4),1);

Figure 4: The first four eigenvectors for hydrogen atom eigenproblem
(1)-(4)-(5) with $l:=1$ computed by mapped ChC with scaling factor $c:=2$ and
$N:=512.$

The first four vectors of the problem are displayed in Fig. 4. It is clear
that they satisfy both boundary conditions but in the right neighborhood of
origin they have a totally different behavior from the eigenvectors of regular
problems, i.e., they vanish out on continuous portions and not in discrete
points. However, they clearly approximate square-integrable eigenfunctions and
thus confirm some theoretical results proved in [12].

Figure 5: The coefficients of the first four eigenvectors for hydrogen atom
eigenproblem (1)-(4)-(5) with $l:=1$ computed by mapped ChC with scaling
factor $c:=2$ and $N:=1600.$

Their Chebyshev coefficients obtained using FCT (fast Chebyshev transform-see
[30] for details) are displayed in Fig.5. They decrease sharply and smoothly
to some limits, followed by a wide rounding-off plateau. For the first vector
(the rightmost one) this limit is around $10^{-13}$. It increases with the
index of the vector.

Figure 6: The absolute drift with respect to $N$ of the first $50$ negative
eigenvalues of hydrogen atom eigenproblem (1)-(4)-(5) computed by mapped ChC
with scaling factor $c:=2$ and $N_{1}:=1600$ and $N_{1}:=2048$ (red line). The
green circled line signifies the drift of eigenvalues computed by mapped ChC
using $N:=2048$ with respect to the exact eigenvalues
$\lambda_{n}=-1/\left(2n+4\right)^{2},$ $n=0,1,2,\ldots,50$ corresponding to
$l:=1.$

Roughly this means that we cannot hope for a better approximation than
something of the order $10^{-13}$ when computing the first eigenvector.

When $N:=2048$ mapped ChC has found
$\lambda_{0}=\numprint{-6.250000000166379e-02}$ which is a very good
approximation of $-1/16.$ The largest negative eigenvalue has been
$\lambda_{66}=\numprint{-7.865782521027431e-07}$ and the next eigenvalue, the
smallest positive has been computed as
$\lambda_{67}=\numprint{6.517056834998433e-05}.$ To see how the mapped ChC
method simulates the notion of continuous spectrum we will provide in Table 2
some significant eigenvalues.

Table 2: High index eigenvalues of hydrogen atom eigenproblem (1)-(4)-(5). $\;j$ |  $\;\lambda_{j}$ by mapped ChC
---|---
$\;70$ | $\;\numprint{3.346710010488799e-4}$
$\;80$ | $\;2.246853675452558e-3$
$\;90$ | $\;6.551029735506819e-3$
$\;100$ | $\;1.474442763764248e-2$
$\;110$ | $\;2.890880952550005e-2$
$\;120$ | $\;5.183912194697291e-2$

In the monograph [12] Sect. 18 it is proved that if the potential
$q\left(x\right)$ is continuous and has the asymptotic behaviour
$q\left(x\right)=\frac{A}{x}+O\left(\frac{1}{x^{2}}\right)$ as
$x\rightarrow\infty,$ for $\lambda>0$ the spectrum is continuous and the
eigenfunctions are not square-integrable and for $\lambda<0$ the spectrum is
discrete and the eigenfunctions are square-integrable. The numerical results
gathered around this problem plainly confirm this analytical result.

Actually the eigenvector corresponding to the first positive eigenvalue (the
sixty-seventh one) is depicted in the left panel of Fig. 7. It is hard to
believe that this could approximate a square-integrable eigenfunction. In the
right panel of the same figure we displayed the Chebyshev coefficients of the
corresponding eigenvector. The oscillations of these coefficients mimic the
steep oscillations of the eigenvector.

For the real scaling factor $c$ we have to mention that in case of eigenvalue
problems it can be adjusted only on the mathematical basis. Thus, it can be
tuned in order to:

  * •

improve the decaying rate of the coefficients of spectral expansions.

  * •

find as orthogonal as possible eigenvectors to an eigenproblem;

Figure 7: (a) The eigenvector corresponding to the first positive eigenvalue
of hydrogen atom eigenproblem (1)-(4)-(5) computed by mapped ChC with scaling
factor $c:=2$ and $N:=2024$, $l:=1.$ (left panel). (b)In the right panel, in a
log-linear plot we display the Chebyshev coefficients of this eigenvector.

####  4.2.2 Potential with a Coulomb type decay

Let’s consider now the Schrödinger equation (1) equipped with the potential
(8)

$q\left(x\right):=-\frac{1-5\exp\left(-2x\right)}{x}+\frac{l\left(l+1\right)}{x^{2}},\
l\in{\mathbb{R}}.$ (8)

and supplied with boundary conditions (5). Some eigenvalues computed by ChC
method when $c:=2$ and $N:=512$ are compared in Table 3 with their
counterparts computed by LGRC and perturbation methods. If for the first
eigenvalues the coincidence is excellent, the same does not happen for higher
indices. That is why in Fig. 8 in a log-linear plot we illustrate the absolute
drift in the case of the ChC method. It is observed that, for instance for the
tenth mode, the difference between the values of this eigenvalue calculated
with two very different orders of approximation N is of the order $10^{-10}$.
This leads us to believe that the eigenvalues computed with ChC are the most
accurate. It is also clear from this figure that in the second case, i.e.,
$N_{1}:=760$ and $N_{2}:=1024$ the drift oscillates less than in the first
case reported in this figure.

Table 3: The first five eigenvalues and the tenth one of the Schrödinger problem (1)-(5) when the potential has a Coulomb-type decay (8), computed by three different methods. $\;j$ |  $\;\lambda_{j}$ computed by LGRC in [8] |  $\;\lambda_{j}$ computed in [28] |  $\;\lambda_{j}$ by mapped ChC
---|---|---|---
$\;1$ | $\;\numprint{-0.0616818466333}$ | $\;-0.061681846633$ | $\;-0.06168184663316705$
$\;2$ | $\;-0.0274980999429$ | $\;-0.027498099943$ | $\;-0.02749809994382280$
$\;3$ | $\;-0.0155015616910$ | $\;-0.015501561691$ | $\;-0.01550156169420540$
$\;4$ | $\;-0.0099354968508$ | $\;-0.009935496851$ | $\;-0.009935496853885005$
$\;5$ | $\;-0.0069067013822$ | $\;-0.006906701382$ | $\;-0.006906701375461869$
$\;10$ | $\;-0.001963685230$ | $\;-0.001736111111$ | $\;-0.002059879612641054$
Figure 8: The absolute drift with respect to $N$ of the first $25$ eigenvalues
for the Schrödinger problem (1)-(5) with Coulomb-type decay potential (8),
computed by mapped ChC with scaling factor $c:=2$; red dotted line for the
case $N_{1}:=512$ and $N_{2}:=1024$ and green circled line for the case
$N_{1}:=760$ and $N_{2}:=1024.$

###  4.3 A singular Schrödinger eigenproblem on the real line

Let’s consider now the equation (1) on the real line with the potential

$q\left(x\right):=x^{2}+\frac{\nu x^{2}}{1+\mu x^{2}},$ (9)

where $\mu$ and $\nu$ are real parameters. This is a more general (anharmonic)
oscillator than the simpler harmonic one.

Two behavioral boundary conditions requiring the boundedness of the solutions
at large distance, i.e., $x\rightarrow\pm\infty$ are attached to this
equation. Actually we impose the conditions

$u\left(x\right)\rightarrow 0\;\textrm{as}\;x\rightarrow\pm\infty.$ (10)

We have to observe that these conditions are automatically satisfied in
spectral collocation based on Laguerre, Hermite and sinc functions.

SiC with $N:=500$ and scaling factor $h:=0.1$ has been used in order to
produce the following results. The fist four eigenvectors of Schrödinger
eigenproblem (1)-(10) with potential (9) are displayed in the left panel of
Fig. 9. Their coefficients are illustrated in the right panel of the same
figure. These coefficients symmetrically decrease to approximations of at
least $10^{-12}$ which means a reasonable accuracy of the method.

Figure 9: (a) Zoom in the first four eigenvectors computed by SiC for the
Schrödinger eigenproblem (1)-(10) with potential (9) (left panel). (b) In the
right panel we display the sinc coefficients for these vectors; red for the
first, green for the second, black for the third and magenta for the fourth.

We also have computed some high-lying eigenvalues namely $\lambda_{100}$,
$\lambda_{150}$ and $\lambda_{200}$, for $\nu:=1$ and $\mu:=500.$ They have
respectively the numerical values

$199.001994801512,\;301.001995781805,\;\textrm{and}\;403.001995224433.$

They are five digits approximation for the corresponding harmonic oscillator
eigenvalues $\lambda_{n}=2n-1,$ $n=1,2,\ldots.$

In [31] the author solved this problem numerically for general $\nu$ and
$\mu$. He used a very elegant variational argument of Ritz type, based on
Hermite functions, and observed that for large $\mu$ and fixed $\nu$ the
eigenvalues of this problem approximate those of the harmonic oscillator. We
have confirmed this observation even for higher index eigenvalues. In the left
panel of Fig. 10 we display the relative drift of the first $200$ eigenvalues
of the problem. It is clear that approximately the first $70$ eigenvalues are
calculated with an accuracy better than $10^{-12}$. Eigenvalues with index up
to $150$ remain at the same accuracy as $N\approx 500.$ Above this index the
accuracy decreases to $10^{-3}$. Actually, over years, a lot of literature has
gathered on this problem. Low order eigenvalues have accurately computed using
Runge-Kutta type methods for instance in [32] and [33]. In [34] the author
uses Hermite collocation in order to find only the first eigenvalues for some
Schrödinger eigenproblem. With respect to the departure from orthogonality,
i.e. the distance from zero of the scalar product of two eigenvectors, we
display in the right hand side of Fig. 10, in a log linear plot, the absolute
values of the of the scalar products of $u_{1}$ and $u_{j}$, $j=2,\ldots,200.$
It is again remarkable that this departure is less than $10^{-15}.$

It is of some importance to justify our choice for SiC. Unlike all other
methods of spectral collocation, where the differentiation matrices are highly
non-normal, in SiC these matrices are symmetric or skew-symmetric for even
respectively odd values of cut-off parameter $N$. And this is an important
numerical advantage. For instance in the MATLAB code below we use the routine
eigs instead of eig which would have been more expensive and slower.

The following very simple MATLAB code has been used.

    N=500;                                      % order of approximation
    h=0.1;                                      % spacing (scaling factor)
    [x,D]=sincdif(N,2,h); D2=D(:,:,2);          % 2nd SiC differentiation
    nu=1;mu=500;                                % parameter of the problem
    A=-D2+diag((x.^2)+nu*(x.^2)./(1+mu*(x.^2)));% the matrix
    [V,D] = eigs(A,250,0); D=diag(D);           % call MATLAB code eigs
    [t,o]=sort(real(D)); D=D(o); V=V(:,o);      % sort eigenvalues

This code is fairly similar with that from Section 4.2 with two differences.
First, the Chebyshev differentiation matrices are replaced by sinc
differentiation matrices and then the boundary conditions are absent. More
exactly, the discrete sinc functions, on which the unknown solution is
expanded, are defined by

$S_{k}\left(x,h\right):=\frac{\sin\left[\frac{\pi}{h}\left(x-x_{k}\right)\right]}{\frac{\pi}{h}\left(x-x_{k}\right)},\
k=1,2,\ldots,N,$ (11)

and satisfy the boundary conditions (10). In (11) the nodes $x_{k}$ are
equidistant with spacing $h$ and symmetric with respect to the origin.

Figure 10: (a) The relative drift of the first $200$ eigenvalues of problem
(1)-(10) with potential (9), $N_{1}:=500$ and $N_{2}:=400$-red dotted line and
$N_{1}:=300$ and $N_{2}:=400$\- green line.(left panel). (b) In the right
panel we display the orthogonality errors (deficiency) of the first
eigenvector of the same problem with respect to the subsequent $199$ of it.

##  5 Concluding remarks and open problems

Regarding the regular problems, Chebfun is unbeatable in terms of accuracy,
computation speed, and the information they provide on the accuracy of
computational process. It displays the optimal approximation order of unknowns
(eigenvectors) and how and to what extent their Chebyshev coefficients
decrease. It also specifies the degree to which some boundary conditions are
satisfied.

As for the singular problems, the situation is not so offering. In this case,
Chebfun leaves enough room for the application of usual (classical) spectral
methods. For singular eigenproblems on unbounded domains mapped ChC on the
half line and SiC, for problems on the real line, perform better than Chebfun
even for Chebfun applied on a mapped domain.

However, approaching these problems in parallel, with Chebfun as well as with
the classical spectral methods, one can get greater confidence in the accuracy
of numerical results. Regarding ChC method, using the so-called absolute or
relative drift in terms of some parameters we ensure the numerical stability
of the numerical process and can eliminate (numerically) spurious eigenvalues.
In other words we get automatically and precisely the accuracy at which a
specified set of eigenvalues is computed. The departure of orthogonality of a
set of eigenvectors provides a useful hint for the accuracy of the numerical
process. We have managed to correctly identify sets of multiple and high
indices eigenvalues (triples) and even give some numerical meaning to the
notion of continuous spectrum. This second issue obviously remains an open
one.

All in all we can say that Chebfun as well as classical spectral methods in
various forms, produce more accurate and comprehensive outcomes when solving
singular and of course regular eigenproblems than classical methods. However,
there remains an open problem, namely that of establishing (possibly
automatically) the scaling factor.

The following abbreviations are used in this manuscript:

c | scaling factor for unbounded domains
---|---
f.e. | finite element method
f. d. | finite difference method
GPS | generalized pseudospectral method
SL | Sturm-Liouville eigenproblem
ChC | Chebyshev collocation
LGRC | Laguerre Gauss Radau Collocation
SiC | sinc collocation
MATSLISE | Matlab package for SL and Schrödinger equations
MEP | Multiparameter Eigenvalue Problem
N | the order of approximation of spectral method (cutting-off parameter)
SLEDGE | SL Estimates Determined by Global Errors
SLEIGN | FORTRAN package for numerical solution to SL eigenproblem
SLDRIVER | interactive package for the previous packages

## References

  * [1] Driscoll, T. A.; Bornemann, F.; Trefethen, L.N. The CHEBOP System for Automatic Solution of Differential Equations. BIT 2008, 48 701–723
  * [2] Driscoll, T. A.; Hale, N.; Trefethen, L.N. Chebfun Guide. Pafnuty Publications, Oxford 2014
  * [3] Driscoll, T. A.; Hale, N.; Trefethen, L.N. Chebfun-numerical computing with functions. http://www.chebfun.org. Accessed 15 November 2019
  * [4] Olver, S.; Townsend, A. A fast and well-conditioned spectral method. SIAM Rev. 2013, 55, 462–489
  * [5] Trefethen, L.N.; Birkisson, A.; Driscoll, T. A. Exploring ODEs. SIAM, Philadelphia 2018
  * [6] Trefethen, L.N.: Approximation Theory and Approximation Practice, Extended Edition. SIAM, Philadelphia 2019
  * [7] Gheorghiu, C.I. Spectral Methods for Non-Standard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond. Springer-Verlag, Cham Heidelberg New-York Dondrecht London 2014
  * [8] Gheorghiu, C.I. Spectral Collocation Solutions to Problems on Unbounded Domains. Casa Cărţii de Ştiinţă Publishing House, Cluj-Napoca, Romania 2018
  * [9] Weideman, J. A. C.; Reddy, S. C., A MATLAB Differentiation Matrix Suite. ACM T. Math. Software 2000, 26, 465–519
  * [10] Roy, A.K. The generalized pseudospectral approach to the bound states of the Hulthén and the Yukawa potentials. PRAMANA-journal of physics 2005, 65, 1–15
  * [11] Shizgal, B.D., Pseudospectral Solution of the Fokker-Planck Equation with Equilibrium Bistable States: the Eigenvalue Spectrum and the Approach to Equilibrium. J. Stat. Phys. 2016. https://doi.org/10.1007/s10955-016-1594-9
  * [12] Birkhoff, G.; Rota, G-C. Ordinary Differential Equations, Fourth Edition, John Willey and Sons, New York Chichester Brisbane Toronto Singapore 1989; pp. 336–343
  * [13] Pruess, S.; Fulton, C. T. Mathematical Software for Sturm-Liouville Problem. ACM T. Math. Software 1993, 19 360–376
  * [14] Pruess, S.; Fulton, C. T.; Xie, Y. An Asymptotic Numerical Method for a Class of Singular Sturm-Liouville Problems. SIAM J. Numer. Anal. 1995, 32 1658–1676
  * [15] Marletta, M.; Pryce, J.D. LCNO Sturm-Liouville problems computational difficulties and examples. Numer. Math. 1995, 69 303–320
  * [16] Pryce, J.D.; Marletta, M. A new multi-purpose software package for Schrödinger and Sturm–Liouville computations. Comput. Phys. Comm. 1991, 62 42–54
  * [17] Bailey, P. B.; Everitt, W. N.; Zettl, A. Computing Eigenvalues of Singular Sturm-Liouville Problems. Results Math. 1991, 20 391–423
  * [18] Bailey, P.B.; Garbow, B.; Kaper, H.; Zettl, A. Algorithm 700: A FORTRAN software package for Sturm-Liouville problems. ACM T. Math. Software 1991, 17 500–501
  * [19] Ledoux, V.; Van Daele, M.; Vanden Berghe, G. MATSLISE: A MATLAB Package for the Numerical Solution of Sturm-Liouville and Schrödinger Equations. ACM T. Math. Software. 2005, 31 532–554
  * [20] Solomonoff, A.; Turkel, E. Global Properties of Pseudospectral Methods. J. Comput. Phys. 1989, 81 230–276
  * [21] Hoepffner, J. Implementation of boundary conditions. http://www.lmm.jussieu.fr/$\sim$ hoepffner/boundarycondition.pdf. Accessed 25 Aug. 2012
  * [22] Gheorghiu, C.I., Pop, I.S. A Modified Chebyshev-Tau Method for a Hydrodynamic Stability Problem. Proceedings of the International Conference on Approximation and Optimization (Romania) - ICAOR Cluj-Napoca, July 29 - August 1, Vol. II, 119–126 (1996); http://ictp.acad.ro
  * [23] Gheorghiu, C.I. On the numerical treatment of the eigenparameter dependent boundary conditions. Numer. Algor. 2018, 77 77–93
  * [24] Gheorghiu, C.I.; Hochstenbach, M.E.; PLestenjak, B.; Rommes, J. Spectral collocation solutions to multiparameter Mathieu’s system. Appl. Math. Comput. 2012, 218 11990–12000
  * [25] Plestenjak, B.; Gheorghiu, C.I.; Hochstenbach, M.E. Spectral collocation for multiparameter eigenvalue problems arising from separable boundary value problems. J. Comput. Phys. 2015, 298 585–601
  * [26] Boyd, J. P. Traps and Snares in Eigenvalue Calculations with Application to Pseudospectral Computations of Ocean Tides in a Basin Bounded by Meridians. J. Comput. Phys. 1996, 126 11–20
  * [27] Ledoux, V.; Van Daele, M.; Vanden Berghe, G. Efficient computation of high index Sturm-Liouville eigenvalues for problems in physics. Comput. Phys. Commun. 2009, 180 241–250
  * [28] Ledoux, V.; Ixaru, L.Gr.; Rizea, M.; Van Daele, M.; Vanden Berghe, G. Solution of the Schrödinger equation over an infinite integration interval by perturbation methods, revisited. Comput. Phys. Commun. 2006, 175 612–619
  * [29] Schonfelder, J.L. Chebyshev Expansions for the Error and Related Functions. Math. Comput. 1978, 32 1232–1240
  * [30] von Winckel, G. Fast Chebyshev Transform (1D).
https://www.mathworks.com/matlabcentral/fileexchange/4591-fast-chebyshev-
transform-1d. Accessed 15 May 2015

  * [31] Mitra, A.K. On the interaction of the type $\frac{\nu x^{2}}{1+\mu x^{2}}$ . J. Math. Phys. 1978, 19 2018–2022
  * [32] Simos, T.E. Some embedded modified Runge-Kutta methods for the numerical solution of some specific Schrödinger equations. J. Math. Chem. 1998, 24 23–37
  * [33] Simos, T.E. An accurate finite difference method for the numerical solution of the Schrödinger equation. J. Comput Appl. Math. 1998, 91 47–61
  * [34] Trif, D. Matlab package for the Schrödinger equation. J. Math. Chem. 2008, 43 1163–1176