# Self-supervision meets kernel graph neural models:
From architecture to augmentations

Jiawang Dan21, Ruofan Wu21, Yunpeng Liu23, Baokun Wang2 Changhua Meng2,
Tengfei Liu2, Tianyi Zhang2, Ningtao Wang2, Xing Fu2, Qi Li3, Weiqiang Wang2 1
Equal contribution 2Ant Group  3Tsinghua University  {yancong.djw, ruofan.wrf,
lyp402072, yike.wbk, aaron.ltf<EMAIL_ADDRESS><EMAIL_ADDRESS>{changhua.mch, zty113091, zicai.fx,
<EMAIL_ADDRESS>
###### Abstract

Graph representation learning has now become the de facto standard when
handling graph-structured data, with the framework of message-passing graph
neural networks (MPNN) being the most prevailing algorithmic tool. Despite its
popularity, the family of MPNNs suffers from several drawbacks such as
transparency and expressivity. Recently, the idea of designing neural models
on graphs using the theory of graph kernels has emerged as a more transparent
as well as sometimes more expressive alternative to MPNNs known as kernel
graph neural networks (KGNNs). Developments on KGNNs are currently a nascent
field of research, leaving several challenges from algorithmic design and
adaptation to other learning paradigms such as self-supervised learning. In
this paper, we improve the design and learning of KGNNs. Firstly, we extend
the algorithmic formulation of KGNNs by allowing a more flexible graph-level
similarity definition that encompasses former proposals like random walk graph
kernel, as well as providing a smoother optimization objective that alleviates
the need of introducing combinatorial learning procedures. Secondly, we
enhance KGNNs through the lens of self-supervision via developing a novel
structure-preserving graph data augmentation method called _latent graph
augmentation (LGA)_. Finally, we perform extensive empirical evaluations to
demonstrate the efficacy of our proposed mechanisms. Experimental results over
benchmark datasets suggest that our proposed model achieves competitive
performance that is comparable to or sometimes outperforming state-of-the-art
graph representation learning frameworks with or without self-supervision on
graph classification tasks. Comparisons against other previously established
graph data augmentation methods verify that the proposed LGA augmentation
scheme captures better semantics of graph-level invariance.

##  I Introduction

Recent years have witnessed surging developments in graph representation
learning (GRL) which utilizes neural models over graph-structured data for
downstream tasks like node classification [18] and graph classification [40].
Among these the task of graph classification has shown promising results in
applications like drug discovery [34] and protein function prediction [13].
The design of GRL algorithms has attracted significant interest, with the idea
of message-passing graph neural networks [12] (hereafter abbreviated as MPNNs)
being the de facto choice. In its most standard form, MPNNs obtain node
representations via aggregating neighborhood information for each node in a
recursive manner, with an optional combination mechanism at each step. The
node representations are further aggregated into a graph-level representation
for downstream tasks like graph classification [40]. In the seminal work [40],
the authors explored the expressivity limits of message-passing protocols and
concluded that the expressivity limit of such kinds of procedures is bounded
by first-order Weisfeler-Leman isomorphism tests. While there exist more
theoretically expressive variants such as higher-order graph neural networks
[23], these alternatives are typically computationally expensive and lose the
scalability of MPNNs.

In a recent line of works [24, 8, 25], the authors have studied an alternative
way of defining neural models over graphs which are inspired by the theory of
graph kernels [30], we will call such variants kernel graph neural networks
(KGNNs). KGNNs use a set of trainable hidden graphs and use the similarity of
the input graph with the hidden graphs as the building block of graph
representations. KGNNs provide a powerful alternative to MPNNs in the
following aspects: Firstly, KGNNs allow more transparent graph modeling, since
the obtained hidden graphs themselves are easily visualized, they provide
interpretations of the learned model. Secondly, KGNNs share several important
properties with MPNNs such as permutation invariance and scalability. Finally,
it was proved in previous works [8] that KGNNs are provably more expressive
than MPNNs. Despite these advantages, direct optimization over hidden graphs
is computationally intractable. In the original proposal [24], the authors
represent hidden graphs as continuous-valued matrices, resulting in subtle
definitions. It is thus of interest to further explore this modeling paradigm
via reinspecting the similarity definition, and its relation to the theory of
graph kernels.

Another recent trend in the area of GRL is the development of self-supervision
techniques [39], as they are particularly effective in label-scarce scenarios,
providing a principled way of facilitating unlabelled data. In self-
supervision, one typically constructs a pretext task whose supervisory signals
are derived through a certain sense of _invariance_ , and utilizes some
invariance-promoting objectives to obtain pre-trained models such as
contrastive objectives [43] and non-contrastive objectives [37]. While there
has been abundant work on self-supervision for MPNNs [39], so far as we have
noticed there have been no studies on self-supervision over KGNNs. Applying
self-supervision to the (pre)training of KGNN is challenging since KGNNs do
not produce _intermediate node-level representations_ which are necessary for
many self-supervised GRL frameworks [35, 36, 16]. To perform (invariance-
based) self-supervised learning with KGNNs, we need to construct our
supervision signals via graph-level invariant transforms instead of node-level
operations. Although common practices in graph data augmentations [43] also
apply to KGNNs, they might not be able to express the right semantics of
invariance that could be captured by KGNN architectures: As KGNNs are
intuitively understood to focus more on graph structures, we hypothesize that
successful pretraining of KGNN networks require _structure-preserving_ graph
data augmentations, which has yet been under-explored in the current
literature.

In this paper, we make several contributions to the design and learning of
KGNNs which are elaborated on as follows:

  * •

Our first contribution is a novel extension of the current KGNN design that
allows a more flexible definition than contemporary KGNNs which are mostly
derived through random walk graph kernels. We also instantiate a concrete
model that utilizes the technique of graph diffusion [11]. The proposed model
is intuitively understood as a smoothed version of contemporary KGNNs.

  * •

Our second contribution is a novel structure-preserving graph data
augmentation method that serves as a building block for self-supervised pre-
training using KGNNs. The semantics of structure-preserving is expressed
through a wide range of random graph models, as well as being adaptive to
individual instances.

  * •

Finally, we conduct extensive empirical evaluations under the task of graph
classification over benchmark datasets. Experimental results suggest our
proposed model improves over contemporary KGNNs, and sometimes achieves
performance on par with or even outperforms state-of-the-art methods.

##  II Methodology

Figure 1: A concise illustration of the proposed self-supervised graph
representation learning framework that consists of two main components: The
first one is the SWAG architecture that utilizes graph diffusion transform to
allow continuous optimization. The second one is the LGA method for graph data
augmentation that reflects the notion of structural invariance which is
beneficial for the learning of KGNNs.

###  II-A Kernel graph neural networks (KGNNs)

Let $G=(V,E)$ be an undirected graph associated with node features
$X\in\mathbb{R}^{d}$. We will use $V_{G}$ and $E_{G}$ to denote the node set
and edge set of graph $G$, further denote $[N]$ as the set $\\{1,\ldots,N\\}$.
Given a training set $\mathbf{G}=\\{G_{1},\ldots,G_{N}\\}$ of $N$ graphs that
are of varying sizes with each one having a label $Y_{i},i\in[N]$, we are
interested in learning a graph representation extractor $h$ that maps graph
objects to euclidean vectors, as well as learning a task-specific downstream
model $g$ (sometimes referred to as _predictor_) that maps graph
representations to predictions. In supervised learning paradigms, $h$ and $g$
are trained in an end-to-end fashion, with $h$ usually being certain kinds of
graph neural networks and $g$ being an MLP. In self-supervised learning
paradigms, $h$ and $g$ are trained separately, with $h$ being trained using
supervision signals that are not directly associated with the downstream task
label. After obtaining the pre-trained encoder, we either train $g$ using
probing (i.e., freezing the pre-trained graph representations) or using fine-
tuning techniques. In this paper, we will focus on self-supervised approaches
with $f$ being parameterized using KGNNs. Specifically, a KGNN model typically
involves a collection of $M$ _learnable hidden graphs_ as trainable
parameters, i.e., $\mathcal{H}=\\{H_{m}\\}_{m\in[M]}$. Graph representations
are obtained via computing some similarity metric between the input graph and
all the hidden graphs, with the prevailing practice being inspired by the
theory of graph kernels [30]. With the most representative architecture
derived from the theory of _random walk kernels_ [24, 8]. Given input graphs
$G$ and $G^{\prime}$, the random walk kernel evaluation $K(G,G^{\prime})$
computes a weighted sum of the counts of random walks of varying lengths that
the two graphs have in common. The kernel definition allows incorporating node
features [24] via treating pairwise feature similarity as additional weights.
Formally, given a walk length $p\geq 1$, the $p$-th constituent of random walk
kernel between $G$ and $G^{\prime}$ is computed as

$\displaystyle K^{(p)}(G,G^{\prime})=\sum_{k,l\in V_{G^{\prime}}}\sum_{i,j\in
V_{G}}\langle x_{i},x^{\prime}_{k}\rangle
A^{p}_{ij}{A^{\prime}}^{p}_{kl}\langle x_{j},x^{\prime}_{l}\rangle,$ (1)

where we define $A$ and $A^{\prime}$ as the adjacency matrix of $G$ and
$G^{\prime}$, and $x^{\prime}_{u},u\in V_{G^{\prime}}$ as the node feature of
graph $G^{\prime}$. To evaluate graph kernels, we may compute a weighted sum
(ideally infinite) over all $K^{(p)}(G,G^{\prime}),p>1$. In practice, it was
observed that computing $p$ up to a small to moderate integer (i.e., $p\leq
3$) suffices for most applications, and the representations are constructed
via concatenating all $K^{(p)}(G,G^{\prime})$ with $1\leq p\leq P$ and
$G^{\prime}\in\mathcal{H}$.

###  II-B Extensions and smoothed random walks

According to the original definition (1), the trainable parameters involve a
collection of (hidden) feature matrices and a collection of (hidden) adjacency
matrices. However, adjacency matrices are binary-valued and thus do not fit
into the standard gradient-based learning framework in neural models. Besides,
directly optimizing potentially large binary matrices will incur a
combinatorial level of complexity and is computationally intractable. To fix
this issue, in [24] the authors suggested replacing $A^{\prime}$ in (1) by
$\text{ReLU}(W)$ with $W$ being an arbitrary continuous-valued matrix. While
this fix addresses the computational issue, it breaks the semantics of graph
kernels (as the resulting formula does not count the common random walks).
However, since exact kernel evaluations may not contribute significantly to
the performance in downstream tasks, we might directly extend the definition
of (1) so that it allows smoother optimization formulations, in particular, we
propose the following extension

(2)

where $s_{K}$ denotes some similarity function that applies to euclidean
vectors, $B$ and $B^{\prime}$ are some _continuous-valued characteristics_ of
graphs $G$ and $G^{\prime}$, respectively. The above extension naturally
serves as a concrete similarity measure, while at the same time allowing
direct gradient-based optimization. In this paper, we will stick to the
sleekest similarity formulation of $s_{K}$ which is the inner product
function. What remains is to find a proper characteristic $B$, which in this
paper we choose to be the diffusion matrix of graph $G$ [11].

$\displaystyle\mathcal{D}_{\beta}(G)=\sum_{j=0}^{\infty}\beta_{j}T^{j},$ (3)

where $\beta_{j}$ stands for some coefficient configurations that make the
summation (3) convergent and $T$ is a transition matrix defined as
$D^{-1/2}AD^{-1/2}$, with $D$ being the diagonal matrix with diagonal values
corresponding to node degrees. For $B^{\prime}$ in (2), we may simply
parameterize it as a continuous-valued symmetric matrix with values in
$[0,1]$. We term the resulting encoder as Smoothed random WAlk Graph neural
network (SWAG), whose encoding function with maximum walk length $P>1$ and the
number of hidden graphs $M$ is given by:

$\displaystyle
h_{\text{SWAG}}(G)=\operatorname*{\mathchoice{\Big{\|}}{\big{\|}}{\|}{\|}}_{m\in[M],p\in[P]}\widetilde{K}^{(p)}(G,H_{m}),$
(4)

where we use $\operatorname*{\mathchoice{\Big{\|}}{\big{\|}}{\|}{\|}}$ to
denote the concatenation operation. We use $d_{\text{o}}=MP$ to denote the
dimension of output graph representation.

###  II-C Self-supervision via struture-preserving graph data augmentation

In this paper, we are mainly interested in self-supervised learning primitives
using _augmentations_ , with two canonical frameworks being contrastive
learning and non-contrastive learning. Both frameworks are based on an
augmentation mechanism. Under the context of GRL, for any input graph $G$
(sometimes referred to as _anchor graph_), we augment it with a positive
sample $G^{+}\sim\mathcal{A}(\cdot|G)$, where we use $\mathcal{A}(\cdot|G)$ to
denote the (conditional) augmentation distribution given the anchor graph $G$
that is supported on the universe of graphs $\mathcal{G}$.

Constrastive pre-training In contrastive frameworks, the learning objective is
derived via contrasting augmented views (positive samples) with a set of
negative samples $\\{G^{-}_{1},\ldots,G^{-}_{K}\\}$, each being sampled from
the (conditional) distribution $\mathcal{N}(\cdot|G)$. We will adopt the
InfoNCE objective [26]:

(5)

with $s_{L}$ being some similarity function. We will be using the form
$s_{L}(x,y)=\langle\phi(x),\phi(y)\rangle$ with $\phi$ being a two-layer MLP.

Non-contrastive pre-training Non-contrastive frameworks alleviate the need for
negative sampling by directly maximizing the similarity between input graphs
and their augmentation. We will utilize the simplest objective [4]

$\displaystyle\mathcal{L}_{\text{non-
cst}}=\frac{1}{N}\sum_{G\in\mathbf{G}}\langle\psi(h(G)),\text{sg}(\psi(h(G^{+}))\rangle,$
(6)

where $\psi$ is a _prediction head_ that takes its form as a two-layer MLP and
sg stands for the stop-gradient operation. Both components have been shown
necessary for avoiding collapsed solutions [5].
The prevailing paradigm in graph data augmentations are of the _perturbation_
style, i.e., perturbing node attributes or perturbing the graph structure [6].
It has been observed in previous works [43] that the augmentation quality in
self-supervised graph representation learning plays a vital role in downstream
performance, most likely due to the different types of _invariance_ semantics
encoded by the augmentation strategy. In particular, the standard practice of
graph structure perturbation techniques that adds or drops either nodes or
edges essentially expresses a strong belief in invariance with respect to a
small perturbation measured in terms of graph edit distance [10]. However,
this may not provide reasonable invariance if the underlying graph is well-
structured. For example, in the case of stochastic block model [14], randomly
dropping or adding a significant proportion of edges changes the block
structure and hence makes the precise semantics of the invariance
questionable. Moreover, the inductive bias of the encoder also has an
intriguing relationship between the augmentation strategy as well as
downstream performances of self-supervised pre-training [29]. Under the
context of KGNNs, we intuitively expect that the learned hidden graphs shall
somehow extract the _signals_ in the underlying graph instead of noise.
Therefore, we require the graph data augmentation strategy to be _structure-
preserving_. To begin with, we will use the following random graph model [14]
as a starting point that describes the generating process of the underlying
graph $G$, with each entry of its adjacency matrix sampled independently from
a Bernoulli distribution,

$\displaystyle
A_{uv}\overset{\text{i.i.d}}{\sim}\text{Ber}(\theta_{uv}),\forall u,v\in
V(G).$ (7)

The above formulation is nonparametric and very general, and covers many
important generative mechanisms in random graph models, to name a few:

Graphon [9]

Suppose that each node $v\in V(G)$ is associated with a scalar
$\xi_{u}\overset{\text{i.i.d}}{\sim}\text{Unif}(0,1)$ that is independently
drawn from a uniform distribution. And the generating probability is
constructed via the _graphon_ $f:[0,1]^{2}\mapsto[0,1]$, with
$\theta_{uv}=f(\xi_{u},\xi_{v})$.

Latent space model [14]

We may alternatively let the unobserved node-wise characteristic be an
arbitrary $d$-dimensional vector $\\{\lambda_{v}\\}_{v\in V(G)}$, and the
generating probability is formulated as
$\theta_{uv}=f(\lambda_{u},\lambda_{v})$, where with a slight abuse of
notation we let $f$ be some continuous function that takes two $d$-dimensional
vector arguments.

A notable result in modern random graph theory is that while precisely
estimating the function $f$ in graphon or latent space models may encounter
identifiability issues, the generating probabilities
$\Theta=\\{\theta_{uv}\\}_{u,v\in V(G)}$ could be recovered in a provable
fashion using the idea of _universal singular value thresholding_ (USVT) [3].
Inspired by this phenomenon, we use the estimated generating matrix
$\widehat{\Theta}$ as a reference of _graph structure_ , and use this matrix
to generate augmented views, which share (in an asymptotic fashion) similar
graph structures with the input graph $G$. Since the process applies to many
latent graph generative models, we term our augmentation method Latent Graph
Augmentation (LGA), detailed in algorithm 1.

Algorithm 1 LGA using universal singular value thresholding (USVT)

1:Input graph $G=(V,E)$ with adjacency matrix $A$, threshold $\tau>0$.

2:Let $m$ be the number of nodes in $G$. Obtain the singular value
decomposition of $A$: $A=\sum_{i}^{m}\sigma_{i}u_{i}v_{i}^{T}$

3:Select indices based on the threshold
$\mathcal{S}=\\{i:\sigma_{i}\geq\tau\sqrt{m}\\}$.

4:Compute $\widetilde{\Theta}=\sum_{i\in\mathcal{S}}\sigma_{i}u_{i}v_{i}^{T}$

5:Obtain the recovered random graph probability matrix $\widehat{\Theta}$ via
clipping $\widetilde{\Theta}$:
$\displaystyle\widehat{\Theta}_{uv}=\begin{cases}\widetilde{\Theta}_{uv}\quad&\text{If
$\widetilde{\Theta}_{uv}\in[0,1]$ }\\\ 0&\text{If $\widetilde{\Theta}_{uv}<0$
}\\\ 1&\text{If $\widetilde{\Theta}_{uv}>1$ }\end{cases}$ (8)

6:Generate matrix $A^{+}$ with each
$A^{+}_{uv}\overset{\text{i.i.d}}{\sim}\text{Ber}(\widehat{\Theta}_{uv})$ with
$u,v\in[m]$.

7:return Positive sample $G^{+}$ derived from adjacency matrix $A^{+}$, i.e.,
with node features left untouched.

The most important step in LGA is to select a sufficiently large singular
value based on a predefined (relative) threshold $\tau$, with its magnitude
reflecting the trade-off between signal extraction and noise removal. In this
paper, we cast the threshold $\tau$ as a tunable _hyperparameter_ which is
regarded as a crude estimate of the upper bound of the matrix $A-\Theta$
(normalized by the square root of graph size) that is shown to theoretically
characterize the recovery under dense graph regimes [3].

###  II-D Complexity considerations

Now we present a brief analysis of the computational complexity of the
proposed approach. Firstly, the complexity of obtaining the graph
representation could be analyzed in a similar fashion as in [24]: Ff we use
$M$ hidden graphs with each involved in computing similarities up to order $P$
will incur a computational cost of $O(Md(nm+Pm^{2}+Pn^{2}))$ for a dense
sample graph and $O(Md(nm+Pm^{2}+Ps))$ for a sparse sample graph with $s$
edges. The resulting complexity is comparable to a $P$-layer MPNN model with
each layer producing hidden representations of dimension $m$. Next, we
consider the cost of obtaining the characteristic matrix $B$ using graph
diffusion. While we may use a truncated matrix sum with a moderate number of
steps to reasonably approximate the diffusion matrix, we found in our
experiments that typically a few diffusion steps (i.e., no more than three
steps) suffice for performance improvements. Consequently, the resulting
computational cost is dominated by that of representation computation.
Finally, the computation cost of LGA (dominated by that of SVD factorization
with $O(n^{3})$) is only incurred once for each sample graph during training
time. Under moderate-sized graphs, we argue that the computational cost is
controllable.

##  III Related works

###  III-A Neural models for graph representation learning

As graphs are generally believed to contain rich structural knowledge, neural
architectures over graph-structured objects require different inductive biases
compared to neural architectures over independently identically distributed
(i.i.d.) data. The authors in [1] proposed the concept of _relational
inductive bias_ encoded by MPNNs that are closely related to dynamic
programming. However, MPNNs are criticized for their limited expressivity [23,
40], we refer the readers to the article [23] for a thorough overview. A more
recent line of works derives neural architectures over graphs from the theory
of graph kernels [20, 24, 25]. In comparison to MPNNs, kernel approaches
compute convolutions over graphs in a different way that might be regarded as
a more transparent generalization of the convolution operation used in image
modeling [8]. Besides, kernel approaches allow intuitive visualization of the
learned hidden graphs (graph filters), thereby providing a natural sense of
_interpretability_.

###  III-B Self-supervision on graphs

Early developments on self-supervised learning (SSL) on graphs considered
using link prediction as the pretext task for downstream applications [19].
Later developments have been mostly focusing on applying SSL frameworks based
on _graph data augmentations (GDA)_ [38, 43, 37] , with the prevailing
practices mostly involving editing graph structure (i.e., adding or dropping
nodes or edges) and node features [6] either randomly [43] or non-randomly
[36, 21, 22]. These modifications are based on the belief that the invariance
relation is captured through certain edits of graphs.

##  IV Experiments

In this section, we report empirical evaluations on benchmark datasets. The
evaluations are conducted for both supervised learning tasks which focus on
the empirical performance of the proposed SWAG network design, and self-
supervised learning tasks which further investigate the effectiveness of the
LGA mechanism. We also provide visualizations of the learned hidden graphs as
an indirect way of comparing canonical patterns in the data, as well as an
ablation study that inspects the effects of important hyperparameters.

###  IV-A Datasets

We use $8$ public datasets [41], with $4$ bio/chemo-informatics network
datasets: MUTAG, D&D, NCI1 and PROTEINS, and $4$ social interaction network
datasets: IMDB-BINARY(IMDB-B), COLLAB, REDDIT-BINARY(RDT-B) and REDDIT-
MULTI-5K(RDT-M5K). Across all the datasets, the evaluation task is graph
classification with varying numbers of classes. For a more detailed
introduction of the datasets as well their summary statistics, we refer to
[41] for details.

###  IV-B Baseline comparisons

As our experiments involve both supervised learning and self-supervised
learning, we describe our baseline comparison strategy corresponding to both
learning schemes as follows:
Supervised learning baselines We pick three representative graph kernel
methods: shortest path kernel (SP) [2], graphlet kernel [32] and Weisfeiler-
leman subtree kernel [31]. For MPNN methods, choose five influential MPNNs
that applies to graph classification problems, namely DGCNN [44], DiffPool
[42], ECC [33], GIN [40] and GraphSAGE [15]. Additionally, we compare our
proposed model with two recent KGNN models, RWGNN [24] and GRWNN [25].
Self-supervised learning baselines For a thorough investigation, we make two
types of comparisons:

Comparison with MPNN encoders

We will compare our proposed SWAG encoder with LGA augmentation scheme with
various state-of-the-art self-supervised GRL frameworks: GraphCL [43],
InfoGraph [35], GCC [28], AD-GCL [36], RGCL [21], GraphMAE [16] and GCL-SPAN
[22], with their underlying encoder being a five-layer GIN [40]. We note that
such kind of comparisons are rarely done in previous works on KGNNs.

Comparison with KGNN encoders

Among the aforementioned baseline methods, four of them indeed apply to KGNN
encoders. Including GraphCL [43], GCC [28], RGCL [21] and GCL-SPAN [22]. We
additionally apply these augmentation strategies with a standard RWGNN encoder
as a comparison of self-supervised KGNN models.

TABLE I: Experimental results on supervised graph classification over $4$ molecular network datasets and $4$ social network datasets, reported with format ${\text{{mean}}}_{\pm\text{{std}}}$, with mean and std (abbreviation for standard deviation) computed under $10$ trials for each setting. We use bolded red to highlight the best performance and underlined blue to highlight the second best performance both in the sense of mean value. | MUTAG | D&D | NCI1 | PROTEINS | IMDB-B | COLLAB | RDT-B | RDT-M5K
---|---|---|---|---|---|---|---|---
SP | ${80.2}_{\pm 6.5}$ | ${78.1}_{\pm 4.1}$ | ${72.7}_{\pm 1.4}$ | ${75.3}_{\pm 3.8}$ | ${57.7}_{\pm 4.1}$ | ${79.9}_{\pm 2.7}$ | ${89.0}_{\pm 1.0}$ | ${51.1}_{\pm 2.2}$
GR | ${80.8}_{\pm 6.4}$ | ${75.4}_{\pm 3.4}$ | ${61.8}_{\pm 1.7}$ | ${71.6}_{\pm 3.1}$ | ${63.3}_{\pm 2.7}$ | ${71.1}_{\pm 1.4}$ | ${76.6}_{\pm 3.3}$ | ${38.1}_{\pm 2.3}$
WL | ${84.6}_{\pm 8.3}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{78.1}_{\pm 2.4}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{84.8}_{\pm 2.5}}$ | ${73.8}_{\pm 4.4}$ | ${72.8}_{\pm 4.5}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{78.0}_{\pm 2.0}}$ | ${74.9}_{\pm 1.8}$ | ${49.6}_{\pm 2.0}$
DGCNN | ${84.0}_{\pm 6.7}$ | ${76.6}_{\pm 4.3}$ | ${76.4}_{\pm 1.7}$ | ${72.9}_{\pm 3.5}$ | ${69.2}_{\pm 3.0}$ | ${71.2}_{\pm 1.9}$ | ${87.8}_{\pm 2.5}$ | ${49.2}_{\pm 1.2}$
DiffPool | ${79.8}_{\pm 7.1}$ | ${75.0}_{\pm 3.5}$ | ${76.9}_{\pm 1.9}$ | ${73.7}_{\pm 3.5}$ | ${68.4}_{\pm 3.3}$ | ${68.9}_{\pm 2.0}$ | ${89.1}_{\pm 1.6}$ | ${53.8}_{\pm 1.4}$
GIN | ${84.7}_{\pm 6.7}$ | ${75.3}_{\pm 2.9}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{80.0}_{\pm 1.4}}$ | ${73.3}_{\pm 4.0}$ | ${71.2}_{\pm 3.9}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{75.6}_{\pm 2.3}}$ | ${89.9}_{\pm 1.9}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{56.1}_{\pm 1.7}}$
GraphSAGE | ${83.6}_{\pm 9.6}$ | ${72.9}_{\pm 2.0}$ | ${76.0}_{\pm 1.8}$ | ${73.0}_{\pm 4.5}$ | ${68.8}_{\pm 4.5}$ | ${73.9}_{\pm 1.7}$ | ${84.3}_{\pm 1.9}$ | ${50.0}_{\pm 1.3}$
RWGNN | ${89.2}_{\pm 4.3}$ | ${77.6}_{\pm 4.7}$ | ${73.9}_{\pm 1.3}$ | ${74.7}_{\pm 3.3}$ | ${70.6}_{\pm 4.4}$ | ${71.9}_{\pm 2.5}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{90.4}_{\pm 1.9}}$ | ${53.4}_{\pm 1.6}$
GRWNN | ${83.4}_{\pm 5.6}$ | ${75.6}_{\pm 4.6}$ | ${67.7}_{\pm 2.2}$ | ${74.9}_{\pm 3.5}$ | ${72.8}_{\pm 4.2}$ | ${72.1}_{\pm 1.9}$ | ${90.0}_{\pm 1.8}$ | ${54,4}_{\pm 1.7}$
SWAG(SL) | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{89.4}_{\pm 6.0}}$ | ${77.8}_{\pm 3.8}$ | ${73.6}_{\pm 1.8}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{76.3}_{\pm 4.4}}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{73.1}_{\pm 3.3}}$ | ${74.1}_{\pm 1.6}$ | ${90.3}_{\pm 1.8}$ | ${55.4}_{\pm 2.2}$
SWAG(SSL+FT) | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{90.3}_{\pm 6.2}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{79.5}_{\pm 2.3}}$ | ${76.2}_{\pm 2.7}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{77.5}_{\pm 4.5}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{74.4}_{\pm 3.3}}$ | ${74.6}_{\pm 2.5}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{91.5}_{\pm 1.7}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{56.2}_{\pm 1.7}}$

###  IV-C Experimental setup

Evaluation protocol We perform $10$-fold cross-validation to obtain testing
performance across all the models, with the same splitting configuration as in
[7]. Within each fold, we use $10\%$ of the data as a hold-out validation set
that enables model selection. We use accuracy as the evaluation metric across
all the datasets.
Training configurations Across all the experimental trials, we use a SWAG
network architecture as the representation extractor $h$, and a two-layer MLP
with hidden dimension $32$ as the downstream predictor. As not all of the
chosen datasets contain node features, for those without node features we use
node degree as the node feature. For the graph diffusion transform, we use the
personalized PageRank (PPR) configuration $\beta_{j}=\alpha(1-\alpha)^{j}$
with teleport probability $\alpha=0.15$ as suggested in [11]. We train for
$500$ epochs using the Adam optimizer [17] with a learning rate $0.01$, under
a batch size of $64$. Across all the experimental trials involving KGNN, we
tune the number of hidden graphs $M$ in the range $\\{8,16\\}$, with the
number of nodes in each hidden graph tuned in the range $\\{5,10\\}$. For the
node features corresponding to hidden graphs, we tune the feature dimension in
the range $\\{32,64\\}$ for chemo-informatic network datasets, and $\\{4,8\\}$
for the social network datasets. We choose the walk length, as well as the
number of graph diffusion steps to be no greater than $3$ for computational
efficiency. During self-supervision, we tune the thresholding parameter $\tau$
in the USVT scheme inside the range $[0.3,4.2]$, with a detailed study on this
parameter elaborated in section IV-G. We adopt both contrastive pre-training
objective (5) with negative samples chosen as the remaining ones in the same
batch for any anchor graph, and the non-contrastive pretraining objectives
(6). For both objectives, the additionally involved neural functions $\phi$
and $\psi$, as explained in section II-C, are parameterized by a two-layer MLP
with $32$ hidden units. During downstream task adaptation, we tested both
probing and fine-tuning and report the better-performing one.
Implementations We mostly base our model implementation on PyTorch [27]. For
all the baseline methods, we adopt open-source implementations available in
the corresponding papers and inherit their default training configurations.

###  IV-D Performance on supervised graph classification

TABLE II: Experimental results on self-supervised graph classification over $4$ molecular network datasets and $4$ social network datasets, reported with format ${\text{{mean}}}_{\pm\text{{std}}}$, with mean and std (abbreviation for standard deviation) computed under $10$ trials for each setting. We use bolded red to highlight the best performance and underlined blue to highlight the second best performance both in the sense of mean value. | MUTAG | D&D | NCI1 | PROTEINS | IMDB-B | COLLAB | RDT-B | RDT-M5K
---|---|---|---|---|---|---|---|---
With MPNN encoder
GraphCL | ${81.5}_{\pm 7.0}$ | ${70.4}_{\pm 3.8}$ | ${72.4}_{\pm 3.0}$ | ${69.3}_{\pm 3.6}$ | ${64.7}_{\pm 5.1}$ | ${74.4}_{\pm 1.3}$ | ${83.4}_{\pm 3.0}$ | ${49.9}_{\pm 1.3}$
GCC | ${83.6}_{\pm 7.4}$ | ${67.9}_{\pm 2.3}$ | ${68.0}_{\pm 2.0}$ | ${67.2}_{\pm 2.4}$ | ${74.4}_{\pm 4.0}$ | ${77.5}_{\pm 1.7}$ | ${85.5}_{\pm 3.1}$ | ${48.7}_{\pm 1.2}$
RGCL | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{89.9}_{\pm 6.2}}$ | ${78.8}_{\pm 3.2}$ | ${78.4}_{\pm 2.1}$ | ${75.5}_{\pm 3.8}$ | ${73.2}_{\pm 4.0}$ | ${70.8}_{\pm 1.7}$ | ${89.5}_{\pm 2.5}$ | ${56.2}_{\pm 2.1}$
GCL-SPAN | ${85.7}_{\pm 7.9}$ | ${79.0}_{\pm 3.0}$ | ${73.8}_{\pm 2.0}$ | ${74.2}_{\pm 4.4}$ | ${71.5}_{\pm 3.1}$ | ${71.3}_{\pm 2.0}$ | ${71.9}_{\pm 3.0}$ | ${53.7}_{\pm 1.7}$
InfoGraph | ${87.7}_{\pm 6.7}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{79.5}_{\pm 2.1}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{80.3}_{\pm 1.8}}$ | ${74.1}_{\pm 3.6}$ | ${71.4}_{\pm 2.9}$ | ${75.4}_{\pm 2.1}$ | ${91.1}_{\pm 1.5}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{56.5}_{\pm 1.2}}$
AD-GCL | ${88.1}_{\pm 1.4}$ | ${74.8}_{\pm 0.7}$ | ${69.4}_{\pm 0.6}$ | ${73.6}_{\pm 0.6}$ | ${71.2}_{\pm 0.5}$ | ${73.3}_{\pm 0.5}$ | ${84.9}_{\pm 1.4}$ | ${54.9}_{\pm 0.7}$
GraphMAE | ${85.6}_{\pm 1.4}$ | ${72.9}_{\pm 1.6}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{79.3}_{\pm 0.5}}$ | ${75.7}_{\pm 0.3}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{75.4}_{\pm 0.6}}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{80.2}_{\pm 0.2}}$ | ${84.5}_{\pm 0.5}$ | ${52.9}_{\pm 0.3}$
With RWGNN encoder
GraphCL | ${86.2}_{\pm 6.3}$ | ${75.1}_{\pm 3.0}$ | ${67.7}_{\pm 2.8}$ | ${71.0}_{\pm 4.5}$ | ${67.7}_{\pm 4.4}$ | ${73.1}_{\pm 1.8}$ | ${82.1}_{\pm 4.9}$ | ${52.9}_{\pm 2.1}$
GCC | ${84.7}_{\pm 6.7}$ | ${74.9}_{\pm 4.8}$ | ${67.0}_{\pm 1.9}$ | ${72.1}_{\pm 3.0}$ | ${71.1}_{\pm 4.7}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{78.3}_{\pm 1.6}}$ | ${85.1}_{\pm 2.8}$ | ${49.0}_{\pm 1.1}$
RGCL | ${83.5}_{\pm 8.3}$ | ${74.8}_{\pm 3.8}$ | ${64.2}_{\pm 2.1}$ | ${72.7}_{\pm 2.6}$ | ${70.9}_{\pm 3.9}$ | ${68.2}_{\pm 2.4}$ | ${82.9}_{\pm 2.7}$ | ${49.5}_{\pm 1.9}$
GCL-SPAN | ${85.2}_{\pm 8.5}$ | ${76.2}_{\pm 2.8}$ | ${62.6}_{\pm 2.7}$ | ${75.2}_{\pm 4.5}$ | ${71.4}_{\pm 2.7}$ | ${68.4}_{\pm 1.9}$ | ${79.4}_{\pm 1.9}$ | ${49.0}_{\pm 2.3}$
SWAG+LGA(infonce) | ${89.1}_{\pm 6.7}$ | ${78.0}_{\pm 2.4}$ | ${75.6}_{\pm 3.6}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{76.3}_{\pm 4.3}}$ | ${74.0}_{\pm 2.9}$ | ${74.6}_{\pm 2.5}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{91.5}_{\pm 2.5}}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{56.2}_{\pm 1.7}}$
SWAG+LGA(simsiam) | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{90.3}_{\pm 6.2}}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{79.5}_{\pm 2.3}}$ | ${76.2}_{\pm 2.7}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{77.5}_{\pm 4.5}}$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\underline{74.4}_{\pm 3.3}}$ | ${74.1}_{\pm 2.0}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{91.5}_{\pm 1.7}}$ | ${56.0}_{\pm 2.1}$

We present experimental evaluations on supervised graph classification in
table I. We have the following observations:

  * •

Our proposed SWAG model is demonstrated to achieve comparable or superior
performance to the two KGNN baselines, measured in _mean_ performance. On
datasets with more dense graphs present such as IMDB-B and COLLAB, we found
the improvements to be more evident, which might be attributed to the
topological information cascading provided by the diffusion operation in the
SWAG mechanism.

  * •

When compared against MPNN and graph kernel baselines, SWAG achieves
comparable performance with SOTA models on $6$ of the $8$ datasets, with the
results in MUTAG, PROTEINS, and IMDB-B datasets relatively convincing. Another
notable fact is that for certain datasets like NCI$1$ and COLLAB, using WL
kernels leads to better results than neural GRL primitives, which is
understandable since graph-level discrimination tasks are more sensitive to
the expressivity limit of MPNNs, which was explicitly achieved using WL
kernels.

  * •

We additionally report the performance of SWAG under self-supervision with LGA
augmentation, the result shows that for all the $8$ datasets, self-supervision
is beneficial to model performance.

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 2: Visualization of $4$ randomly chosen learned subgraphs. Figure 2(a)
through figure 2(d) are visualizations of hidden graphs learned from the
PROTEINS dataset; Figure 2(e) through figure 2(h) are learned from the REDDIT-
BINARY(RDT-B) dataset.

###  IV-E Performance on self-supervised graph classification

We present a detailed report comparing our approach and self-supervised GRL
baselines in table II, where we report separately the results under two
different objectives. We summarize our findings as follows:

  * •

When compared against baselines with MPNN encoders, the proposed SWAG model
pre-trained with LGA augmentation achieves comparable results on $6$ of the
$8$ datasets.

  * •

When compared against baselines with KGNN encoders, our proposed model
exhibits a clear dominance over $7$ of the $8$ datasets, with the COLLAB
dataset being the only exception where GCC is shown to be a better-performing
augmentation method. We found this to be also consistent with comparisons
against MPNN encoders where GCC is also quite effective.

  * •

We found two distinct self-supervision objectives under both contrastive and
non-contrastive frameworks to behave similarly, therefore partially verifying
that the performance gain of the proposed framework is mostly attributed to
the LGA instead of self-supervision objectives.

###  IV-F Visualizations of learned hidden graphs

In this section, we assess the learned SWAG model from the viewpoint of
_interpretability_. In particular, we choose two datasets, PROTEINS and RDT-B
which are representative of bio/chemo-informatics networks as well as social
networks. We randomly extract $4$ learned hidden graph from the learned SWAG
model, and plot them in figure 2. From the illustrations, it is visually clear
that hidden graphs learned from distinct datasets exhibit different types of
structural characteristics: Those learned from PROTEINS (the first row) show a
relatively frequent pattern of _rings_ with typical instantiations being
$3$-rings and $4$-rings. Meanwhile, those learned from RDT-B suggest an
alternative common pattern that is somewhat _chain-like_ , i.e., figure 2(e)
and figure 2(g) are both (isomorphic) $5$-chains. Therefore, the proposed SWAG
model offers reasonable transparency and interpretability that characterizes
learned information, which is not possible for MPNNs.

###  IV-G Ablation study

(a) (b)

Figure 3: Investigation of model performance on PROTEINS (left) and RDT-B
(right) when varying the USVT thresholding parameter $\tau$, mean performance
are plotted along with shades indicating standard deviation under $10$-fold
CV.

In this section, we assess the effect of two critical hyperparameters in the
proposed framework, namely the thresholding parameter $\tau$ involved in the
USVT procedure, as well as the number of hidden graphs $M$ in the SWAG
architecture. All the experiments are conducted on two representative
datasets, PROTEINS and RDT-B.
Effects of $\tau$ We draw insights from the renowned theoretical threshold
$2.02$ which provably recovers dense random graphs as shown in [3], and pick
the largest $\tau$ value to be slightly bigger than twice the theoretical
value as $4.2$. To accommodate for sparser graphs, we choose the smallest
$\tau$ value to be slightly greater than $0$ for thorough evaluation, and
report the results in figure 3. The results demonstrate only a mild level of
fluctuation when varying $\tau$, suggesting that the structural invariance
might be captured even with low-rank approximations of the input graph’s
adjacency matrix.

(a) (b)

Figure 4: Investigation of model performance on PROTEINS (left) and RDT-B
(right) when varying the number of hidden graphs $M$, mean performance are
plotted along with shades indicating standard deviation under $10$-fold CV.

Effects of $M$ As shown in section IV-F, the learned hidden graphs by the SWAG
architecture may exhibit recurring patterns (ignoring edge weights).
Therefore, it is of interest to investigate whether a small number of hidden
graphs suffice and whether too many hidden graphs hurts performance as they
might end up being isomorphic structures. We tested the number of hidden
graphs $M$ with a minimum of $2$ and a maximum of $24$ on the two datasets.
The results are illustrated in figure 4. The results suggest that while a
relatively large number of hidden graphs (i.e., $M>20$) may produce many
similar learned graphs in terms of graph topology, the overall predictive
performance is not significantly affected by varying $M$.

##  V Conclusion and future works

Several improvements are proposed regarding the design and learning of kernel
graph neural networks (KGNNs). We extend previous formulations of random-walk
graph neural networks into a more flexible framework and derive a novel neural
architecture SWAG that allows smoother optimization. We also initiate
explorations on applying self-supervision to KGNNs, developing the structural-
preserving graph data augmentation LGA which is beneficial to self-supervised
KGNN learning. We believe this is just the initial step toward the development
of better-performing KGNN models: It will be of great value if we may further
bridge the theory and practice of MPNNs and KGNNs together, and create better
supervised and self-supervised learning frameworks, we leave such promising
research into future works.

## References

  * [1] Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., et al. Relational inductive biases, deep learning, and graph networks.  arXiv preprint arXiv:1806.01261 (2018).
  * [2] Borgwardt, K. M., and Kriegel, H.-P. Shortest-path kernels on graphs.  In Fifth IEEE international conference on data mining (ICDM’05) (2005), IEEE, pp. 8–pp.
  * [3] Chatterjee, S. Matrix estimation by universal singular value thresholding.  The Annals of Statistics (2015), 177–214.
  * [4] Chen, T., Kornblith, S., Norouzi, M., and Hinton, G. A simple framework for contrastive learning of visual representations.  In International conference on machine learning (2020), PMLR, pp. 1597–1607.
  * [5] Chen, X., and He, K. Exploring simple siamese representation learning.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (2021), pp. 15750–15758.
  * [6] Ding, K., Xu, Z., Tong, H., and Liu, H. Data augmentation for deep graph learning: A survey.  ACM SIGKDD Explorations Newsletter 24, 2 (2022), 61–77.
  * [7] Errica, F., Podda, M., Bacciu, D., and Micheli, A. A fair comparison of graph neural networks for graph classification.  arXiv preprint arXiv:1912.09893 (2019).
  * [8] Feng, A., You, C., Wang, S., and Tassiulas, L. Kergnns: Interpretable graph neural networks with graph kernels.  In Proceedings of the AAAI Conference on Artificial Intelligence (2022), vol. 36, pp. 6614–6622.
  * [9] Gao, C., Lu, Y., and Zhou, H. H. Rate-optimal graphon estimation.  The Annals of Statistics (2015), 2624–2652.
  * [10] Gao, X., Xiao, B., Tao, D., and Li, X. A survey of graph edit distance.  Pattern Analysis and applications 13 (2010), 113–129.
  * [11] Gasteiger, J., Weißenberger, S., and Günnemann, S. Diffusion improves graph learning.  Advances in neural information processing systems 32 (2019).
  * [12] Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. Neural message passing for quantum chemistry.  In Proceedings of the 34th International Conference on Machine Learning (International Convention Centre, Sydney, Australia, 06–11 Aug 2017), D. Precup and Y. W. Teh, Eds., vol. 70 of Proceedings of Machine Learning Research, PMLR, pp. 1263–1272.
  * [13] Gligorijević, V., Renfrew, P. D., Kosciolek, T., Leman, J. K., Berenberg, D., Vatanen, T., Chandler, C., Taylor, B. C., Fisk, I. M., Vlamakis, H., et al. Structure-based protein function prediction using graph convolutional networks.  Nature communications 12, 1 (2021), 3168.
  * [14] Goldenberg, A., Zheng, A. X., Fienberg, S. E., Airoldi, E. M., et al. A survey of statistical network models.  Foundations and Trends® in Machine Learning 2, 2 (2010), 129–233.
  * [15] Hamilton, W., Ying, Z., and Leskovec, J. Inductive representation learning on large graphs.  Advances in neural information processing systems 30 (2017).
  * [16] Hou, Z., Liu, X., Cen, Y., Dong, Y., Yang, H., Wang, C., and Tang, J. Graphmae: Self-supervised masked graph autoencoders.  In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (2022), pp. 594–604.
  * [17] Kingma, D. P., and Ba, J. Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 (2014).
  * [18] Kipf, T. N., and Welling, M. Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907 (2016).
  * [19] Kipf, T. N., and Welling, M. Variational graph auto-encoders.  CoRR abs/1611.07308 (2016).
  * [20] Lei, T., Jin, W., Barzilay, R., and Jaakkola, T. Deriving neural architectures from sequence and graph kernels.  In International Conference on Machine Learning (2017), PMLR, pp. 2024–2033.
  * [21] Li, S., Wang, X., Zhang, A., Wu, Y., He, X., and Chua, T.-S. Let invariant rationale discovery inspire graph contrastive learning.  In International Conference on Machine Learning (2022), PMLR, pp. 13052–13065.
  * [22] Lin, L., Chen, J., and Wang, H. Spectral augmentation for self-supervised learning on graphs.  arXiv preprint arXiv:2210.00643 (2022).
  * [23] Morris, C., Lipman, Y., Maron, H., Rieck, B., Kriege, N. M., Grohe, M., Fey, M., and Borgwardt, K. Weisfeiler and leman go machine learning: The story so far.  arXiv preprint arXiv:2112.09992 (2021).
  * [24] Nikolentzos, G., and Vazirgiannis, M. Random walk graph neural networks.  Advances in Neural Information Processing Systems 33 (2020), 16211–16222.
  * [25] Nikolentzos, G., and Vazirgiannis, M. Geometric random walk graph neural networks via implicit layers.  In International Conference on Artificial Intelligence and Statistics (2023), PMLR, pp. 2035–2053.
  * [26] Oord, A. v. d., Li, Y., and Vinyals, O. Representation learning with contrastive predictive coding.  arXiv preprint arXiv:1807.03748 (2018).
  * [27] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems 32 (2019).
  * [28] Qiu, J., Chen, Q., Dong, Y., Zhang, J., Yang, H., Ding, M., Wang, K., and Tang, J. Gcc: Graph contrastive coding for graph neural network pre-training.  In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining (2020), pp. 1150–1160.
  * [29] Saunshi, N., Ash, J., Goel, S., Misra, D., Zhang, C., Arora, S., Kakade, S., and Krishnamurthy, A. Understanding contrastive learning requires incorporating inductive biases.  In International Conference on Machine Learning (2022), PMLR, pp. 19250–19286.
  * [30] Schraudolph, N. N., Kondor, R., and Borgwardt, K. M. Graph kernels.  Journal of Machine Learning Research 11 (2010), 1201–1242.
  * [31] Shervashidze, N., Schweitzer, P., Van Leeuwen, E. J., Mehlhorn, K., and Borgwardt, K. M. Weisfeiler-lehman graph kernels.  Journal of Machine Learning Research 12, 9 (2011).
  * [32] Shervashidze, N., Vishwanathan, S., Petri, T., Mehlhorn, K., and Borgwardt, K. Efficient graphlet kernels for large graph comparison.  In Artificial intelligence and statistics (2009), PMLR, pp. 488–495.
  * [33] Simonovsky, M., and Komodakis, N. Dynamic edge-conditioned filters in convolutional neural networks on graphs.  In Proceedings of the IEEE conference on computer vision and pattern recognition (2017), pp. 3693–3702.
  * [34] Stokes, J. M., Yang, K., Swanson, K., Jin, W., Cubillos-Ruiz, A., Donghia, N. M., MacNair, C. R., French, S., Carfrae, L. A., Bloom-Ackermann, Z., et al. A deep learning approach to antibiotic discovery.  Cell 180, 4 (2020), 688–702.
  * [35] Sun, F.-Y., Hoffmann, J., Verma, V., and Tang, J. Infograph: Unsupervised and semi-supervised graph-level representation learning via mutual information maximization.  arXiv preprint arXiv:1908.01000 (2019).
  * [36] Suresh, S., Li, P., Hao, C., and Neville, J. Adversarial graph augmentation to improve graph contrastive learning.  Advances in Neural Information Processing Systems 34 (2021), 15920–15933.
  * [37] Thakoor, S., Tallec, C., Azar, M. G., Munos, R., Veličković, P., and Valko, M. Bootstrapped representation learning on graphs.  In ICLR 2021 Workshop on Geometrical and Topological Representation Learning (2021).
  * [38] Velickovic, P., Fedus, W., Hamilton, W. L., Liò, P., Bengio, Y., and Hjelm, R. D. Deep graph infomax.  ICLR (Poster) 2, 3 (2019), 4.
  * [39] Wu, L., Lin, H., Tan, C., Gao, Z., and Li, S. Z. Self-supervised learning on graphs: Contrastive, generative, or predictive.  IEEE Transactions on Knowledge and Data Engineering (2021).
  * [40] Xu, K., Hu, W., Leskovec, J., and Jegelka, S. How powerful are graph neural networks?  arXiv preprint arXiv:1810.00826 (2018).
  * [41] Yanardag, P., and Vishwanathan, S. Deep graph kernels.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (New York, NY, USA, 2015), KDD ’15, Association for Computing Machinery, p. 1365–1374.
  * [42] Ying, Z., You, J., Morris, C., Ren, X., Hamilton, W., and Leskovec, J. Hierarchical graph representation learning with differentiable pooling.  Advances in neural information processing systems 31 (2018).
  * [43] You, Y., Chen, T., Sui, Y., Chen, T., Wang, Z., and Shen, Y. Graph contrastive learning with augmentations.  Advances in neural information processing systems 33 (2020), 5812–5823.
  * [44] Zhang, M., Cui, Z., Neumann, M., and Chen, Y. An end-to-end deep learning architecture for graph classification.  In Proceedings of the AAAI conference on artificial intelligence (2018), vol. 32.


# On the implementation of Adaptive and Filtered MHE

Federico Oliva
Dipartimento di Ing. Civile e Ing. Informatica
University of Rome “Tor Vergata”<EMAIL_ADDRESS>
&Daniele Carnevale
Dipartimento di Ing. Civile e Ing. Informatica
University of Rome “Tor Vergata”<EMAIL_ADDRESS>

###### Abstract

Optimisation-based algorithms known as Moving Horizon Estimator (MHE) have
been developed through the years. This paper illustrates the implementation of
the policy introduced in the companion paper submitted to the 18th IFAC
Workshop on Control Applications of Optimization (Oliva and Carnevale, 2022),
in which we propose two techniques to reduce the computational cost of MHEs.
These solutions mainly rely on output filtering and adaptive sampling. The use
of filters reduces the total amount of data used by MHE, shortening the length
of the moving window (buffer) and consequently decreasing the time consumption
for plant dynamics integration. Meanwhile, the proposed adaptive sampling
policy discards those sampled data that do not allow a sensible improvement of
the estimation error. Algorithms and numerical simulations are provided to
show the effectiveness of the proposed strategies.

_Keywords_ Adaptive observer, Nonlinear systems, Moving Horizon Estimator

##  1 Introduction

Observer design represents a well-known task in control theory, playing a
pivotal role in the controller design. A classic example of the
interconnections of these two problems is described by the Separation
Principle stated by Kalman in its most famous and elegant definition. As far
as observers are concerned, many different approaches have been studied and
developed over the years; from the most used linear algorithms like the
Luenberger and Kalman observers (Davis, 2002; Thrun et al., 2006), up to more
involved solutions for nonlinear systems (Thrun et al., 2006; Reif et al.,
1999; Karagiannis et al., 2008; Karagiannis and Astolfi, 2005). Another topic
of great interest in the field of observers is the adaptive approach, whose
task consists in estimating model parameters too (Luders and Narendra, 1974;
Krener and Isidori, 1983; Marino, 1990; Tyukin et al., 2013; Marine et al.,
2001; Marino and Tomei, 1992). A fascinating approach addressing a general
observation problem (both adaptive and non-adaptive) is the one proposed by
the Moving Horizon Estimators (MHE) (Michalska and Mayne, 1995; P.E. Moraal,
1995; Kang, 2006; Wynn et al., 2014; Sui et al., 2010; Suwantong et al., 2014;
Schiller et al., 2021). These are very powerful optimisation-based observers,
but they suffer from computational cost issues. These issues make their online
implementation quite hard to reach. To increase the computational speed, we
propose in (Oliva and Carnevale, 2022) to exploit output-measurement filters
to reduce the number of model integrations needed to solve the optimisation
problem at the basis of standard MHE algorithms. Moreover, we propose an
adaptive-sampling policy, making the observer capable of autonomously choosing
the "best" measurements to be considered in the moving window buffer. The
strategies proposed hinge upon MHE convergence analysis presented in (Aeyels,
1981; Kang, 2006; Glad, 1983; Menini et al., 2019, 2022).

This short work describes in detail the implementation of solutions proposed
in (Oliva and Carnevale, 2022), presenting the Algorithms to implement the
Filtered MHE and the Adaptive MHE and providing as well some simulation
results. More specifically, in 2 the Standard MHE frameworks is described,
while 3 treats both the Filtered MHE and Adaptive MHE cases. Conclusions are
drawn in 4, as well as future developments.

##  2 Moving Horizon Estimators

Consider the general framework described in (Oliva and Carnevale, 2022),
namely a nonlinear system in the following form

$\displaystyle\dot{\bm{\xi}}$ $\displaystyle=f(\bm{\xi},\bm{u})$ (1a)
$\displaystyle\bm{y}$ $\displaystyle=h(\bm{\xi},\bm{u}).$ (1b)

In this set of equations $\bm{\xi}\in\mathbb{R}^{n}$ is the state vector,
$\bm{u}\in\mathbb{R}^{m}$ is the control input, and $\bm{y}\in\mathbb{R}^{p}$
is the measured output. From this general framework, in order to define the
MHE, system (1) shall be considered in its N-lifted form. Accordingly to (P.E.
Moraal, 1995) and (Tousain et al., 2001), and following the notation
introduced in (Oliva and Carnevale, 2022), the N-lifted system generated from
(1) is described by an output operator in the following form:

$\displaystyle\text{H}_{k}(\bm{\xi}_{k-N+1},u(\cdot))\\!=\\!\\!\\!\begin{bmatrix}\\!h(\bm{\xi}_{k-N+1},\bm{u}_{k-N+1})\\\
\\!h(\phi(t_{k-N+2},t_{k-N+1},\bm{\xi}_{k-N+1},\bm{u}(\cdot)),\bm{u}_{k-N+2})\\\
\\!\vdots\\\
\\!h(\phi(t_{k},t_{k-1},\bm{\xi}_{k-1},\bm{u}(\cdot)),\bm{u}_{k})\end{bmatrix},$
(2)

where $(\text{Y}_{k},\text{U}_{k})$ are the output and input sample buffers of
length $N$ and down-sampling $N_{T_{s}}$. Moreover,
$\phi(t_{k};t_{k-1},\bm{\xi}_{k-1},\bm{u}(\tau))=\bm{\xi}(t_{k})=\bm{\xi}_{k}$
describes the solution of (1) at time $t_{k}$ with $\tau\in[t_{k-1},t_{k}]$,
and $\bm{u}(t)$ is assumed to be a piece-wise constant function of the time,
within sampling times of length $T_{s}$. Therefore, The MHE associated to (1)
through its N-lifted system can be defined as

$\displaystyle\bm{{\zeta}}_{i+1}$
$\displaystyle=\bm{\Psi}(\text{Y}_{k},\text{U}_{k},\text{H}_{k}(\bm{{\zeta}}_{i},\bm{u}(\cdot)),\bm{\zeta}_{i}),$
(3a) $\displaystyle\bm{\zeta}_{0}$
$\displaystyle=\hat{\bm{\xi}}_{k-N+1},\text{ with, }i=0\dots K,\
K\in\mathbb{N},$ (3b) $\displaystyle\hat{\bm{\xi}}_{k}$
$\displaystyle=\phi(t_{k},t_{k-N+1},\bm{\zeta}_{K},\bm{u}(\cdot)),$ (3c)

In this set of equations we consider $\hat{\bm{\xi}}_{k}$ as the estimated
state vector of $\bm{\xi}_{k}$. Instead, $i$ is the iteration number of the
optimisation algorithm defining the updated value of $\bm{\zeta}$. This value
is computed through $\Psi$, that describes a general algorithm designed to
solve optimisation problems (e.g. simplex or gradient-based solutions). This
optimisation algorithm is iteratively launched within the intervals
$[t_{k-1},t_{k}]$ for $K$ times. To sum up, the estimation problem at the
basis of MHEs consists in finding the solution to the following minimisation
problem:

$\underset{\hat{\bm{\xi}}_{k-N+1}}{\text{min}}\text{V}_{k}(\text{Y}_{k},\text{U}_{k},\hat{\text{H}}_{k},\hat{\bm{\xi}}_{k-N+1}).$
(4)

The cost $V_{k}$ has been defined as a quadratic function, namely

$V_{k}(\text{Y}_{k},\text{U}_{k},{\text{H}}_{k},\zeta)\triangleq\sum\limits_{j=1}^{N}(\text{Y}_{k}^{j}-\hat{\text{H}}_{k}^{j})^{T}W_{j}(\text{Y}_{k}^{j}-\hat{\text{H}}_{k}^{j}),$
(5)

where we considered the notation
$\hat{H}_{k}=\text{H}_{k}(\hat{\bm{\xi}}_{k-N+1},$ $\bm{u}(\cdot))$. Moreover,
$W_{i}\in\mathbb{R}^{p\times p}$ are symmetric and positive definite weight
matrices and $(\text{Y}_{k}^{j},\hat{\text{H}}_{k}^{j})$ are the j-th rows of
the matrices $Y_{k}$ and $\hat{H}_{k}$, respectively. For a detailed
discussion on the Standard MHE structure of the optimisation problem solution
and on its convergence properties, please refer to (Oliva and Carnevale,
2022). In addition to these considerations it is important to remark that if
the plant dynamics (3a) are linear, then it is possible to consider the Newton
algorithm yielding zero estimation error with just $K=2$ since $V_{k}$ is
quadratic, and therefore resulting in the dead-beat observer described in
(P.E. Moraal, 1995)111The discrete LTI plant of (1) need to be evaluated
selecting $N_{T_{s}}=1$..

Lastly, with reference to the Standard MHE proposed in (Oliva and Carnevale,
2022), a more clear description of the algorithm is proposed in pseudo-code in
1. Remember that $K$ represents the number of optimisation steps performed by
the optimisation algorithm $\bm{\Psi}$ chosen to solve (4). Moreover, note
that the integral performed in the last algorithm step on the system model is
computed numerically, as the analytical solution $\phi$ is not in general
available in closed-form.

Algorithm 1 Standard MHE

for $t_{k}=k\cdot N_{T_{s}}\cdot T_{s}\in[t_{0},t_{f}]$ do

$\triangleright$ Update buffers

$\text{Y}_{k}[1:N-1]\leftarrow\text{Y}_{k}[2:N]$,
$\text{U}_{k}[1:N-1]\leftarrow\text{U}_{k}[2:N]$

$\text{Y}_{k}[N]\leftarrow y_{k}$, $\text{U}_{k}[N]\leftarrow u_{k}$

$\triangleright$ Initialization of (3a)

$\zeta_{0}\leftarrow\hat{\bm{\xi}}_{k-N+1}$

$\triangleright$ Optimisation procedure

for $i=0,\dots,K-1$ do

$\zeta_{i+1}\leftarrow\bm{\Psi}(\text{Y}_{k},\text{U}_{k},\hat{\text{H}}_{k},\zeta_{i})$

end for

$\hat{\bm{\xi}}_{k-N+1}\leftarrow\zeta_{K}$

$\triangleright$ Propagate the estimate evaluating (3c)

for $i\leq N-1$ do

$\hat{\bm{\xi}}_{k-N+1+i}\leftarrow\int\limits_{t_{k-N+i}}^{t_{k-N+1+i}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

end for

end for

$\triangleright$ Estimate at $t=t_{k+qT_{s}}\neq t_{k}$, with
$q\in\mathbb{Z}$, $|q|<N_{T_{s}}$

$\hat{\bm{\xi}}(t)\leftarrow\int\limits_{t_{k}}^{t_{k}+qT_{s}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

##  3 Strategies for computational efficiency

This section considers the Filtered MHE and Adaptive MHE described in (Oliva
and Carnevale, 2022), and aims at better describing the code implementation.
Moreover, some considerations regarding the estimation convergence are
reported. Generally speaking the focus is on the role of parameters $N$
(buffer length) and $N_{T_{s}}$ (down-sampling). Some considerations will be
explained by simulation results on two different models, the Van der Pol
oscillator and a nonlinear system describing the dynamics of plasma waves in
presence of runaway electrons (Buratti et al., 2021; Karagiannis et al.,
2008).

### Filtered MHE

As reported in (Oliva and Carnevale, 2022) a first solution to decrease the
computational cost necessary to solve (4) consists in adding a filtered
version of the system output measurements $y_{k}$. Roughly speaking, by
increasing the conditions available at each time instants, the buffer length
$N$ can be reduced. Results in this sense have been presented on a Van der Pol
oscillator.

For the model setup and main results refer to (Oliva and Carnevale, 2022). The
major effect of adding a filtering action on the output $y$ is that the filter
state vector shall be propagated, updated and stored too during the
optimisation process. Indeed, a general filter can be described as a discrete-
time system, namely

$\displaystyle\bm{\xi_{f}}(k+1)$ $\displaystyle=\gamma(\xi_{f}(k),y(k)),$ (6a)
$\displaystyle y_{f}(k+1)$ $\displaystyle=\beta(\bm{\xi_{f}}(k+1)),$ (6b)

where $\bm{\xi_{f}}\in\mathbb{R}^{l}$ is the filter state, $y\in\mathbb{R}$ is
the signal measurement to be filtered, and $y_{f}\in\mathbb{R}$ is the actual
filtered signal. The pseudo-code relative to the Filtered MHE is presented in
2. The notation is the same used for Standard MHE.

Algorithm 2 Filtered MHE

for $t_{k}=k\cdot N_{T_{s}}\cdot T_{s}\in[t_{0},t_{f}]$ do

$\triangleright$ Output filter

$\bm{\xi_{f}}(k)\leftarrow\gamma(\bm{\xi_{f}}(k-1),y_{k})$

$y_{k}^{f}\leftarrow\beta(\bm{\xi_{f}}(k))$

$\overline{y}_{k}\leftarrow[y_{k}\ y_{k}^{f}]^{T}$

$\triangleright$ Update buffers

$\text{Y}_{k}[1:N-1]\leftarrow\text{Y}_{k}[2:N]$,
$\text{U}_{k}[1:N-1]\leftarrow\text{U}_{k}[2:N]$

$\text{Y}_{k}[N]\leftarrow\overline{y}_{k}$, $\text{U}_{k}[N]\leftarrow u_{k}$

$\triangleright$ Initialization of (3a)

$\zeta_{0}\leftarrow\hat{\bm{\xi}}_{k-N+1}$

$\triangleright$ Optimisation procedure

for $i=0,\dots,K-1$ do

$\zeta_{i+1}\leftarrow\mathbf{\Psi}(\text{Y}_{k},\text{U}_{k},\hat{\text{H}}_{k},\zeta_{i})$

end for

$\hat{\bm{\xi}}_{k-N+1}\leftarrow\zeta_{K}$

$\triangleright$ Propagate the estimate evaluating (3c)

for $i\leq N-1$ do

$\hat{\bm{\xi}}_{k-N+1+i}\leftarrow\int\limits_{t_{k-N+i}}^{t_{k-N+1+i}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

$\hat{\bm{\xi}}_{\bm{f}}(k-N+1+i)\leftarrow\sum\limits_{i=k-N+i}^{k-N+1+i}\gamma(\hat{\bm{\xi}}_{\bm{f}}(i),\hat{y}_{i})$

end for

end for

$\triangleright$ Estimate at $t=t_{k+qT_{s}}\neq t_{k}$, with
$q\in\mathbb{Z}$, $|q|<N_{T_{s}}$

$\hat{\bm{\xi}}(t)\leftarrow\int\limits_{t_{k}}^{t_{k}+qT_{s}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

$\hat{\bm{\xi}}_{\bm{f}}(k+q)\leftarrow\sum\limits_{i=k}^{k+q}\gamma(\hat{\bm{\xi}}_{\bm{f}}(i),\hat{y}_{i})$

### Adaptive MHE

The second solution proposed in (Oliva and Carnevale, 2022) to speed up MHEs
implementation consists in defining a policy according to which the observer
can automatically select the output measurements that are most informative for
the solution of (4). This selection is done accordingly to the following
indices $\delta_{k}$ and $d_{\text{V}}$, that are a good representative of the
output signal richness as well as the precision of the estimate:

$\displaystyle\Sigma_{k}$
$\displaystyle=\\{\sigma_{i}\triangleq\left\lVert\text{Y}_{k}^{i+1}-\text{Y}_{k}^{i}\right\rVert\
\text{ s.t. }i\in[1,\dots,N-1]\\},$ (7a) $\displaystyle\delta_{k}$
$\displaystyle=\sum\limits_{i=1}^{N-1}\sigma_{i}+\left\lVert
y_{k}-\text{Y}_{k}^{1}\right\rVert,$ (7b) $\displaystyle d_{\text{V}}$
$\displaystyle=\left\lVert\text{Y}_{k}-\hat{\text{H}}_{k}\right\rVert.$ (7c)

The adaptive policy is implemented by thresholding both $\delta_{k}$ and
$d_{\text{V}}$ and deciding whether to run or not the estimation accordingly.
Moreover, if no estimation is performed for a long time due to a high
precision reached, the output buffer is re-initialised to avoid excessively
time-consuming model integrations (Parameter $N_{\text{max}}$ (Oliva and
Carnevale, 2022)). Again, the pseudo-code relative to the Adaptive MHE is
presented in 3. Results have been presented on the following model, describing
the dynamics of plasma waves amplitude ($\xi_{1}$) and anisotropy of the
runaway electrons velocity distribution ($\xi_{2}$) given by (Buratti et al.,
2021):

$\displaystyle\dot{\xi}_{1}$ $\displaystyle=\epsilon(-2\xi_{1}\xi_{2}-2S+Q),$
(8a) $\displaystyle\dot{\xi}_{2}$
$\displaystyle=\epsilon(-\nu\xi_{2}+\xi_{3}(\xi_{1}\xi_{2}+S)-\gamma_{1}\dfrac{\xi_{2}}{1+\dfrac{\xi_{2}}{W_{t}}}),$
(8b) $\displaystyle\dot{\xi_{3}}$ $\displaystyle=0,$ (8c) $\displaystyle y$
$\displaystyle=\xi_{2},$ (8d)

where $y\in\mathbb{R}$ is the output, $\bm{\xi}\in\mathbb{R}^{3}$ is the state
vector. Figure 1 shows the actual measurement samples, when different values
of fixed $N_{T_{s}}$ are considered, namely $N_{T_{s}}=5$ and $N_{T_{s}}=38$.
As reported in (Oliva and Carnevale, 2022), when $N_{T_{s}}=38$, index
$\delta_{k}$ is low and the state is not correctly estimated. Indeed, this
means that the sampled trajectory is not highly informative. This can be
clearly seen in Figure 1 where the blue circles have a nearly null value
compared to the general trend of the trajectory. Clearly, an informative set
of measurements shall be sampled over the trajectory peaks, rather than on the
flat and dense intervals between them.

Figure 1: This figure compares the system trajectory (dashed line) and the
sampled measurements when $N_{T_{s}}=5$ (squares) and $N_{T_{s}}=38$
(circles).

The results of the actual Adaptive MHE are presented in Figure 2. The
simulation follows the very same setup considered in (Oliva and Carnevale,
2022). Note how the sampling is now concentrated at the beginning of the
simulation and on the peaks, in order to quickly reduce the estimation error
and then correct the estimation when the measurements were informative enough.
Whenever the estimation error increases, the observer samples the output in
correspondence of those trajectory regions in which the signal richness is
higher. The estimation precision is comparable to the Standard MHE, yet
significantly reducing the computational cost.

Figure 2: This figure shows the true state trajectory for system (8) (black
dashed line), and the sampled measures provided to the observer when the
adaptive sampling policy is used (blue big circles). Algorithm 3 Adaptive MHE

for $t_{k}=k\cdot N_{T_{s}}\cdot T_{s}\in[t_{0},t_{f}]$ do

$\triangleright$ Adaptive sampling

if AC1 $\And$ AC2 then

$\triangleright$ Check max spacing

if $T(y_{k})-T(\text{Y}_{k}^{N})\geq N_{\max}$ then

$\triangleright$ Re-initialise buffer

for $i\leq N$ do

$\text{Y}_{k}^{i}=y_{k-N+i}$

end for

else

$\triangleright$ Update buffers

$\text{Y}_{k}[1:N-1]\leftarrow\text{Y}_{k}[2:N]$

$\text{U}_{k}[1:N-1]\leftarrow\text{U}_{k}[2:N]$

$\text{Y}_{k}[N]\leftarrow y_{k}$, $\text{U}_{k}[N]\leftarrow u_{k}$

end if

$\triangleright$ Initialization of (3a)

$\zeta_{0}\leftarrow\hat{\bm{\xi}}_{k-N+1}$

$\triangleright$ Optimisation procedure

for $i=0,\dots,K-1$ do

$\zeta_{i+1}\leftarrow\mathbf{\Psi}(\text{Y}_{k},\text{U}_{k},\hat{\text{H}}_{k},\zeta_{i})$

end for

$\hat{\bm{\xi}}_{k-N+1}\leftarrow\zeta_{K}$

$\triangleright$ Propagate the estimate evaluating (3c)

for $i\leq N-1$ do

$\hat{\bm{\xi}}_{k-N+1+i}\leftarrow\int\limits_{t_{k-N+i}}^{t_{k-N+1+i}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

end for

end if

end for

$\triangleright$ Estimate at $t=t_{k+qT_{s}}\neq t_{k}$, with
$q\in\mathbb{Z}$, $|q|<N_{T_{s}}$

$\hat{\bm{\xi}}(t)\leftarrow\int\limits_{t_{k}}^{t_{k}+qT_{s}}f(\hat{\bm{\xi}}(\tau),u(\tau))d\tau$

### Filtered and Adaptive MHE

This last section focuses on the joint action of Filtered MHE and Adaptive
MHE. In fact, Filtered MHE reduces the computational speed of each
optimisation solution of (4), while Adaptive MHE reduces the number of
optimisation required to reach a good estimation precision. Therefore, if both
2 and 3 were used in combination, namely exploiting both output filtering and
adaptive sampling, we would both decrease the single optimisation time and the
total number of performed optimisations. Indeed, the speed up factor would be
even greater.

In order to test the performance of both output filtering and adaptive
sampling we combined the Filtered MHE and the Adaptive MHE on a double
pendulum, considering $K=40$ and two output filter, namely a dirty-derivative
and an integrator with loss [cfr. Oliva et al. 2022]. The model considered for
this set of simulations consists of a double pendulum. The general structure
of the equations resembles the usual model used in robotics, namely

$\mathbb{M}(q)\ddot{q}+\mathbb{V}(q,\dot{q})\dot{q}+\mathbb{G}(q)=\tau,$ (9)

where $q=(\theta_{1},\theta_{2})$ are the angular positions of the links, the
$\mathbb{M}$ term describes the system inertia, the $\mathbb{V}(q,\dot{q})$
term the friction and Coriolis effects, and $\mathbb{G}(q)$ the gravitational
force. Each pendulum was considered with length $L=1m$, and mass $M=1kg$. The
system state is
$\bm{\xi}=[q_{1},q_{2},\dot{q}_{1},\dot{q}_{2}]\in\mathbb{R}^{4}$, while the
output measurement is $y=q_{1}\in\mathbb{R}$. The estimation error norm is
presented in logarithmic scale in Figure 3.

Figure 3: This figure shows the estimation error on both the filtered (solid
line) and filtered-adaptive sampling (dashed line) cases on system (9)
estimation.

The registered total computation time was $25.4$s for the joint action of
Filtered MHE and Adaptive MHE which is nearly $60\%$ lower than the simple
Filtered MHE case reported in (Oliva and Carnevale, 2022).

##  4 Conclusions and future work

This work builds on the results presented in (Oliva and Carnevale, 2022). In
the first section 2, the Standard MHE structure was recalled and presented
more accurately in pseudo-code layout. Moreover, some remarks on the
convergence properties introduced in (Oliva and Carnevale, 2022) are reported.
More specifically, the special case of an MHE turning into a dead-beat
observer has been considered. In 3 the solutions proposed in (Oliva and
Carnevale, 2022) to speed up MHE are considered. More specifically, the
importance of filter propagation is stressed, and the Filtered MHE is
presented in pseudo-code. As far as Adaptive MHE is concerned, the simulation
results on plasma dynamics model described in (Oliva and Carnevale, 2022) are
consolidated with some considerations on both the fixed and adaptive sampling
policies. Lastly, the Adaptive MHE is presented in pseudo-code as well.

Future developments of this tool will consider a detailed analysis of the
effect of measurement noise and model uncertainties on the observer
performance and selection of filters. Another interesting point to be
developed is using this tool for control laws design, relating it to the
general MPC solution. Indeed, output tracking problems could be solved by
including control parameters in the augmented state. For instance, classic
MRAC or state feedback frameworks could be reproduced. Clearly, this
application shall be considered as separated from the state estimation
problem, because no separation principle is available for nonlinear systems.

#### Acknowledgemts

Thanks to Corrado Possieri and Mario Sassano for their comments and
fundamental insights into the development of this work.

## References

  * Oliva and Carnevale [2022] Federico Oliva and Daniele Carnevale.  Moving horizon estimator with filtering and adaptive sampling.  _submitted to the 18th IFAC Workshop Control Applications for Optimisation_ , 2022.
  * Davis [2002] Jon H. Davis.  _Luenberger Observers_ , pages 245–254.  Birkhäuser Boston, 2002.
  * Thrun et al. [2006] S. Thrun, W. Burgard, and D. Fox.  _Probabilistic Robotics_.  The MIT Press, 2006.
  * Reif et al. [1999] K. Reif, S. Gunther, E. Yaz, and R. Unbehauen.  Stochastic stability of the discrete-time extended kalman filter.  _IEEE Transactions on Automatic Control_ , 1999.
  * Karagiannis et al. [2008] Dimitrios Karagiannis, Daniele Carnevale, and Alessandro Astolfi.  Invariant manifold based reduced-order observer design for nonlinear systems.  _IEEE Transactions on Automatic Control_ , 2008.
  * Karagiannis and Astolfi [2005] D. Karagiannis and A. Astolfi.  Nonlinear observer design using invariant manifolds and applications.  In _Proceedings of the 44th IEEE Conference on Decision and Control_ , 2005.
  * Luders and Narendra [1974] G. Luders and K. Narendra.  A new canonical form for an adaptive observer.  _IEEE Transactions on Automatic Control_ , 1974.
  * Krener and Isidori [1983] Arthur J. Krener and Alberto Isidori.  Linearization by output injection and nonlinear observers.  _Systems & Control Letters_, 1983.
  * Marino [1990] R. Marino.  Adaptive observers for single output nonlinear systems.  _IEEE Transactions on Automatic Control_ , 1990.
  * Tyukin et al. [2013] Ivan Y. Tyukin, Erik Steur, Henk Nijmeijer, and Cees van Leeuwen.  Adaptive observers and parameter estimation for a class of systems nonlinear in the parameters.  _Automatica_ , 2013.
  * Marine et al. [2001] R. Marine, G.L. Santosuosso, and P. Tomei.  Robust adaptive observers for nonlinear systems with bounded disturbances.  _IEEE Transactions on Automatic Control_ , 2001.
  * Marino and Tomei [1992] R. Marino and P. Tomei.  Global adaptive observers for nonlinear systems via filtered transformations.  _IEEE Transactions on Automatic Control_ , 1992.
  * Michalska and Mayne [1995] H. Michalska and D.Q. Mayne.  Moving horizon observers and observer-based control.  _IEEE Transactions on Automatic Control_ , 1995.
  * P.E. Moraal [1995] J.W. Grizzle P.E. Moraal.  Observer design for noninear systems with discrete-time measurements.  _IEEE Transactions on Automatic Control_ , 1995.
  * Kang [2006] Wei Kang.  Moving horizon numerical observers of nonlinear control systems.  _IEEE Transactions on Automatic Control_ , 2006.
  * Wynn et al. [2014] Andrew Wynn, Milan Vukov, and Moritz Diehl.  Convergence guarantees for moving horizon estimation based on the real-time iteration scheme.  _IEEE Transactions on Automatic Control_ , 2014.
  * Sui et al. [2010] Dan Sui, Tor Arne Johansen, and Le Feng.  Linear moving horizon estimation with pre-estimating observer.  _IEEE Transactions on Automatic Control_ , 2010.
  * Suwantong et al. [2014] Rata Suwantong, Sylvain Bertrand, Didier Dumur, and Dominique Beauvois.  Stability of a nonlinear moving horizon estimator with pre-estimation.  In _2014 American Control Conference_ , 2014.
  * Schiller et al. [2021] Julian D. Schiller, Sven Knüfer, and Matthias A. Müller.  Robust stability of suboptimal moving horizon estimation using an observer-based candidate solution.  _IFAC-PapersOnLine_ , 2021.
  * Aeyels [1981] D. Aeyels.  On the number of samples necessary to achieve observability.  _Systems & Control Letters_, 1981.
  * Glad [1983] S. T. Glad.  Observability and nonlinear dead beat observers.  In _The 22nd IEEE Conference on Decision and Control_ , 1983.
  * Menini et al. [2019] Laura Menini, Corrado Possieri, and Antonio Tornambè.  Observers for linear systems by the time integrals and moving average of the output.  _IEEE Transactions on Automatic Control_ , 2019.
  * Menini et al. [2022] Laura Menini, Corrado Possieri, and Antonio Tornambe.  On the use of the time-integrals of the output in observer design for nonlinear autonomous systems.  _IEEE Transactions on Automatic Control_ , 2022.
  * Tousain et al. [2001] R. Tousain, E. van der Meche, and O. Bosgra.  Design strategy for iterative learning control based on optimal control.  _Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228)_ , 2001.  doi:10.1109/CDC.2001.980905.
  * Buratti et al. [2021] P Buratti, W Bin, A Cardinali, D Carnevale, C Castaldo, O D’Arcangelo, F Napoli, G L Ravera, A Selce, L Panaccione, and et al.  Fast dynamics of radiofrequency emission in ftu plasmas with runaway electrons.  _Plasma Physics and Controlled Fusion_ , 2021.