# On distributional graph signals

Feng Ji, Xingchao Jian and Wee Peng Tay The authors are with the School of
Electrical and Electronic Engineering, Nanyang Technological University,
639798, Singapore (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>wptay@ntu.edu.sg).

###### Abstract

Graph signal processing (GSP) studies graph-structured data, where the central
concept is the vector space of graph signals. To study a vector space, we have
many useful tools up our sleeves. However, uncertainty is omnipresent in
practice, and using a vector to model a real signal can be erroneous in some
situations. In this paper, we want to use the Wasserstein space as a
replacement for the vector space of graph signals, to account for signal
stochasticity. The Wasserstein is strictly more general in which the classical
graph signal space embeds isometrically. An element in the Wasserstein space
is called a distributional graph signal. On the other hand, signal processing
for a probability space of graphs has been proposed in the literature. In this
work, we propose a unified framework that also encompasses existing theories
regarding graph uncertainty. We develop signal processing tools to study the
new notion of distributional graph signals. We also demonstrate how the theory
can be applied by using real datasets.

###### Index Terms:

Graph signal processing, Wasserstein metric, distributional graph signals,
signal adaptive graph structures

##  I Introduction

Graph signal processing (GSP) is a rapidly growing field that studies signals
defined on graphs [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Many real-world
phenomena can be naturally represented as graphs, such as social networks,
transportation systems, and sensor networks. In GSP, the central concept is
the vector space of graph signals, and a graph signal assigns a number to each
node of a given graph. Being a vector space, we can use linear
transformations, such as the graph Fourier transform and graph filters, to
analyze graph signals and study relations among them.

However, in many practical applications, uncertainty is ubiquitous, and using
a vector to model a real signal can be erroneous. The vector space of graph
signals assumes that the signal is known exactly, but this is often not the
case in real-world scenarios. For example, in a social network, the exact
values of the attributes such as user ratings of each user may not be known
[14], or in a sensor network, the sensor readings may be uncertain due to
measurement errors or sensor variability [15]. Moreover, it is studied in [16]
that in graph neural networks (GNNs), interpreting class labels of nodes as a
graph signal can easily ignore label prediction uncertainty and the resulting
step graph signal can be highly non-smooth.

To address this issue, we propose to use the Wasserstein space [17] as a
replacement for the vector space of graph signals. An element in the
Wasserstein space is a probability distribution on the classical graph signal
space. We call such a distribution a distributional graph signal. Therefore,
uncertainty is encoded in a distributional graph signal. This provides a more
flexible and realistic approach to modeling signals on graphs, which can
account for uncertainty and stochasticity. Moreover, the Wasserstein space is
strictly more general than the classical vector space of graph signals in
which it embeds isometrically. This means that distributional graph signals
can accommodate all signals that can be represented as a vector in the
classical sense, and more, which can be represented by a probability
distribution. By considering distributional graph signals, we can develop a
more comprehensive and accurate framework for studying graph-structured data
that accounts for uncertainty. In the context of GNN, the notion of
distributional graph signals is introduced in [16] and further studied in
[18]. The distributional version of total variation and signal non-uniformity
are introduced to enhance the performance of GNNs. In this paper, we want to
propose a signal processing framework for distributional graph signals.

On the other hand, [19] proposes a signal processing framework for a
probability space of graph shift operators (see also [20] for an overview), to
address the issue that there may not be a single fixed graph topology in many
applications. Therefore, in addition to introducing the use of the Wasserstein
space for modeling graph signals, we also propose a unified framework that
encompasses existing theories regarding graph uncertainty. For this, we
introduce the notion of signal adaptive graph structures that associates a
distribution of graphs with any graph signal, so that we can construct
transformations between distributional graph signals.

In summary, we replace classical graph signals with distributional graph
signals and substitute graph topology with signal adaptive graph structures.
As a result, we have a flexible framework to deal with uncertainties in both
signals and graphs. In terms of methodology, we have to part from linear
algebra and make more use of analysis and probability theory. Therefore, our
approach has the flavor of classical Fourier theory [21] rather than that of
algebraic signal processing [22].

Our main contributions are as follows:

  * •

We introduce distributional graph signals and signal adaptive graph
structures. We develop a signal processing framework by focusing on filter
construction.

  * •

We relate the framework and the notion of conditional expectation. This allows
us to justify some key concepts introduced in [19].

  * •

We explain how classical GSP notions, such as the graph Fourier transform [1,
2, 11], convolution [2, 11], and sampling [23, 24, 25, 26, 27, 28, 29, 30],
can be interpreted using the new framework. We use examples to demonstrate
that the classical notions are special cases of their counterparts introduced
in the paper.

  * •

We demonstrate the practical utility of our proposed framework by using real
datasets. We show how the proposed approach can be used to analyze and process
graph signals with uncertain or stochastic properties. We provide experimental
results that demonstrate the effectiveness of the proposed framework.

The rest of the paper is organized as follows: In Section II, we introduce the
Wasserstein space and define the concept of distributional graph signals. In
Section III, we first introduce the notion of signal adaptive graph structures
to account for uncertainty in graph topology. Then we explain how they can be
used to define transformation between distributional graph signals. The
framework is related to the theory of conditional expectation in Section IV.
In Section V, we review classical GSP theories and describe why they are
special cases of the proposed framework. We present numerical results Section
VI and conclude in Section VII. Proofs of all results are deferred to Appendix
A.

_Notations:_ We use $\circ$ to denote function composition. Let $\mathbb{R}$
denote the set of real numbers and $M_{n}(\mathbb{R})$ be the space of
$n\times n$ real matrices. $\mathbb{E}$ is the expectation operator. Letters
$\mu,\nu,\gamma$ are used for probability distributions, while $\delta$ is for
delta distributions. We use $\mathcal{A}$ for signal adaptive graph structures
(SAGS) introduced in the paper. $G$ is used exclusively for graphs and $f$ is
used exclusively for filters. Letters in fraktur font such as
$\mathfrak{c},\mathfrak{p}$ are used to denote a pair of SAGS and a filter.
Linear operators and vectors are boldfaced.

##  II Wasserstein space and distributional graph signals

Let $V$ be a set of nodes in a network of size $|V|=n$. A classical signal on
$V$ assigns a number to each node of $V$. If an ordering of nodes in $V$ is
fixed as $V=\\{v_{1},\ldots,v_{n}\\}$, then a classical signal can be
identified with $\mathbf{x}\in\mathbb{R}^{n}$ with the $i$-th component the
number assigned to $v_{i}$. In this paper, we are interested in a
probabilistic framework. A natural way to interpret a classical signal
$\mathbf{x}$ is to view it as $\delta_{\mathbf{x}}$, the delta distribution on
$\mathbf{x}$. This prompts the following generalization of classical signals
in terms of the Wasserstein space [17].

######  Definition 1.

Let $X$ be a metric space. Define the _Wasserstein space_ $\mathcal{P}(X)$ to
be the space of (Borel) probability distributions on $X$ with finite mean and
variance. If $X=\mathbb{R}^{n}$, the space of classical graph signals on $V$,
then $\mathcal{P}(X)$ is called the space of _distributional graph signals_ on
$V$.

The main insight is that a distributional signal encodes uncertainties due to
reasons such as limitations in measurement precision, forecasting errors, and
data labeling mistakes. Hence, using distributional signals can be more
realistic than classical signals. The trade-off is that simple and effective
tools such as linear algebra are no longer available. In this paper, we shall
develop signal processing tools using mainly probability theory and analysis.
In view of this, we give $\mathcal{P}(X)$ a metric [17].

######  Definition 2.

Let $X$ be a metric space (with metric $d_{X}$) and $\mathcal{P}(X)$ be the
associated Wasserstein space. Given $\mu_{1},\mu_{2}$ in $\mathcal{P}(X)$, the
_Wasserstein metric_ 111Strictly speaking, the metric considered is the
$2$-Wasserstein metric and $2$ accounts for the power in the integral. As this
is the only version used in the paper, we omit the quantifier $2$.
$W(\mu_{1},\mu_{2})$ between $\mu_{1},\mu_{2}$ is defined by

$\displaystyle
W(\mu_{1},\mu_{2})^{2}=\inf_{\gamma\in\Gamma(\mu_{1},\mu_{2})}\int
d(x,y)^{2}\mathrm{d}\gamma(x,y),$

where $\Gamma(\mu_{1},\mu_{2})$ is the set of _couplings_ of
$\mu_{1},\mu_{2}$, i.e., the collection of probability measures on $X\times X$
whose marginals are $\mu_{1}$ and $\mu_{2}$, respectively.

Intuitively, the Wasserstein metric is the minimum amount of “work” required
to transform one probability distribution into the other, where the “work” is
the sum of the product of the amount of probability mass to be moved and the
distance that it must be moved. It is well-known that $W(\cdot,\cdot)$ makes
$\mathcal{P}(X)$ a metric space [17]. For distributional graph signals,
$\mathcal{P}(\mathbb{R}^{n})$ is complete and separable with the Wasserstein
metric. It is usually challenging to compute the Wasserstein metric for
arbitrary $\mu_{1},\mu_{2}$. However, in special cases, we have closed-form
formulas as in the following examples.

######  Example 1.

  1. (a)

If $\mu_{2}=\delta_{y}$, the delta distribution on $y\in X$, then we have the
explicit formula

$\displaystyle W(\mu_{1},\delta_{y})^{2}=\int
d_{X}(x,y)^{2}\mathrm{d}\mu_{1}(x).$

As a special case, if $\mu_{1}=\delta_{x}$ is also a delta distribution, then
$W(\delta_{x},\delta_{y})=d_{X}(x,y)$. This implies that the space of
classical graph signals $\mathbb{R}^{n}$ embeds isometrically in the space of
distributional graph signals $\mathcal{P}(\mathbb{R}^{n})$.

  2. (b)

Let $\mu_{1}=\mathcal{N}(x_{1},\Sigma_{1})$ and
$\mu_{2}=\mathcal{N}(x_{2},\Sigma_{2})$ be two non-degenerate normal
distributions on $\mathbb{R}^{n}$ with mean $x_{1},x_{2}$ and covariance
matrices $\Sigma_{1},\Sigma_{2}$ respectively. Then the Wasserstein metric is
given by

$\displaystyle W(\mu_{1},\mu_{2})^{2}$
$\displaystyle={\left\lVert{x_{1}-x_{2}}\right\rVert}^{2}$
$\displaystyle+\text{trace}\big{(}\Sigma_{1}+\Sigma_{2}-2(\Sigma_{2}^{1/2}\Sigma_{1}\Sigma_{2}^{1/2})^{1/2}\big{)}.$

Therefore, fitting data in terms of the Wasserstein metric requires one to
consider fitting covariance in addition to fitting the mean.

We have introduced the fundamental object $\mathcal{P}(\mathbb{R}^{n})$ to be
studied in the paper. However, we have not yet described contributions from
graphs. We explain how graphs enter into the overall picture in the next
section.

##  III The Bayesian perspective of distributional graph operators

In this section, we want to introduce a signal processing framework for
distributional graph signals that generalizes classical GSP. There are two
aspects of the framework: (1) to encode graph structural information, and (2)
to describe (distributional) signal transformations. Each topic occupies one
of the following subsections. For concreteness, we do not present the theory
in full generality. A more general framework is briefly outlined in Appendix
B.

###  III-A Signal adaptive graph structures

Let $\mathcal{G}_{n}$ be the set of undirected graphs without multiple edges
on (ordered) $n$ vertices $V=\\{v_{1},\ldots,v_{n}\\}$. The graphs can be
weighted. Therefore, there is an embedding of $\mathcal{G}_{n}$ in
$M_{n}(\mathbb{R})$, the space matrices of size $n$. More specifically, the
embedding associates a $G\in\mathcal{G}_{n}$ with its weighted adjacency
matrix $\mathbf{A}_{G}$, where the $(i,j)$-th entry of $\mathbf{A}_{G}$ is the
weight between $v_{i}$ and $v_{j}$. As $M_{n}(\mathbb{R})$ is measurable with
Lebesgure $\sigma$-algebra, it induces a $\sigma$-algebra on
$\mathcal{G}_{n}$. Moreover, $\mathcal{G}_{n}$ is equipped with the subspace
topology.

The key insight is that we allow the graph structure to depend on the signal,
moreover, it can be random. We formally introduce the following notion.

######  Definition 3.

A _signal adaptive graph structures_ (SAGS) assigns to each
$\mathbf{x}\in\mathbb{R}^{n}$ a probability distribution $\nu_{\mathbf{x}}$ on
$\mathcal{G}_{n}$. Denote it by
$\mathcal{A}=(\nu_{\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$.

To associate the notion with the Bayesian theory, consider the product space
$\mathbb{R}^{n}\times\mathcal{G}_{n}$. It is a measurable space with the
product $\sigma$-algebra. The probability distribution $\nu_{\mathbf{x}}$ can
be interpreted as the (conditional) distribution on $\mathcal{G}_{n}$ given
$\mathbf{x}\in\mathbb{R}^{n}$. Therefore, for any distributional graph signal
$\mu\in\mathcal{P}(\mathbb{R}^{n})$, we have the an associated distribution
$\mathcal{A}^{*}(\mu)$ on $\mathbb{R}^{n}\times\mathcal{G}_{n}$ defined by

$\displaystyle\mathcal{A}^{*}(\mu)(g)=\int\int
g(\mathbf{x},G)\mathrm{d}\nu_{\mathbf{x}}(G)\mathrm{d}\mu(\mathbf{x}),$

for any compactly supported continuous function $g$ on
$\mathbb{R}^{n}\times\mathcal{G}_{n}$. The distribution $\mathcal{A}^{*}(\mu)$
is uniquely determined by the integral formula by the Riesz–Markov–Kakutani
representation theorem [21]. The expression reminds us of the law of total
probability if $\nu_{\mathbf{x}}$ is interpreted as the conditional
distribution. We now give some examples.

######  Example 2.

  1. (a)

If $\nu_{\mathbf{x}}=\nu$, i.e., independent of $\mathbf{x}$, then we have the
setup of [19]. We call it a _constant SAGS_. Moreover, if $\nu=\delta_{G}$ for
a single $G\in\mathcal{G}_{n}$, we recover the classical GSP. A further
generalization is given next.

  2. (b)

A SAGS $\mathcal{A}=(\nu_{\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$ is
_locally constant_ for almost every $\mathbf{x}\in\mathbb{R}^{n}$, there is an
open neighborhood $U_{\mathbf{x}}$ of $\mathbf{x}$ such that for every
$\mathbf{y}\in U_{\mathbf{x}}$, we have $\nu_{\mathbf{y}}=\nu_{\mathbf{x}}$.
Intuitively, for such an $\mathcal{A}$, the signal space $\mathbb{R}^{n}$ can
be (almost) partitioned into open subsets on each of which $\mathcal{A}$ is a
constant.

Analogous to these examples, a SAGS $\mathcal{A}$ encodes the graph structural
information. It tells us for a given signal $\mathbf{x}$, the most suitable
graph structures on $V$, according to $\nu_{\mathbf{x}}$, to process
$\mathbf{x}$. In the next subsection, we describe how distributional signal
transformation is performed in this framework.

###  III-B Distributional signal transformations

Recall that in classical GSP, given a graph $G$, one constructs linear
transformations or filters by using the structure of $G$. For example, one may
first fix a graph shift operator GSO $\mathbf{S}$ such as the adjacency matrix
or the Laplacian of $G$. Then one applies an algebraic construction such as
taking polynomials in $\mathbf{S}$ to construct desired filters $\mathbf{F}$.
The entire process $G\mapsto\mathbf{F}$ can be summarized as a map from
$\mathcal{G}_{n}$ to $M_{n}(\mathbb{R})$ if we omit the intermediate steps.

Based on this prototype, we call any measurable function $f:\mathcal{G}_{n}\to
M_{n}(\mathbb{R})$ a _pre-filter_ or a _pre-transformation_. It induces a
measurable function
$\widetilde{f}:\mathbb{R}^{n}\times\mathcal{G}_{n}\to\mathbb{R}^{n}$ by
$\widetilde{f}(\mathbf{x},G)=f(G)(\mathbf{x})$, using the fact that $f(G)\in
M_{n}(\mathbb{R})$ and $f(G)(\mathbf{x})$ is the ordinary matrix operation.

Given any probability distribution $\mu$ on
$\mathbb{R}^{n}\times\mathcal{G}_{n}$, the map $\widetilde{f}$ induces the
_pushforward_ distribution $f_{*}(\mu)$ on $\mathbb{R}^{n}$. More
specifically, for any measurable subset $U$ of $\mathbb{R}^{n}$, we have

$\displaystyle f_{*}(\mu)(U)=\mu(\widetilde{f}^{-1}(U)).$ (1)

We do not yet call $f$ a filter or a transformation because we want to impose
more constraints on $f$ regarding the distributional graph signals
$\mathcal{P}(\mathbb{R}^{n})$.

######  Definition 4.

Given an SAGS $\mathcal{A}$, a measurable $f:\mathcal{G}_{n}\to
M_{n}(\mathbb{R})$ is a _filter_ or a _transformation_ with respect to
(w.r.t.) $\mathcal{A}$ if for any distributional graph signal
$\mu\in\mathcal{P}(\mathbb{R}^{n})$, the distribution
$f_{*}\circ\mathcal{A}^{*}(\mu)$ is also a distributional graph signal, i.e.,
$f_{*}\circ\mathcal{A}^{*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$. For
convenience, we use $\mathfrak{c}$ to denote the pair $(\mathcal{A},f)$ and
write $\mathfrak{c}_{*}(\mu)$ for $f_{*}\circ\mathcal{A}^{*}(\mu)$, if no
confusions arise.

The map
$\mathfrak{c}_{*}(\mu):\mathcal{P}(\mathbb{R}^{n})\to\mathcal{P}(\mathbb{R}^{n})$
satisfies the following explicit integral formula:

$\displaystyle\mathfrak{c}_{*}(\mu)(g)=\int\int
g\big{(}f(G)(\mathbf{x})\big{)}\mathrm{d}\nu_{\mathbf{x}}(G)\mathrm{d}\mu(\mathbf{x}),$

for any compactly supported continuous function on $\mathbb{R}^{n}$.

As we have mentioned, the space distributional graph signals
$\mathcal{P}(\mathbb{R}^{n})$ is not linear and we want to focus on the
analytic perspective of filters. Recall that one of the most desired analytic
properties of a linear map (in functional analysis) is continuity, or
equivalently boundedness [31]. In our framework, we also want to study when
the map
$\mathfrak{c}_{*}(\mu):\mathcal{P}(\mathbb{R}^{n})\to\mathcal{P}(\mathbb{R}^{n})$
induced by a filter $f$ is continuous.

For this, we notice that $\mathcal{A}$ and $f$ give rise to a probability
distribution $f_{\mathcal{A}}(\mathbf{x})$ on $M_{n}(\mathbb{R})$ given
$\mathbf{x}\in\mathbb{R}^{n}$. Recall
$\mathcal{A}=(\nu_{\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$, and the
distribution $f_{\mathcal{A}}(\mathbf{x})$ is given by

$\displaystyle
f_{\mathcal{A}}(\mathbf{x})(U)=f_{*}(\nu_{\mathbf{x}})(U)=\nu_{\mathbf{x}}(f^{-1}(U)),$

for any measurable subset $U$ of $M_{n}(\mathbb{R})$. Intuitively, the SAGS
$\mathcal{A}$ associates a family of probable graphs (according to
$\nu_{\mathbf{x}}$) to each $\mathbf{x}\in\mathbb{R}^{n}$ and the filter $f$
turns them into a family of probable linear maps that in terms of
$f_{\mathcal{A}}(\mathbf{x})$. We endow $M_{n}(\mathbb{R})$ with the operator
norm.

######  Theorem 1.

Let $K$ be a compact subset of $\mathbb{R}^{n}$. If $f_{\mathcal{A}}$ (when
restricted to $K$) is a continuous function from $K$ to
$\mathcal{P}(M_{n}(\mathbb{R}))$, then restricted to $\mathcal{P}(K)$,
$\mathfrak{c}_{*}:\mathcal{P}(K)\to\mathcal{P}(\mathbb{R}^{n})$ is uniformly
continuous.

In many practical situations, it is reasonable to assume that signals belong
to a compact and hence bounded subset of $\mathbb{R}^{n}$. In such a case, the
condition of the theorem is not restrictive. We also have a version of the
continuity result without the compactness assumption.

######  Theorem 2.

If $f_{\mathcal{A}}$ is Lipschitz continuous, then $\mathfrak{c}_{*}(\mu)$ is
continuous at any $\mu\in\mathcal{P}(\mathbb{R}^{n})$ with finite $6$-th
moments.

We remark that the condition on finite $6$-th moment can be further improved.
However, it is sufficient for us as it already includes essential cases such
as compactly supported distributions and (mixed) Gaussian distributions.

We have the following consequence of the result. It is known (e.g., [17]) that
finite point distributions are dense in $\mathcal{P}(\mathbb{R}^{n})$, i.e,
for any distributional graph signals $\mu\in\mathcal{P}(\mathbb{R}^{n})$,
there is a sequence $(\mu_{i})_{i\geq 1}$ of distributional graph signals each
supported on finitely many points such that $\mu_{i}\to\mu,i\to\infty$. If
$\mu$ has bounded $6$-th moment and $\mathcal{X}$ satisfies the conditions of
Theorem 2, then by continuity, we have
$\mathfrak{c}_{*}(\mu_{i})\to\mathfrak{c}_{*}(\mu),i\to\infty$. This means
that knowledge of the filter at delta distributions tells us a lot about the
filter at more general distributions.

##  IV Conditional expectations

In this section, we propose construction based on a $\mathfrak{c}$ that is
related to conditional expectations [32]. The approximation result Theorem 3
justifies many constructions in [19].

As we have seen in the previous section, given a pair of SAGS and a filter
$\mathfrak{c}=(\mathcal{A},f)$, we have
$\mathfrak{c}_{*}=f_{*}\circ\mathcal{A}^{*}:\mathcal{P}(\mathbb{R}^{n})\to\mathcal{P}(\mathbb{R}^{n})$.
On the other hand, for any measurable function
$g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $\mu\in\mathcal{P}(\mathbb{R}^{n})$,
pushforward (cf. (1)) induces a probability distribution $g_{*}(\mu)$ on
$\mathbb{R}^{n}$. We call $g$ _bounded_ if
$g_{*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$ for any
$\mu\in\mathcal{P}(\mathbb{R}^{n})$. Denote the set of bounded measurable
functions by $B(\mathbb{R}^{n})$. For example, a linear transformation is
bounded and hence belongs to $B(\mathbb{R}^{n})$.

In classical GSP when $\mathcal{A}$ is a constant delta distribution,
$\mathfrak{c}_{*}$ is induced by the pushforward of a linear transformation.
It is easier to study such a map coming directly from a function on the more
familiar space $\mathbb{R}^{n}$. However, for a general
$\mathfrak{c}=(\mathcal{A},f)$, it is not always true that
$\mathfrak{c}_{*}=g_{*}$ for some $g\in B(\mathbb{R}^{n})$. Nevertheless, it
is possible to find good approximations of $\mathfrak{c}_{*}$. For this, we
introduce a function $e_{\mathfrak{c}}$ as follows.

To construct $e_{\mathfrak{c}}$, assume that
$\mathcal{A}=(\nu_{\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$. For
$\mathbf{x}\in\mathbb{R}^{n}$, we define

$\displaystyle e_{\mathfrak{c}}(\mathbf{x})=\int
f(G)(\mathbf{x})\mathrm{d}\nu_{\mathbf{x}}(G)=\int\mathbf{y}\mathrm{d}\mathfrak{c}_{*}(\delta_{\mathbf{x}})(\mathbf{y}).$
(2)

It is related to conditional expectation as follows. Let
$p:\mathbb{R}^{n}\times\mathcal{G}_{n}\to\mathbb{R}^{n}$ be the projection to
the first component. Recall that for any $\mu\in\mathcal{P}(\mathbb{R}^{n})$,
we have constructed the distribution $\mathcal{A}^{*}(\mu)$ on
$\mathbb{R}^{n}\times\mathcal{G}_{n}$. In this respect, both
$p:\mathbb{R}^{n}\times\mathcal{G}_{n}\to\mathbb{R}^{n}$ and
$f:\mathbb{R}^{n}\times\mathcal{G}_{n}\to\mathbb{R}^{n}$ can be viewed as
random variables on the sample space $\mathbb{R}^{n}\times\mathcal{G}_{n}$. It
is well known that there is a condition expectation
$e_{f}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ such that
$e_{f}(\mathbf{x})=e_{\mathfrak{c}}(\mathbf{x})$ up to a set with $\mu$
measure $0$. Due to this fact, the promised approximation property of
$e_{\mathfrak{c}}$ reads as follows.

######  Theorem 3.

For $\mathfrak{c}=(\mathcal{A},f)$, the function $e_{\mathfrak{c}}$ is
measurable and belongs to $B(\mathbb{R}^{n})$. Moreover, for any $g\in
B(\mathbb{R}^{n})$ and subset $S\subset\mathbb{R}^{n}$, the following holds:

$\displaystyle\sup_{\text{supp}(\mu)\subset
S}W\big{(}e_{\mathfrak{c},*}(\mu),\mathfrak{c}_{*}(\mu)\big{)}\leq\sup_{\text{supp}(\nu)\subset
S}W\big{(}g_{*}(\nu),\mathfrak{c}_{*}(\nu)\big{)},$

where $W$ is the Wasserstein metric and the supreme is taken over $\mu$ (resp.
$\nu$) in $\mathcal{P}(\mathbb{R}^{n})$ supported in $S$.

We give some examples.

######  Example 3.

If $\mathcal{A}$ is a constant SAGS with the common probability measure $\nu$
(cf. Example 2a), then $e_{\mathfrak{c}}$ is the linear transformation given
by the operator

$\displaystyle e_{\mathfrak{c}}(\cdot)=\int f(G)(\cdot)\mathrm{d}\nu(G).$

Similarly, if $\mathcal{A}$ is a locally constant SAGS (cf. Example 2b), then
outside a subset of measure $0$, the function $e_{\mathfrak{c}}$ is piecewise
linear, i.e., for each $\mathbf{x}$ there is an open neighborhood of
$\mathbf{x}$ on which $e_{\mathfrak{c}}$ is linear.

The construction of $e_{\mathfrak{c}}$ enjoys other analytic properties.

######  Lemma 1.

Consider a sequence $\mathfrak{c}_{i}=(\mathcal{A}_{i},f_{i}),i\geq 1$. If
there is a $\mathfrak{c}=(\mathcal{A},f)$ such that
$f_{i,\mathcal{A}_{i}}(\mathbf{x})\to f_{\mathcal{A}}(\mathbf{x})$ as
$i\to\infty$, then $e_{\mathfrak{c}_{i}}(\mathbf{x})\to
e_{\mathfrak{c}}(\mathbf{x})$.

Intuitively, the lemma says that the construction $e_{\mathfrak{c}}$ is
“continuous” in $\mathfrak{c}$.

For the rest of this section, we discuss some algebraic properties of
$e_{\mathfrak{c}}$. Unlike classical GSP, Wasserstein spaces are not linear.
However, we can still define binary operations such as addition, analogous to
the sum of random variables.

Let $f_{1},f_{2}$ be filters w.r.t. SAGSs
$\mathcal{A}_{1}=(\nu_{1,\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$ and
$\mathcal{A}_{2}=(\nu_{2,\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{n}}$
respectively and denote $(\mathcal{A}_{i},f_{i})$ by $\mathfrak{c}_{i},i=1,2$.
We define the addition $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ by the
property

$\displaystyle\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)(g)=$
$\displaystyle\int\int
g\big{(}f_{1}(G_{1})(\mathbf{x})+f_{2}(G_{2})(\mathbf{x})\big{)}\mathrm{d}\nu_{1,\mathbf{x}}\times\nu_{2,\mathbf{x}}(G_{1},G_{2})\mathrm{d}\mu(\mathbf{x}),$

for any continuous function $g$ with compact support on $\mathbb{R}^{n}$. The
addition $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ allows us to combine
filters w.r.t. different SAGSs. From the expression, we see its similarity to
the sum of random variables. Analogous to (2), its associated “conditional
expectation” $e_{\mathfrak{c}_{1}\boxplus\mathfrak{c}_{2}}$ is given by the
integral

$\displaystyle
e_{\mathfrak{c}_{1}\boxplus\mathfrak{c}_{2}}(\mathbf{x})=\int\mathbf{y}\mathrm{d}\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\delta_{\mathbf{x}})(\mathbf{y}).$

Scalar multiplication is simpler: given $r\in\mathbb{R}$ and
$\mathfrak{c}=(\mathcal{A},f)$, then $r\mathfrak{c}$ denote the pair
$(\mathcal{A},rf)$. The construction of $e_{\mathfrak{c}}$ from $\mathfrak{c}$
respects addition and scalar multiplication.

######  Lemma 2.

The addition $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ is well defined.
Moreover, for any $r\in\mathbb{R}$, we have
$e_{r\mathfrak{c}_{1}\boxplus\mathfrak{c}_{2}}=re_{\mathfrak{c}_{1}}+e_{\mathfrak{c}_{2}}.$

In the next section, we revisit some key concepts of graph signal processing
theories and interpret them with the new framework.

##  V GSP theories revisited

In this section, we describe how we may understand some of the most important
GSP concepts in the proposed work, in view of how they are perceived before.
We mainly base on [1] and [19], which are briefly reviewed in Example 2 and
Example 3.

Given a graph $G\in\mathcal{G}_{n}$, recall that the Fourier transform in the
classical GSP is defined as the orthogonal base change w.r.t. an eigenbasis
$\mathbf{U}_{G}$ of a prescribed graph shift operator $\mathbf{S}_{G}$ (e.g.,
the adjacency or Laplacian matrices). An interpretation is that each
eigenvector of $\mathbf{S}_{G}$ accounts for a level of signal smoothness
quantified by its eigenvalue.

Given a SAGS $\mathcal{A}$, if we want to imitate the classical construction,
we may copy the classical recipe and define the filter222As $\mathbf{U}_{G}$
defines an orthogonal transformation that is norm preserving, we have
$\mathfrak{c}_{*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$ for
$\mu\in\mathcal{P}(\mathbb{R}^{n})$. $\phi:\mathcal{G}_{n}\to
M_{n}(\mathbb{R}),G\mapsto\mathbf{S}_{G}\mapsto\mathbf{U}_{G}$. Let
$\mathfrak{c}=(\mathcal{A},\phi)$. Such a transform allows us to probe signal
smoothness by incorporating probabilistic information. If $\mathcal{A}$ is a
constant SAGS, then the Fourier transform introduced in [19] is nothing but
$e_{\mathfrak{c}}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ (cf. Theorem 3) when the
notion of distributional graph signal is not yet introduced.

From this explicit construction, we see the route to follow. Suppose a
classical construction can be described by a function $f:\mathcal{G}_{n}\to
M_{n}(\mathbb{R})$. It also defines a filter if
$c_{*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$ for
$\mu\in\mathcal{P}(\mathbb{R}^{n})$, where $\mathfrak{c}=(\mathcal{A},f)$. For
another important example, if $f:\mathcal{G}_{n}\to M_{n}(\mathbb{R})$ is a
filter such that $f(G)$ is a polynomial in (a prescribed GSO)
$\mathbf{S}_{G}$, then
$\mathfrak{c}_{*}:\mathcal{P}(\mathbb{R}^{n})\to\mathcal{P}(\mathbb{R}^{n})$
is a _convolution_. Similarly to Fourier transform, the notion of convolution
introduced in [19] is nothing but $e_{\mathfrak{c}}$ for constant
$\mathcal{A}$. A special family of convolutions leads to the theory of
sampling. Such a convolution takes the form $\rho:\mathcal{G}_{n}\to
M_{n}(\mathbb{R})$, where each $\rho(G)$ is the orthogonal projection matrix
to the direct sum of a subcollection of eigenspaces of $\mathbf{S}_{G}$. Let
$\mathfrak{p}=(\mathcal{A},\rho)$. Inspired by [33] and [19], for
$\epsilon>0$, a distributional graph signal $\mu$ is called
_$(\epsilon,\mathfrak{p})$ -invariant_ if
$W(\mu,\mathfrak{p}_{*}(\mu))<\epsilon$. Recovery requires one to estimate
such a $\mu$ based on its partial sampled observations.

######  Example 4.

If $\mathcal{A}$ is constant and
$\mu=\delta_{\mathbf{x}},\mathbf{x}\in\mathbb{R}^{n}$, then
${\left\lVert{\mathbf{x}-e_{\mathfrak{p}}(\mathbf{x})}\right\rVert}\leq
W(\mu,\mathfrak{p}_{*}(\mu))$ by the Jensen inequality [34] and Example 1a.
Therefore, if $\mu$ is $(\epsilon,\mathfrak{p})$-invariant, then $\mathbf{x}$
is $(\epsilon,e_{\mathfrak{p}})$-bandlimited in the sense of [19], where an
explicit recovery scheme is given. If $\mathcal{A}$ is locally constant, a
brief discussion is given in Appendix C.

Though it is impossible to discuss all important GSP concepts exhaustively,
some essential ones have been covered. In the next section, we use numerical
experiments to demonstrate how the framework of the paper can be applied in
practice.

##  VI Experimental results

###  VI-A MNIST: examples of distributional graph signals

In this experiment, we showcase visualizations of distributional graph signals
by summarising samples of each digit from $0$ to $9$ in the MNIST
dataset333http://yann.lecun.com/exdb/mnist/ as a distributional signal. We
preprocess the sample images by introducing i.i.d Gaussian noise to each
pixel. The graph $G$ used is $28\times 28$ 2D-lattice.

We consider $2$ different approaches.

  1. (I)

Edgewise Gaussian $\mu_{E}$ (abbreviated as “Edgewise”): We learn from samples
the joint Gaussian distribution of pairs of pixel values for each edge of the
graph $G$. To draw a sample, we give $G$ an acyclic orientation with a single
root. We draw a pixel value at the root using its marginal. For any directed
edge, if the pixel value at the tail is already known, then the value at the
head is drawn according to the conditional distribution derived from the joint
distribution of the edge. The pixel values are averaged if a node is the head
of multiple directed edges. The approach captures more refined pairwise signal
relations in closed vicinity.

  2. (II)

Joint Gaussian $\mu_{G}$ (abbreviated as “Joint”): It is the joint Gaussian
distribution of values at all the pixels that fits the samples. To draw a
sample, we just draw from the joint distribution. The approach is based on a
global perspective on the entire graph.

We draw samples from both $\mu_{E}$ and $\mu_{G}$. From the sample images
shown in (the right half of) Fig. 1 and Fig. 2, we see that non of the
approaches generate images with reasonable equality. For example, for
$\mu_{E}$, the digits are not even recognizable. However, this does not
necessarily mean that the distributions contain no useful information. We
apply a thresholding function. The resulting samples are also shown in (the
left half of) Fig. 1 and Fig. 2. We see that now the digits are clearly
recognizable. Moreover, $\mu_{E}$, the only distribution that leverages the
graph structure, generates arguably the sharpest image of digits.

Fig. 1: Samples drawn from $\mu_{E}$. The first half of the images are
obtained by thresholding the second half of the images. Fig. 2: Samples drawn
from $\mu_{G}$. The first half of the images are obtained by thresholding the
second half of the images.

We further investigate by resorting to the primary purpose that the dataset
created: digit recognition. We take a base neural network model and perform
the following two tasks.

  1. (a)

In the first task, we train the network with varying sizes of training sets.
Then we test with the original test data (of size $10000$), as well as test
data generated from distributional signals (Edgewise and Joint). The
distributional signals are obtained using the original test data. The results
are shown in Fig. 3. From the results, we see that the accuracy of samples
from distributional signals: Edgewise is the highest in all the cases. This
may suggest hidden statistical features might be captured by the
distributional signals.

  2. (b)

In the second task, we consider augmentation by distributional graph signals.
We train the network with a small training set (of varying size $\leq 2000$).
Moreover, we augment the dataset with $10000$ samples generated from
distributional signals (Edgewise and Joint). Unlike the previous task, to get
the distributional signals, we make use of a small portion of the original
training dataset. We show the test results (in Fig. 4) on the original test
dataset both with and without the augmentation. From the results, we notice
that augmentation with distributional graph signals does significantly improve
the test accuracy when the number of training samples is small. Moreover,
using augmented samples: Joint has a better overall performance.

Fig. 3: Test accuracy of the original test data and samples from the
distributional graph signals.

Fig. 4: Test accuracy using training samples of small size (blue curve) and
augmented training samples (red curve).

The investigations suggest that samples from distributional signals: Edgewise
might capture more details of the digits while using distributional signals:
Joint can be more robust. To verify the last claim, we consider neural network
adversarial attacks FGSM and PGD [35, 36]. More specifically, we use the
original test dataset, while for the training dataset, we either use $10000$
samples from the original training dataset or $10000$ samples drawn from
distributional signals: Joint. Test accuracies are shown in Table I. We see
that in general, the distributional approach can better resist adversarial
attacks.

TABLE I:

Perturbation | $0.1$ | $0.2$ | $0.3$
---|---|---|---
Original | $34.3\%$ | $16.5\%$ | $11.9\%$
Joint | $44.2\%$ | $25.2\%$ | $20.9\%$
(a) FGSM attack

Perturbation | $0.05$ | $0.1$ | $0.2$ | $0.3$
---|---|---|---|---
Original | $72.8\%$ | $9.85\%$ | $0.16\%$ | $0.16\%$
Joint | $96.0\%$ | $71.8\%$ | $4.74\%$ | $0\%$
(b) PGD attack

###  VI-B Weather dataset: filters and prediction

In this example, we consider filter learning for signal prediction. We use the
US weather station network444http://www.ncdc.noaa.gov/data-access/ with 194
nodes, and they are connected by a $20$-NN graph $G$. Signals are temperature
reading over a year. We want to learn a convolution filter $\mathbf{F}$ in the
normalized Laplacian $\widetilde{\mathbf{L}_{G}}$ of degree up to $2$ that
predicts temperature $4$ days or $7$ days in the future. In the setting of the
paper, $\mathbf{F}$ is $f(G)$ as in Definition 4. We compare two approaches.

  1. (a)

Classical GSP: we estimate $\mathbf{F}$ that best predicts readings $4$ days
(or $7$ days) in the future for $7$ consecutive days via a least mean square
optimization.

  2. (b)

Distributional signals: we summarize readings in $7$ consecutive days as a
(Guassian) distribution. The filter $\mathbf{F}$ fits the distribution with
the distribution of readings $4$ days (or $7$ days) in the future, by
minimizing Wasserstein distance (Example 1). In particular, the variance of
$7$ days reading at each station is taken into consideration.

We perform the experiments for the readings in the $1$st half and $2$nd half
of the years separately. We (uniformly) randomly sample a small fraction of
groups of signals, with each group consisting of readings from $7$ consecutive
days. For each group $\mathfrak{g}$ of readings, a filter
$\mathbf{F}_{\mathfrak{g}}$ is estimated using one of the two approaches
described above. More specifically, let $a_{\mathfrak{g}}$ be the average
reading (over all the stations) on the first day of the group $\mathfrak{g}$.
An insight of Example 2b is that filters may change with signals. In this
spirit, we may estimate $\mathbf{F}_{\mathfrak{g}}$ that depends on
$\mathfrak{g}$, e.g., the coefficients of $\mathbf{F}_{\mathfrak{g}}$ (in
$\widetilde{\mathbf{L}_{G}}$) are themselves polynomials in
$a_{\mathfrak{g}}$. In summary, given any number $a$, we can output a filter
$\mathbf{F}$ that is a degree $2$ polynomial in $\widetilde{\mathbf{L}_{G}}$,
whose coefficients are (learned and hence known) functions in $a$.

In testing, given any signal $\mathbf{x}$, we compute $a_{\mathbf{x}}$ as the
average reading (over all the stations) of $\mathbf{x}$. We hence obtain a
filter $\widetilde{\mathbf{F}_{\mathbf{x}}}$ using $a_{\mathbf{x}}$. The
filter $\widetilde{\mathbf{F}_{\mathbf{x}}}$ is used for prediction and the
performance is evaluated by the SNR of the predicted signal against the actual
reading in the future. In the experiments below, we may consider either degree
$2$ or degree $0$ polynomials in $a_{\mathfrak{g}}$ for filter coefficients.
Degree $0$ is equivalent to the filter unchanged for different signals.

In summary, we may propose approaches that consider distributional graph
signals or classical (statistic) graph signals, denoted by (d) or (s)
respectively for convenience. Moreover, the filter coefficients can either
vary as polynomials in the mean of the signals or remain constant. The two
situations are denoted by (p) and (c) for convenience. Altogether, we have
four different combinations of approaches (d)(p), (d)(c), (s)(p), (s)(c).
Their performance, with $20\%$ of training samples, is shown in Fig. 5. We see
that the distributional signal approaches have much better performance in all
the cases.

Fig. 5: The performance of the four different approaches with $20\%$ of
training samples.

To further compare (d)(p) and (d)(c), we vary the fraction of training samples
and compute and record (in Fig. 6) the average SNR for the two approaches. We
see that for the $4$ days prediction when the predictions are supposed to be
more accurate (as compared with $7$ days), (d)(p) is better than (d)(c) by a
small margin but with a clear overall trend, i.e., it is preferable to let the
filters change according to signals in the spirit of Example 2b. On the other
hand, for the $7$ days prediction, (d)(p) and (d)(c) have comparable
performance. In summary, using (d)(p) is at least as effective as the other
approaches and can even be beneficial in some cases.

Fig. 6: The average SNR comparison between (d)(p) and (d)(c) against different
fractions of training samples.

###  VI-C Brain ECoG dataset: anomaly detection

In this experiment, we apply the framework of the paper to anomaly detection.
We use the brain ECoG
dataset.555https://math.bu.edu/people/kolaczyk/datasets.html For each of the
eight subjects in the dataset, there are $76$ sensors recording (normalized)
brain ECoG signals in a time-series of $4000$ time-stamps. There are two
signal types: pre-ict and ict signals. We consider ict signals abnormal.

We segment the entire time-series into sub-intervals of size $10$ each. The
$10$ time stamps can be modeled by the path graph $P$ on $10$ nodes. Suppose
there is a connection $H$ among the sensors. Then there is the graph
$G=H\times P$ of size $760$, with each node $v$ of $G=(V,E)$ corresponding to
a pair $(s,t)$ where $s$ is a sensor and $t$ is a time-stamp. Different $H$
results in different $G$. A graph signal $\mathbf{x}$ consists of sensor
readings for $10$ consecutive time-stamps.

We consider $\mathfrak{c}_{j}=(\mathcal{A}_{j},f_{j}),j=1,2$ with SAGS
$\mathcal{A}_{j}=(\mu_{j,\mathbf{x}})_{\mathbf{x}\in\mathbb{R}^{760}}$ defined
as follows. We assume that for different subjects, their signals are disjoint,
i.e., no two patients can have the same ECoG signal. Therefore,
$\mathbb{R}^{n}$ is decomposed as $\mathbb{R}^{n}=\cup_{1\leq i\leq
8}C_{i}\cup C^{\prime}$, where $C_{i}$ are all possible signals of the $i$-th
subject and $C^{\prime}$ is the complement of $\cup_{1\leq i\leq 8}C_{i}$ that
plays no role in the problem. Therefore, we effectively consider locally
constant SAGSs.

For any $G=H\times P$, let $\mathbf{L}_{G}$ be its Laplacian. The filters
$f_{1},f_{2}$ are the high pass filter for the range from $730$ to $760$
w.r.t. $\mathbf{L}_{G}$. Hence, $\mathfrak{c}_{1}$ and $\mathfrak{c}_{2}$
differ only in $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$. To define
$\mathcal{A}_{j}$, first for $\mathbf{x}\in C^{\prime}$, let
$\mu_{j,\mathbf{x}}$ be supported on $G_{0}=H_{0}\times P$ for any fixed
$H_{0}$ for convenience, as it is not used in the sequel. For $\mathbf{x}\in
C_{i}$, the empirical distribution of $H$ is estimated as in [19] Section VII
D using $10\%$ of data as training samples. It is lengthy to give the details
here, we just point out that in the estimation, one needs to specify a (graph)
frequency range of $\mathbf{L}_{G}$. We choose the frequency range from $0$ to
$50$ for $\mathcal{A}_{1}$ and from $50$ to $100$ for $\mathcal{A}_{2}$.

We apply $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ (Section IV) to each
of the training set to obtain an empirical distribution in $\mathbb{R}^{2}$
and fit it with a mixed Gaussian with at most $3$ components. Though a mixed
Gaussian may not be the best choice of distribution, we only need to know the
positions of the peaks. For each subject, we randomly sample $10$ signals
either all pre-ict or ict to form a signal test instance. The map
$\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ is applied to the instance, and
the resulting empirical distribution in $\mathbb{R}^{2}$ is again fitted with
a mixed Gaussian distribution. It is compared with the base distributions.
Examples for subject $1$ are shown in Fig. 7 and Fig. 8) (more are shown in
the supplementary materials). We see that for both ict and pre-ict signals of
each subject, the peak positions of Gaussian obtained from the test instances
match well with those of the base distributions.

Fig. 7: Pre-ict signals, subjects $1$: the base distribution and $3$ test
instances. Fig. 8: Ict signals, subjects $1$: the base distribution and $3$
test instances.

The above observation suggests the following anomaly detection scheme. For the
setup, we randomly choose a subject and a condition. Moreover, from the
corresponding dataset (for the chosen patient and condition), we randomly draw
a small number of samples ($\leq 8$). Let $\mu$ be the discrete distribution
supported on the chosen samples. Using the method described earlier, we
estimate the peak locations of the mixed Gaussian that fits
$\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)$ and compare to peaks
locations of the $16$ base distributions. The comparison uses the Euclidean
norm between the peak locations. The condition, either ict (abnormal) or pre-
ict (normal), is declared using the corresponding condition of the base
distribution models that are closest in average peak distance. We run the
experiments for sample size $1,\ldots,8$ and compute the detection accuracy
based on $200$ runs for each subject and condition. The results are shown in
Fig. 9. We see that the accuracy increases rapidly if we increase the sample
size. With $\geq 6$ samples, the accuracy is already $\approx 95\%$ or higher.

Fig. 9: Accuracy of anomaly detection.

##  VII Conclusions

In this work, we proposed a new approach to modeling and processing graph
signals using distributional graph signals to account for signal
stochasticity. Our proposed framework unifies existing approaches, provides a
more flexible and realistic approach to modeling uncertain graph signals and
graph topologies jointly, and has potential applications in various domains.
The results of our experiments demonstrate the effectiveness of the proposed
approach. We hope that this work can contribute to advancing the field of GSP
by inspiring further research on the use of probability spaces in signal
processing.

##  Appendix A Proofs of theoretical results

###### Proof of Theorem 1.

We first remark that as $K$ is compact, so is $\mathcal{P}(K)$ by the
Prokhorov theorem and the Skorokhod representation theorem [37]. Therefore, on
$K$ and $\mathcal{P}(K)$, any continuous function is also uniformly
continuous. We first show that if $f_{\mathcal{A}}$ is (uniformly) continuous,
then the restriction of $\mathfrak{c}_{*}$ to $K$ is (uniformly) continuous.

Consider $\mathbf{x}_{1},\mathbf{x}_{2}\in K$. Let
$\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}$ be a distribution on
$M_{n}(\mathbb{R})\times M_{n}(\mathbb{R})$ that realizes
$W(f_{\mathcal{A}}(\mathbf{x}_{1}),f_{\mathcal{A}}(\mathbf{x}_{2}))$, i.e.,

$\displaystyle
W(f_{\mathcal{A}}(\mathbf{x}_{1}),f_{\mathcal{A}}(\mathbf{x}_{2}))^{2}=\int{\left\lVert{\mathbf{M}_{1}-\mathbf{M}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2}).$

The condition that $f_{\mathcal{A}}$ is uniformly continuous means that: if
$\mathbf{x}_{1},\mathbf{x}_{2}$ are close enough in Euclidean distance, then
the above integral can be arbitrarily small. Moreover, by uniform continuity,
there is a uniform upper bound on the $W(f_{\mathcal{A}}(\mathbf{x}),0)$ for
$\mathbf{x}\in K$.

To estimate
$W(\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}}),\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}}))$,
consider

$\displaystyle M_{n}(\mathbb{R})\times
M_{n}(\mathbb{R})\to\mathbb{R}^{n}\times\mathbb{R}^{n},(\mathbf{M}_{1},\mathbf{M}_{2})\mapsto(\mathbf{M}_{1}\mathbf{x}_{1},\mathbf{M}_{2}\mathbf{x}_{2})$

and let $\gamma^{\prime}_{\mathbf{x}_{1},\mathbf{x}_{2}}$ be the pushforward
probability distribution of $\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}$ on
$\mathbb{R}^{n}\times\mathbb{R}^{n}$. Moreover, based on the construction, the
marginals of $\gamma^{\prime}_{\mathbf{x}_{1},\mathbf{x}_{2}}$ are
$\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}})$ and
$\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}})$. We have:

$\displaystyle
W(\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}}),\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}}))^{2}$
$\displaystyle\leq$
$\displaystyle\int{\left\lVert{\mathbf{z}_{1}-\mathbf{z}_{2}}\right\rVert}^{2}\mathrm{d}\gamma^{\prime}_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{z}_{1},\mathbf{z}_{2})$
$\displaystyle=$
$\displaystyle\int{\left\lVert{\mathbf{M}_{1}\mathbf{x}_{1}-\mathbf{M}_{2}\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2})$
$\displaystyle\leq$
$\displaystyle\int{\left\lVert{\mathbf{M}_{1}\mathbf{x}_{1}-\mathbf{M}_{2}\mathbf{x}_{1}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2})$
$\displaystyle+\int{\left\lVert{\mathbf{M}_{2}\mathbf{x}_{1}-\mathbf{M}_{2}\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2})$
$\displaystyle\leq$
$\displaystyle(\int{\left\lVert{\mathbf{M}_{1}-\mathbf{M}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2}))({\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2})$
$\displaystyle+(\int{\left\lVert{\mathbf{M}_{2}}\right\rVert}^{2}df_{\mathcal{A}}(\mathbf{x}_{2}))({\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2})$
$\displaystyle=$
$\displaystyle(\int{\left\lVert{\mathbf{M}_{1}-\mathbf{M}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{M}_{1},\mathbf{M}_{2}))({\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2})$
$\displaystyle+W(f_{\mathcal{A}}(\mathbf{x}_{2}),0)^{2}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}.$

The last sum can be arbitrarily small if $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$
are close enough because we have noticed that
$W(f_{\mathcal{A}}(\mathbf{x}),0)$ is uniformly bounded. Moreover, it is
independent of the location of $\mathbf{x}_{1}$ in $K$. Therefore,
$\mathfrak{c}_{*}$ is uniformly continuous when restricted to $K$. This
further implies that for $\epsilon>0$, there is $B_{\epsilon}$ depending only
on $\epsilon$ such that if $\mathbf{x}_{1},\mathbf{x}_{2}\in K$ satisfy
${\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}\geq\epsilon$, then
$W(\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}}),\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}}))^{2}\leq
B_{\epsilon}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}$.
Moreover, there is a $C_{\epsilon}$ also depending only on $\epsilon$ such
that if ${\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}<\epsilon$,
then
$W(\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}}),\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}}))^{2}\leq
C_{\epsilon}$. As $\epsilon\to 0$, $C_{\epsilon}\to 0$. We also remark that
based on the expression, if the compact set $K$ is contained in the ball of
radius $R$ (centered at the origin) in $\mathbb{R}^{n}$, then
$B_{\epsilon}=O(R)$. This will be used in the next proof.

Consider general $\mu,\mu^{\prime}$ on $K\subset\mathbb{R}^{n}$. Let $\eta$ be
a distribution on $\mathbb{R}^{n}\times\mathbb{R}^{n}$ that realizes
$W(\mu,\mu^{\prime})$ and $\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}$ be defined
earlier for
$(\mathbf{x}_{1},\mathbf{x}_{2})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$. Define
a distribution $\eta^{\prime}$ on $\mathbb{R}^{n}\times\mathbb{R}^{n}$ by the
following integral equation. For any compactly supported continuous function
$g$ on $\mathbb{R}^{n}\times\mathbb{R}^{n}$, $\eta^{\prime}$ satisfies:

$\displaystyle\int
g(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta^{\prime}(\mathbf{w}_{1},\mathbf{w}_{2})$
$\displaystyle=$ $\displaystyle\int\int
g(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2}).$

We verify that the marginals of $\eta^{\prime}$ are $\mathfrak{c}_{*}(\mu)$
and $\mathfrak{c}_{*}(\mu^{\prime})$ respectively. Let
$p:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ be the projection to
either component. Consider any compactly supported continuous function $g$ on
$\mathbb{R}^{n}$. We have

$\displaystyle\int
g(\mathbf{w}_{1})\mathrm{d}p_{*}(\eta^{\prime})(\mathbf{w}_{1})=\int
g(\mathbf{w}_{1})\mathrm{d}\eta^{\prime}(\mathbf{w}_{1},\mathbf{w}_{2})$
$\displaystyle=$ $\displaystyle\int\int
g(\mathbf{w}_{1})d\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$ $\displaystyle\int\int
g(\mathbf{w}_{1})d\mathfrak{c}_{*}({\delta_{\mathbf{x}_{1}}})(\mathbf{w}_{1})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$ $\displaystyle\int\int
g(\mathbf{w}_{1})d\mathfrak{c}_{*}({\delta_{\mathbf{x}_{1}}})(\mathbf{w}_{1})\mathrm{d}\mu(\mathbf{x}_{1})$
$\displaystyle=$ $\displaystyle\int
g(\mathbf{w}_{1})\mathrm{d}\mathfrak{c}_{*}(\mu)(\mathbf{w}_{1}).$

This proves the claim.

For any $\epsilon>0$, we estimate

$\displaystyle
W(\mathfrak{c}_{*}(\mu),\mathfrak{c}_{*}(\mu^{\prime}))^{2}\leq\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\eta^{\prime}(\mathbf{w}_{1},\mathbf{w}_{2})$
$\displaystyle=$
$\displaystyle\int\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$
$\displaystyle\int_{{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}\geq\epsilon}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+$
$\displaystyle\int_{{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}<\epsilon}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\leq$ $\displaystyle\int
B_{\epsilon}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})+\int
C_{\epsilon}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$ $\displaystyle=$
$\displaystyle B_{\epsilon}W(\mu,\mu^{\prime})+C_{\epsilon}.$

Therefore, as long as $\epsilon$ (chosen first) is small enough and
$W(\mu,\mu^{\prime})$ is small enough,
$W(\mathfrak{c}_{*}(\mu),\mathfrak{c}_{*}(\mu^{\prime}))^{2}$ can be
arbitrarily small. This proves the theorem. ∎

###### Proof of Theorem 2.

The structure of the proof follows that of the proof of Theorem 1. Following
the argument of Theorem 1 and using Lipschitz continuity of $f_{\mathcal{A}}$,
for any compact subset $K$ of $\mathbb{R}^{n}$, there is a $B_{K}$ depending
only on $K$ such that
$W(\mathfrak{c}_{*}(\delta_{\mathbf{x}_{1}}),\mathfrak{c}_{*}(\delta_{\mathbf{x}_{2}}))^{2}\leq
B_{K}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}$ for any
$\mathbf{x}_{1},\mathbf{x}_{2}\in K$. Moreover, if $K$ is contained in the
ball centered at the origin with radius $R$, then $B_{K}=O(R)$.

Let $\mu\in\mathcal{P}(\mathbb{R}^{n})$ have finite $6$-th moment and
$\epsilon>0$. By the Markov inequality, there is a closed ball $K_{\epsilon}$
(centered at the origin) with radius $R_{\epsilon}=o(1/\epsilon)$ such that
$\int_{\mathbf{x}\notin
K_{\epsilon}}{\left\lVert{\mathbf{x}}\right\rVert}^{2}d\mu(\mathbf{x})\leq\epsilon$.
Let $K^{\prime}_{\epsilon}$ be the ball (centered at the origin) with radius
$2R_{\epsilon}$. Consider any $\mu^{\prime}\in\mathcal{P}(\mathbb{R}^{n})$
such that $W(\mu,\mu^{\prime})\leq\epsilon$. Let $\eta$ be the distribution on
$\mathbb{R}^{n}\times\mathbb{R}^{n}$ that realizes $W(\mu,\mu^{\prime})$. We
estimate:

$\displaystyle\epsilon\geq$
$\displaystyle\int{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\geq$ $\displaystyle\int_{\mathbf{x}_{1}\in
K_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\int_{\mathbf{x}_{1}\notin K_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\geq$ $\displaystyle\frac{1}{4}\int_{\mathbf{x}_{1}\in
K_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\frac{1}{2}\int_{\mathbf{x}_{1}\notin
K_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle-\int_{\mathbf{x}_{1}\notin K_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\geq$ $\displaystyle\frac{1}{4}\int_{\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})-\int_{\mathbf{x}_{1}\notin
K_{\epsilon}}{\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$ $\displaystyle\frac{1}{4}\int_{\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}d\mu^{\prime}(\mathbf{x}_{2})-\int_{\mathbf{x}_{1}\notin
K_{\epsilon}}{\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2}d\mu(\mathbf{x}_{1}).$

Therefore, $\int_{\mathbf{x}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}}\right\rVert}^{2}d\mu^{\prime}(\mathbf{x})\leq
8\epsilon$.

Since it is assumed that $f_{\mathcal{A}}$ is Lipschitz, there is $B_{0}$ such
that
$W(f_{\mathcal{A}}(\mathbf{x}_{1}),f_{\mathcal{A}}(\mathbf{x}_{1}))^{2}\leq
B_{0}{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}$ for every
$\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{n}$. For
$\mu^{\prime}\in\mathcal{P}(\mathbb{R})^{n}$ such that
$W(\mu,\mu^{\prime})\leq\epsilon$, choose $K^{\prime}_{\epsilon}$ and hence
$B_{K^{\prime}_{\epsilon}}=o(1/\epsilon)$ as earlier in the proof. Moreover,
let $\eta$ and
$\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}},\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{n}$
be as in the proof of Theorem 1, the same estimation yields:

$\displaystyle W(\mathfrak{c}_{*}(\mu),\mathfrak{c}_{*}(\mu^{\prime}))^{2}$
$\displaystyle\leq$
$\displaystyle\int\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$ $\displaystyle\int_{\mathbf{x}_{1}\in
K^{\prime}_{\epsilon},\mathbf{x}_{2}\in
K^{\prime}_{\epsilon}}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+$ $\displaystyle\int_{\mathbf{x}_{1}\in
K^{\prime}_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+$ $\displaystyle\int_{\mathbf{x}_{1}\notin
K^{\prime}_{\epsilon},\mathbf{x}_{2}\in
K^{\prime}_{\epsilon}}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$ $\displaystyle\int_{\mathbf{x}_{1}\in
K^{\prime}_{\epsilon},\mathbf{x}_{2}\in
K^{\prime}_{\epsilon}}\int{\left\lVert{\mathbf{w}_{1}-\mathbf{w}_{2}}\right\rVert}^{2}\mathrm{d}\gamma_{\mathbf{x}_{1},\mathbf{x}_{2}}(\mathbf{w}_{1},\mathbf{w}_{2})\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\int_{\mathbf{x}_{1}\in
K^{\prime}_{\epsilon},\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}W(f_{\mathcal{A}}(\mathbf{x}_{1}),f_{\mathcal{A}}(\mathbf{x}_{1}))^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\int_{\mathbf{x}_{1}\notin
K^{\prime}_{\epsilon},\mathbf{x}_{2}\in
K^{\prime}_{\epsilon}}W(f_{\mathcal{A}}(\mathbf{x}_{1}),f_{\mathcal{A}}(\mathbf{x}_{1}))^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\leq$ $\displaystyle
B_{K^{\prime}_{\epsilon}}\int{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\int_{\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}4B_{0}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle+\int_{\mathbf{x}_{1}\notin
K^{\prime}_{\epsilon}}4B_{0}{\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle\leq$ $\displaystyle
B_{K^{\prime}_{\epsilon}}W(\mu,\mu^{\prime})$
$\displaystyle+4B_{0}\Big{(}\int_{\mathbf{x}_{2}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\mu^{\prime}(\mathbf{x}_{2})+\int_{\mathbf{x}_{1}\notin
K^{\prime}_{\epsilon}}{\left\lVert{\mathbf{x}_{1}}\right\rVert}^{2}\mathrm{d}\mu(\mathbf{x}_{1})\Big{)}$
$\displaystyle\leq$
$\displaystyle(B_{K^{\prime}_{\epsilon}}+32B_{0})\epsilon.$

As $(B_{K^{\prime}_{\epsilon}}+32B_{0})\epsilon=o(1)$, the distance
$W(\mathfrak{c}_{*}(\mu),\mathfrak{c}_{*}(\mu^{\prime}))$ can be arbitrarily
small as long as $\epsilon\to 0$ and $W(\mu,\mu^{\prime})<\epsilon$. The
theorem is proved. ∎

###### Proof of Theorem 3.

For $\mathfrak{c}=(\mathcal{A},f)$, we first show that $e_{\mathfrak{c}}$ is
measurable. Let $C$ be any compact subset of $\mathbb{R}^{n}$ and $\mu_{C}$ be
the uniform distribution on $C$. It induces the measure
$\mathcal{A}^{*}(\mu_{C})$ on $Y=\mathbb{R}^{n}\times M_{n}(\mathbb{R})$.
Moreover, it is easy to verify that
$p_{*}\circ\mathcal{A}^{*}(\mu_{C})=\mu_{C}$, where $p:Y\to\mathbb{R}^{n}$ is
the projection.

We view $Y$ as the sample space with probability distribution
$\mathcal{A}^{*}(\mu_{C})$ and measurable functions $p,f$ as random variables.
Let $e_{C}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be the associated conditional
expectation. By the construction, we have $e_{C}=e_{\mathfrak{c}}$ on $C$ and
$e_{C}=0$ on the complement $\mathbb{R}^{n}\backslash C$. Moreover, $e_{C}$ is
measurable w.r.t. the measure $\mu_{C}$. However, as $\mu_{C}$ is uniform,
$e_{C}$ is also measurable w.r.t. the Lebesgue measure.

Let $C_{1}\subset C_{2}\subset\ldots\subset C_{i}\subset\ldots$ be a sequence
of compact subsets of $\mathbb{R}^{n}$ such that $\cup_{i\geq
1}C_{i}=\mathbb{R}^{n}$. Then $(e_{C_{i}})_{i\geq 1}$ is a sequence of
measurable functions whose pointwise limit is $e_{\mathfrak{c}}$. Therefore,
$e_{\mathfrak{c}}$ is also measurable.

As a consequence, given any distribution $\mu$ on
$\mathcal{P}(\mathbb{R}^{n})$, pushforward of $e_{\mathfrak{c}}$ induces a
distribution $e_{\mathfrak{c},*}(\mu)$. We need to show that
$e_{\mathfrak{c},*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$ in order to claim that
$e_{\mathfrak{c},*}$ is well defined as a map
$\mathcal{P}(\mathbb{R}^{n})\to\mathcal{P}(\mathbb{R}^{n})$. Consider the
Jensen inequality [34]:

$\displaystyle{\left\lVert{\mathbb{E}_{G\sim\mu_{\mathbf{x}}}f(G)(\mathbf{x})}\right\rVert}^{2}\leq\mathbb{E}_{G\sim\mu_{\mathbf{x}}}{\left\lVert{f(G)(\mathbf{x})}\right\rVert}^{2},\mathbf{x}\in\mathbb{R}^{n}.$
(3)

For any $\mu\in\mathcal{P}(\mathbb{R}^{n})$, to show that
$e_{\mathfrak{c},*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$, it suffices to check
that

$\displaystyle\int{\left\lVert{\mathbb{E}_{G\sim\mu_{\mathbf{x}}}f(G)(\mathbf{x})}\right\rVert}^{2}d\mu(\mathbf{x})<\infty,$

as finiteness of mean follows from that of $\mathfrak{c}_{*}(\mu)$ and
linearity of expectation. However, by (3), the left-hand side is bounded by

$\displaystyle\int\mathbb{E}_{G\sim\mu_{\mathbf{x}}}{\left\lVert{f(G)(\mathbf{x})}\right\rVert}^{2}d\mu(\mathbf{x})<\infty,$

due to the assumption that
$\mathfrak{c}_{*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$ and its $2$nd moment is
finite.

To show the claimed inequality, let $S$ be a subset $\mathbb{R}^{n}$ and
$\mu\in\mathcal{P}(\mathbb{R}^{n})$ be supported on $S$. Consider $\eta$ the
pushforward measure of $\mathcal{A}^{*}(\mu)$ on
$\mathbb{R}^{n}\times\mathbb{R}^{n}$ via the map:
$Y\to\mathbb{R}^{n}\times\mathbb{R}^{n},(\mathbf{x},G)\mapsto(e_{\mathfrak{c}}\circ
p(\mathbf{x},G),f(G)\mathbf{x})=(e_{\mathfrak{c}}(\mathbf{x}),f(G)\mathbf{x})$.
The marginals of $\eta$ are

$\displaystyle(e_{\mathfrak{c}}\circ
p)_{*}(\mathcal{A}^{*}(\mu))=e_{\mathfrak{c},*}(\mu)\text{ and
}f_{*}\circ\mathcal{A}^{*}(\mu)=\mathfrak{c}(\mu)$

respectively. Therefore,

$\displaystyle W(e_{\mathfrak{c},*}(\mu),\mathfrak{c}(\mu))^{2}$
$\displaystyle\leq$
$\displaystyle\int{\left\lVert{\mathbf{x}_{1}-\mathbf{x}_{2}}\right\rVert}^{2}\mathrm{d}\eta(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=$
$\displaystyle\int{\left\lVert{e_{\mathfrak{c}}(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2}\mathrm{d}\mathcal{A}^{*}(\mu)(\mathbf{x},G)$
$\displaystyle=$
$\displaystyle\mathbb{E}_{(\mathbf{x},G)\sim\mathcal{A}^{*}(\mu)}{\left\lVert{e_{\mathfrak{c}}(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{(\mathbf{x},G)\sim\mathcal{A}^{*}(\mu)}{\left\lVert{g(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2},$

The last inequality holds as $e_{\mathfrak{c}}$ is the conditional expectation
(up to a set of measure $0$) w.r.t. $\mathcal{A}^{*}(\mu)$ on $Y$ [32].

To estimate the right-hand-side, we have

$\displaystyle\mathbb{E}_{(\mathbf{x},G)\sim\mathcal{A}^{*}(\mu)}{\left\lVert{g(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2}$
$\displaystyle\leq$ $\displaystyle\sup_{\mathbf{x}\in
S}\mathbb{E}_{(\mathbf{x},G)\sim\mathcal{A}^{*}(\delta_{\mathbf{x}})}{\left\lVert{g(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2}.$

As $g(\mathbf{x})$ is independent of $G$, we have

$\displaystyle\mathbb{E}_{(\mathbf{x},G)\sim\mathcal{A}^{*}(\delta_{\mathbf{x}})}{\left\lVert{g(\mathbf{x})-f(G)(\mathbf{x})}\right\rVert}^{2}$
$\displaystyle=$ $\displaystyle
W\big{(}\delta_{g(\mathbf{x})},\mathfrak{c}(\delta_{\mathbf{x}})\big{)}^{2}=W\big{(}g_{*}(\delta_{\mathbf{x}}),\mathfrak{c}(\delta_{\mathbf{x}})\big{)}^{2}$
$\displaystyle\leq$ $\displaystyle\sup_{\text{supp}(\nu)\subset
S}W\big{(}g_{*}(\nu),\mathfrak{c}(\nu)\big{)}^{2}.$

The result follows. ∎

###### Proof of Lemma 1.

Suppose $\gamma$ is a distribution on $M_{n}(\mathbb{R})\times
M_{n}(\mathbb{R})$ that realizes
$W(f_{i,\mathcal{A}_{i}}(\mathbf{x}),f_{\mathcal{A}}(\mathbf{x}))$. We
estimate

$\displaystyle
W(f_{i,\mathcal{A}_{i}}(\mathbf{x}),f_{\mathcal{A}}(\mathbf{x}))^{2}{\left\lVert{\mathbf{x}}\right\rVert}^{2}$
$\displaystyle=$
$\displaystyle\int{\left\lVert{\mathbf{M}_{1}-\mathbf{M}_{2}}\right\rVert}^{2}{\left\lVert{\mathbf{x}}\right\rVert}^{2}\mathrm{d}\gamma(\mathbf{M}_{1},\mathbf{M}_{2})$
$\displaystyle\geq$
$\displaystyle\big{(}\int{\left\lVert{(\mathbf{M}_{1}-\mathbf{M}_{2})\mathbf{x}}\right\rVert}\mathrm{d}\gamma(\mathbf{M}_{1},\mathbf{M}_{2})\big{)}^{2}$
$\displaystyle\geq$
$\displaystyle{\left\lVert{\int(\mathbf{M}_{1}-\mathbf{M}_{2})\mathbf{x}\mathrm{d}\gamma(\mathbf{M}_{1},\mathbf{M}_{2})}\right\rVert}^{2}$
$\displaystyle=$
$\displaystyle{\left\lVert{\int\mathbf{M}_{1}\mathbf{x}\mathrm{d}\gamma(\mathbf{M}_{1},\mathbf{M}_{2})-\int\mathbf{M}_{2}\mathbf{x}\mathrm{d}\gamma(\mathbf{M}_{1},\mathbf{M}_{2})}\right\rVert}^{2}$
$\displaystyle=$
$\displaystyle{\left\lVert{\int\mathbf{M}_{1}\mathbf{x}\mathrm{d}f_{i,\mathcal{A}_{i}}(\mathbf{x})(\mathbf{M}_{1})-\int\mathbf{M}_{2}\mathbf{x}\mathrm{d}f_{\mathcal{A}}(\mathbf{x})(\mathbf{M}_{2})}\right\rVert}^{2}$
$\displaystyle=$
$\displaystyle{\left\lVert{e_{\mathfrak{c}_{i}}(\mathbf{x})-e_{\mathfrak{c}}(\mathbf{x})}\right\rVert}^{2}.$

Therefore, if $f_{i,\mathcal{A}_{i}}(\mathbf{x})\to
f_{\mathcal{A}}(\mathbf{x}),i\to\infty$, then
$e_{\mathfrak{c}_{i}}(\mathbf{x})\to e_{\mathfrak{c}}(\mathbf{x})$. ∎

###### Proof of Lemma 2.

To show $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}$ is well-defined, we
want to prove that
$\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)\in\mathcal{P}(\mathbb{R}^{n})$
if $\mu\in\mathcal{P}(\mathbb{R}^{n})$. Let $g$ be the function
$\mathbf{x}\mapsto{\left\lVert{\mathbf{x}}\right\rVert}^{2}$, and we have

$\displaystyle\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)(g)=$
$\displaystyle\int\int{\left\lVert{f_{1}(G_{1})(\mathbf{x})+f_{2}(G_{2})(\mathbf{x})}\right\rVert}^{2}\mathrm{d}\nu_{1,\mathbf{x}}\times\nu_{2,\mathbf{x}}(G_{1},G_{2})\mathrm{d}\mu(\mathbf{x})$
$\displaystyle\leq\int\int
2{\left\lVert{f_{1}(G_{1})(\mathbf{x})}\right\rVert}^{2}\mathrm{d}\nu_{1,\mathbf{x}}\times\nu_{2,\mathbf{x}}(G_{1},G_{2})\mathrm{d}\mu(\mathbf{x})$
$\displaystyle+\int\int
2{\left\lVert{f_{2}(G_{2})(\mathbf{x})}\right\rVert}^{2}\mathrm{d}\nu_{1,\mathbf{x}}\times\nu_{2,\mathbf{x}}(G_{1},G_{2})\mathrm{d}\mu(\mathbf{x})$
$\displaystyle=2\mathfrak{c}_{1,*}(\mu)(g)+2\mathfrak{c}_{2,*}(\mu)(g).$

Therefore, $\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)$ has finite
variance as both $f_{1}$ and $f_{2}$ are finites. Similarly, if we choose $g$
to be the identity function, we see that
$\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)$ has finite mean. Hence,
$\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\mu)$ is in
$\mathcal{P}(\mathbb{R}^{n})$.

To verify the algebraic identity, we compute

$\displaystyle
e_{r\mathfrak{c}_{1}+\mathfrak{c}_{2}}(\mathbf{x})=\int\mathbf{y}\mathrm{d}r\mathfrak{c}_{1,*}\boxplus\mathfrak{c}_{2,*}(\delta_{\mathbf{x}})(\mathbf{y})$
$\displaystyle=$ $\displaystyle\int
rf_{1}(G_{1})(\mathbf{x})+f_{2}(G_{2})(\mathbf{x})d\nu_{1,\mathbf{x}}\times\nu_{2,\mathbf{x}}(G_{1},G_{2})$
$\displaystyle=$ $\displaystyle\int
rf_{1}(G_{1})(\mathbf{x})d\nu_{1,\mathbf{x}}(G_{1})$ $\displaystyle+\int
f_{2}(G_{2})(\mathbf{x})d\nu_{2,\mathbf{x}}(G_{2})$ $\displaystyle=$
$\displaystyle
e_{r\mathfrak{c}_{1}}(\mathbf{x})+e_{\mathfrak{c}_{2}}(\mathbf{x}).$

∎

##  Appendix B A category theoretical perspective

Category theory [38] is a branch of mathematics that deals with the abstract
study of structures and relationships between objects. It provides a framework
for organizing mathematical concepts and objects. Category theory aims to
identify common patterns and structures across different mathematical
disciplines and provide a unified language for talking about these structures.

A category $\mathcal{C}$ consists of the following entities:

  * •

A class $ob(\mathcal{C})$, the _objects_ of $\mathcal{C}$.

  * •

For every pair objects $C_{1},C_{2}$, a class $mor(C_{1},C_{2})$ of
_morphisms_ from $C_{1}$ to $C_{2}$.

  * •

For any triple of objects $C_{1},C_{2},C_{3}$, there is the composition
$\circ:mor(C_{1},C_{2})\times mor(C_{2},C_{3})\to mor(C_{1},C_{3})$ expressed
as $\circ(f,g)=g\circ f$ such that

    1. (a)

$f\circ(g\circ h)=(f\circ g)\circ h$.

    2. (b)

For each object $C$, there is the identity morphism $1_{C}\in mor(C,C)$ such
that $f\circ 1_{C}=f=1_{C}\circ f$.

The most relevant category to traditional GSP is $\text{Vect}_{\mathbb{R}}$,
the category of finite dimensional vector spaces. In
$\text{Vect}_{\mathbb{R}}$, the objects are finite dimensional
$\mathbb{R}$-vector spaces, and the morphisms between a pair of vector spaces
are the linear transformations between them. More generally, we have Meas the
category of measurable spaces. The morphisms between two measurable spaces are
measurable functions.

The framework of the paper can also be described using a category
$\mathcal{C}$. We highlight some essential ideas. The objects are measurable
spaces. A morphism (up to certain equivalence) $\mathfrak{c}=(Y,f_{1},f_{2})$
between two objects $X_{1},X_{2}$ consists of a measurable space $Y$ and
measurable functions $f_{1}:Y\to X_{1}$ and $f_{2}:Y\to X_{2}$ such that the
following holds (cf. [19] Section V):

  * •

For each $x\in X_{1}$, there a probability distribution $\mu_{x}$ on
$f_{1}^{-1}(x)$. The collection of _fiberwise_ distributions $(\mu_{x})_{x\in
X_{1}}$ induces for any probability measure $\mu$ on $X_{1}$, a probability
measure on $f_{1}^{*}(\mu)$ on $Y$.

  * •

Let ${f_{2}}_{*}$ be the pushforward map of probability measures on $Y$. The
composition ${f_{2}}_{*}\circ f_{1}^{*}$ is well-defined as a map
$\mathcal{P}(X_{1})\to\mathcal{P}(X_{2})$.

In the setup, the primary example of $X_{1}$ is the graph signal space
$\mathbb{R}^{n}$. Graph structural information is encoded in
$f_{1}:Y\to\mathbb{R}^{n}$ and $(\mu_{x})_{x\in\mathbb{R}^{n}}$, when $Y$
consists of pairs $(x,G)$ with $G$ a graph of size $n$. The notion of
fiberwise distributions $(\mu_{x})_{x\in\mathbb{R}^{n}}$ corresponds to that
of SAGS in Section III. On the other hand, $f_{2}$ is related to signal
transformation including filtering.

As category theory is out of the scope of the paper, details on the
categorical perspective can be found in [39].

##  Appendix C Remarks on piecewise linear functions

As we have seen in Example 3 that for $\mathfrak{p}=(\mathcal{A},f)$, if
$\mathcal{A}$ is locally constant, then $e_{\mathfrak{p}}$ is piecewise linear
a.e. Let $e_{\mathfrak{p}}$ consist of linear transformations
$(\mathbf{P}_{i})_{i\geq 1}$. To simplify the discussion, we assume that each
$\mathbf{P}_{i}$ is a projection to an $m$-dimensional subspace
$\mathcal{W}_{i}$ of $\mathbb{R}^{n}$ for $m<n$. We want to discuss sampling
and recovery with this setup. A rigorous discussion requires the theory of
Grassmannians to parametrize linear spaces[40], which is out of the scope. We
content to explain the main idea.

In classical GSP when there is only a single linear projection $\mathbf{P}$ to
($m$ dimensional) $\mathcal{W}$, then sampling and recovery amount to find a
set of $m$ coordinates corresponding to $V^{\prime}\subset V$, and identify
the intersection of $\mathcal{W}$ with the signal space $\mathcal{S}$ with
fixed observation on $V^{\prime}$.

This approach does not work for a set of projections $(\mathbf{P}_{i})_{i\geq
1}$ as above. Assume that there is an index $j$ such that
$\mathcal{W}_{j}\cap\mathcal{S}=\\{\mathbf{x}\\}$, which we want to identify.
The challenge is that $\mathcal{S}\cap\mathcal{W}_{i}$ is usually non-empty
for any $i\geq 1$. Therefore, it is not possible to find $\mathbf{x}$ based on
the partial observations at $V^{\prime}$.

However, the issue can be resolved by enlarging $V^{\prime}$ by including one
more sample. For the new $V^{\prime}$ and $\mathcal{S}$, it is usually true
that $\mathcal{S}\cap\mathcal{W}_{i}=\emptyset$ by dimension counting, except
for the single index $j$ that is known (a priori) to satisfy
$\mathcal{W}_{j}\cap\mathcal{S}=\\{\mathbf{x}\\}$.

Though we have been vague in the claims by using “usually” a few times, it is
possible to make them precise by stating “non-empty open set in the
Grassmannian manifold”. However, our key message here is that in almost any
case, it is necessary and sufficient to find $m+1$ samples.

## References

  * [1] D. I. Shuman, B. Ricaud, and P. Vandergheynst, “A windowed graph Fourier transform,” in _Proc. IEEE Workshop on Stats. Signal Process._ , 2012.
  * [2] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Process. Mag._ , vol. 30, no. 3, pp. 83–98, 2013.
  * [3] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [4] ——, “Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure,” _IEEE Signal Process. Mag._ , vol. 31, no. 5, pp. 80–90, 2014.
  * [5] A. Gadde, A. Anis, and A. Ortega, “Active semi-supervised learning using sampling theory for graph signals,” in _Proc. ACM SIGKDD_ , 2014.
  * [6] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning Laplacian matrix in smooth graph signal representations,” _IEEE Trans. Signal Process._ , vol. 64, no. 23, pp. 6160–6173, 2016.
  * [7] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _NeurIPS_ , 2016.
  * [8] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _ICLR_ , 2017.
  * [9] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data under Laplacian and structural constraints,” _IEEE J. Sel. Top. Signal Process._ , vol. 11, no. 6, pp. 825–841, 2017.
  * [10] F. Grassi, A. Loukas, N. Perraudin, and B. Ricaud, “A time-vertex signal processing framework: Scalable processing and meaningful representations for time-series on graphs,” _IEEE Trans. Signal Process._ , vol. 66, no. 3, pp. 817–829, 2018.
  * [11] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [12] B. Girault, A. Ortega, and S. S. Narayanan, “Irregularity-aware graph fourier transforms,” _IEEE Trans. Signal Process._ , vol. 66, no. 21, pp. 5746–5761, 2018.
  * [13] F. Ji and W. P. Tay, “A Hilbert space theory of generalized graph signal processing,” _IEEE Trans. Signal Process._ , vol. 67, no. 24, pp. 6188 – 6203, 2019.
  * [14] H. Huang, N. Sunar, and J. Swaminathan, “Do noisy customer reviews discourage platform sellers? empirical analysis of an online solar marketplace,” _SSRN Electronic Journal_ , 2020.
  * [15] S. Prabhakar and R. Cheng, _Data Uncertainty Management in Sensor Networks_. Springer, 2009.
  * [16] F. Ji, S. H. Lee, Z. Kai, W. P. Tay, and J. Yang, “Distributional signals for node classification in graph neural networks,” _openreview.net/forum?id=eoqfMQJogx0_ , 2023.
  * [17] C. Villani, _Optimal Transport, Old and New_. Springer, 2009.
  * [18] F. Ji, S. H. Lee, H. Meng, Z. Kai, W. P. Tay, and J. Yang, “Leveraging label non-uniformity for node classification in graph neural networks,” _openreview.net/forum?id=HfUWnPeLLH_ , 2023.
  * [19] F. Ji, W. P. Tay, and A. Ortega, “Graph signal processing over a probability space of shift operators,” _arXiv preprint arXiv:2108.09192v2_ , 2022.
  * [20] X. Jian, F. Ji, and W. P. Tay, “Generalizing graph signal processing: High dimensional spaces, models and structures,” _Found. Trends Signal Process._ , vol. 17, no. 3, pp. 209–290, 2023.
  * [21] W. Rudin, _Real and Complex Analysis_. McGraw-Hill, 1987.
  * [22] M. Puschel and J. Moura, “Algebraic signal processing theory: Foundation and 1-D time,” _IEEE Trans. Signal Process._ , vol. 56, no. 8, pp. 3572–3585, 2008.
  * [23] A. Agaskar and Y. M. Lu, “A spectral graph uncertainty principle,” _IEEE Trans. Inf. Theory_ , vol. 59, no. 7, pp. 4338–4356, 2013.
  * [24] S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete signal processing on graphs: Sampling theory,” _IEEE Trans. Signal Process._ , vol. 63, no. 24, pp. 6510–6523, 2015.
  * [25] M. Tsitsvero, S. Barbarossa, and P. Di Lorenzo, “Signals on graphs: Uncertainty principle and sampling,” _IEEE Trans. Signal Process._ , vol. 64, no. 18, pp. 4845–4860, 2016.
  * [26] A. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Sampling of graph signals with successive local aggregations,” _IEEE Trans. Signal Process._ , vol. 64, no. 7, pp. 1832–1843, 2016.
  * [27] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for bandlimited graph signals using graph spectral proxies,” _IEEE Trans. Signal Process._ , vol. 64, no. 14, pp. 3775–3789, 2016.
  * [28] F. Ji, G. Kahn, and W. P. Tay, “Signal processing on simplicial complexes with vertex signals,” _IEEE Access_ , vol. 10, pp. 41 889–41 901, 2022.
  * [29] Y. Tanaka, Y. Eldar, A. Ortega, and G. Cheung, “Sampling signals on graphs: From theory to applications,” _IEEE Signal Process. Mag._ , vol. 37, no. 6, pp. 14–30, 2020.
  * [30] L. Ruiz, L. F. O. Chamon, and A. Ribeiro, “Graphon signal processing,” _IEEE Trans. Signal Process._ , vol. 69, pp. 4961–4976, 2021.
  * [31] P. Lax, _Functional Analysis_ , 1st ed. Wiley-Interscience, 2002.
  * [32] A. Kolmogorov, _Foundations of the Theory of Probability_. New York: Chelsea., 1956.
  * [33] D. Kazhdan, “On the connection of the dual space of a group with the structure of its closed subgroups,” _Funct. Anal._ , vol. 1, no. 1, pp. 63–65, 1967\.
  * [34] R. Durrett, _Probability: Theory and Examples (5th ed.)_. Cambridge University Press, 2019.
  * [35] I. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” in _ICLR_ , 2015.
  * [36] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, “Towards deep learning models resistant to adversarial attacks,” in _ICLR_ , 2018.
  * [37] P. Billingsley, _Convergence of Probability Measures_. New York, NY: John Wiley and Sons, Inc., 1999.
  * [38] T. Hungerford, _Algebra (Graduate Texts in Mathematics) (v. 73)_ , 8th ed. Springer, 2002.
  * [39] F. Ji, X. Jian, and W. P. Tay, “Graph signal processing with categorical perspective,” 2023.
  * [40] J. Milnor and J. Stasheff, _Characteristic classes. Annals of Mathematics Studies. Vol. 76_. Princeton, NJ: Princeton University Press, 1974.