Graph Signal Processing over a Probability Space of Shift Operators
Feng Ji, Wee Peng Tay, Senior Member, IEEE and Antonio Ortega, Fellow, IEEE

F. Ji and W. P. Tay are with the School of Electrical and Electronic Engineering, Nanyang Technological University, 639798, Singapore (e-mail<EMAIL_ADDRESS>wptay@ntu.edu.sg). A. Ortega is with Viterbi School of Engineering, University of Southern California, Los Angeles, CA 90089-2564  (e-mail: aortega@usc.edu).

Graph signal processing (GSP) uses a shift operator to define a Fourier basis for the set of graph signals. The shift operator is often chosen to capture the graph topology. However, in many applications, the graph topology may be unknown a priori, its structure uncertain, or generated randomly from a predefined set for each observation. Each graph topology gives rise to a different shift operator. In this paper, we develop a GSP framework over a probability space of shift operators. We develop the corresponding notions of Fourier transform, convolution, and band-pass filters in this framework, which subsumes classical GSP theory as the special case where the probability space consists of a single shift operator. We show that a convolution filter under this framework is the expectation of random convolutional filters in classical GSP, while the notion of bandlimitedness requires additional wriggle room from being simply a fixed point of a band-pass filter. We develop a mechanism that facilitates mapping from one space of shift operators to another, which allows our framework to be applied to a rich set of scenarios. Numerical results on both synthetic and real datasets verify the superiority of performing GSP over a probability space of shift operators versus restricting to a single shift operator.

Graph signal processing, distribution of operators, Fourier transform, convolution, band-pass, sampling

§ INTRODUCTION

Since its emergence, the theory and applications of graph signal processing (GSP) have rapidly developed <cit.>. GSP theory is based on the choice of a graph shift operator or fundamental graph operator, which is a preferred linear transformation on the vector space of graph signals. Once such an operator is given, there is a systematic way to develop a framework for signal processing tasks. To highlight a few important elements of GSP, the change of basis wrt an eigenbasis of the shift operator defines the graph Fourier transform (GFT) <cit.>. The coefficients of a graph signal in the new basis are the components in the frequency domain. A central theme of GSP theory is the theory of filtering <cit.>, which discusses transformation families. Convolution is a transformation by a diagonal matrix in the frequency domain.
Sampling <cit.> refers to observing only a subset of vertices and reconstructing the signal under some assumptions on the type of signals to be considered, e.g., bandlimited signals.

Suppose $V$ is a sensor network (left). There are different possible connections among $V$, as shown on the right.

Many techniques have been developed in recent years to learn graphs from data <cit.>.
Roughly speaking, some methods consider features associated with each node, and derive graph edges from the distances between feature vectors, while others learn a graph such that a series of signal examples have some desirable properties, such as smoothness <cit.>.
Examples of the former type of methods include $k$ nearest neighbors ($k$-NN) and its variations <cit.>, while the latter include techniques based on smoothness <cit.> or inverse covariance <cit.>.
These methods often require the choice of hyperparameters (e.g., $k$ in $k$-NN illustrated in <ref>), which can be selected so that the performance of a downstream task is optimized.
Additionally, when the graph construction is based on several heterogeneous feature vectors, the conventional approach is to construct a single graph by essentially giving different weights to each of the features. A popular approach based on this principle is the bilateral filter and related methods <cit.>. In the machine learning community, a series of works (cf. <cit.>) is based on the idea of constructing a homogeneous graph by choosing only vertices connected by a specific sequence of edge types (known as a meta-path), before further inference or learning. In the case of information transmission over a network, the effective network topology depends not only on physical connections but also on transmission rates $\zeta$ between members of the network. Certain topology inference methods are sensitive to knowing accurately a priori such rates <cit.>.

Common to all the methods described above is the goal of designing a single graph, even though that may entail choosing a specific design parameter (e.g., a hyperparameter to learn a graph from data) or combining heterogeneous features (e.g., by choosing the relative weights for the different components in the features associated with a node). In this paper, we introduce a novel approach that avoids having to choose a single graph by developing processing on distributions of graphs. Thus, instead of selecting one hyperparameter, it is possible to work with multiple choices of hyperparameters simultaneously. Similarly, instead of combining heterogeneous features into a single graph, it is possible to proceed with multiple graphs in parallel, selecting instead the relative weights between the graphs for the downstream task.

In this paper, we consider a probability space of graph shift operators for signals on a finite vertex set $V$. The probability distribution corresponds to our prior belief of which shift operator is more likely, and can be estimated from available data. The distribution can be associated with hyperparameters such as those in the examples earlier. We develop a “distribution version” of the graph Fourier transform and the associated theory of filtering. Our main contributions are as follows:

* We introduce a novel GSP framework for a probability space $(\mathcal{X},\scA,\mu_{\mathcal{X}})$, where $\mathcal{X}$ is a sample space of shift operators, $\scA$ is a $\sigma$-algebra on $\mathcal{X}$, and $\mu_{\mathcal{X}}$ is a probability measure. We define the Fourier transform wrt this space. This new framework subsumes classical GSP theory as the special case where the probability space consists of a single shift operator.

* We develop the concept of a convolution filter in our framework and show that this is an expectation of convolution filters in classical GSP theory. In the case where $\mathcal{X}$ can be parametrized by a real parameter, a convolution filter can be characterized as a bi-polynomial filter, which is a polynomial filter whose coefficients are themselves polynomials of the parametrization variable of $\mathcal{X}$.

* We develop the notion of bandlimited signals, which in general, do not form a vector space. We develop bounds for recovering such signals from a subset of vertex signals.

* As filters and observed signals may not always be associated with the same probability space of operators, we develop a mechanism that allows us to map from one probability space of operators to another.

The main focus of the paper is on the theory of signal processing with a given distribution of operators on a network. However, learning the distribution of operators is itself an important topic if such prior knowledge is unavailable. For completeness, we give an overview of a Bayesian approach based on Markov chain Monte Carlo (MCMC) <cit.> in <ref>. On the other hand, in the machine learning community, there are end-to-end Bayesian learning approaches <cit.> to address ambiguity in graph structure as well as robustness in graph neural networks.

The rest of this paper is organized as follows. We introduce the basic setup and define the Fourier transform in <ref>. In <ref> and <ref>, we discuss various families of filters. We present sampling theory in <ref> in conjunction with band-pass filters. In <ref>, we discuss base changes that deal with switching from one sample space of operators to another. The framework requires knowledge of the distribution of the shift operators. In <ref>, we describe ways to learn a distribution if such a priori knowledge is unavailable. We present several numerical examples in <ref> and conclude in <ref>. Proofs of all results are deferred to <ref>.

A preliminary version of this work was presented in <cit.>. In this paper, we include more thorough theoretical results and discussions not in the preliminary version. We also include further numerical experiments.

Notations: We use $\circ$ to denote function composition. Let $\Real$ denote the set of real numbers, $\Real_+$ the set of non-negative real numbers, $M_n(\Real)$ be the space of $n\times n$ real matrices, and $[n]$ the
discrete set $\{1,2,\ldots, n\}$. $\E_{\mu}$ is the expectation operator wrt the probability measure $\mu$.
For a Hilbert space $L^2(\Omega)$, we denote its inner product and norm (and associated operator norm) as $\ip{\cdot}{\cdot}_{L^2(\Omega)}$ and $\norm{\cdot}_{L^2(\Omega)}$, respectively. If the Hilbert space is clear from the context, to avoid clutter, we drop the subscripts and use $\ip{\cdot}{\cdot}$ and $\norm{\cdot}$ respectively. We use calligraphic fonts such as $\mathcal{X},\mathcal{Y}, \mathcal{Z}$ for spaces of operators, while the operators are boldfaced.

§ DISTRIBUTION OF SHIFTS AND THE FOURIER TRANSFORM

In this section, we describe the setup of our framework and introduce the Fourier transform and its left inverse within the new setup. We demonstrate with visual examples the spectral plots of signals in our framework.

Let $V$ be the set of vertices of a finite graph $G$, where $|V|=n$. A signal on $V$ is a function $f: V \mapsto \mathbb{R}$, where each signal $f(v)$ associates a real value to a vertex $v\in V$. Denote the Hilbert space of such signals by $L^2(V)$, with $\langle f,f' \rangle = \sum_{v\in V}f(v)f'(v)$. It can be identified with $\mathbb{R}^n$ for a fixed ordering of $V$.

Suppose $\mathcal{X}$ is a metric space of operators on $L^2(V)$, each of which has eigenvectors that form an orthonormal eigenbasis of $\Real^n$, e.g., each $\bX\in \mathcal{X}$ can be represented as an $n \times n$ symmetric matrix. Let $(\mathcal{X}, \scA, \mu_{\mathcal{X}})$ be a probability space. It may be abbreviated as $(\mathcal{X},\mu_{\mathcal{X}})$ if the $\sigma$-algebra $\scA$ is clear from the context.  We call $(\mathcal{X},\mu_{\mathcal{X}})$ the base space. Let $\mathcal{Y}$ be a set and consider the product $\mathcal{X}\times \mathcal{Y}$. Then $\{\bX\}\times \mathcal{Y}$, where $\bX\in \mathcal{X}$, is called the fiber at $\bX$. In some examples, for ease of presentation, we may abuse terminology by letting $\mathcal{X}$ to be a set of objects (e.g., graphs) instead of operators, where each of these objects is associated with a shift operator. Moreover, applying $\bX \in \mathcal{X}$ to a signal $f$ is denoted by $\bX(f)$.

To illustrate the above setup, suppose that the underlying graph $G$ is random and is generated by a distribution of graphs on $V$. Each $G$ drawn from the distribution gives a Laplacian $\bL_G$, and the collection of such $\bL_G$ yields the space $\mathcal{X}$. The probability measure $\mu_{\mathcal{X}}$ is thus induced by the distribution generating $G$. A specific case is when the edges among vertices in $V$ are fixed, and the probability distribution $\mu_{\mathcal{X}}$ of the graph Laplacian is induced by a probability distribution on the edge weights. If a training set is available, the probability measure $\mu_{\calX}$ can be learned via a Bayesian framework, see <ref>.

We next introduce the Fourier transform and its (left) inverse. Our definition is simply the collection of the classical GFT for each shift operator $\bX\in\calX$. This definition however leads to a nontrivial generalization of convolution and band-pass filters in the subsequent sections.

For each operator $\bX \in \mathcal{X}$, let $\lambda_{\bX,i}$ be the $i$-th eigenvalue of the operator $\bX$ (ordered increasingly in absolute value) and $u_{\bX,i}$ be an associated eigenvector with $\set{u_{\bX,i} \given i\in[n]}$ forming an eigenbasis of $\Real^n$. Let $L^2(\mathcal{X}\times [n])$ be the $L^2$ Hilbert space endowed with the product measure $\mu_{\mathcal{X}} \times |\cdot|$ where $|\cdot|$ is the counting measure, and the product $\sigma$-algebra $\scA\times 2^{[n]}$ where $2^{[n]}$ is the power set of $[n]$. Define the Fourier transform wrt $(\mathcal{X},\mu_{\mathcal{X}})$,

\begin{align}
\calF_{\mathcal{X}}: L^2(V) \mapsto L^2(\mathcal{X}\times [n])
\end{align}

by $\calF_{\mathcal{X}}(f) = \hat{f}$ where $\hat{f}(\bX,i) = \ip{f}{u_{\bX,i}}$.

The space $L^2(V)$ is an $n$-dimensional vector space. On the other hand, $L^2(\mathcal{X}\times [n])$ can be infinite dimensional in general. Therefore, it is impossible for $\calF_\mathcal{X}$ to be invertible. However, it has a left inverse:

\begin{align} \label{eq:psx}
\calF^{\dagger}_{\mathcal{X}}: L^2(\mathcal{X}\times [n]) & \mapsto L^2(V), \\
g & \mapsto \int_{\mathcal{X}} \sum_{i=1}^n g(\bX,i)u_{\bX,i} \ud \mu_{\mathcal{X}}(\bX).
\end{align}

Note that $\set{u_{\bX,i}\given (\bX,i)\in \mathcal{X}\times [n]}$ are the kernels of the integration. However, unlike GSP, they are not pair-wise orthogonal for different $\bX$. Intuitively, the Fourier transform should contain spectral information w.r.t. all the operators $\bX$ in the family $\mathcal{X}$. Therefore, in <ref>, we enlarge the “frequency domain” as compared with classical GSP, and hence the Hilbert space $L^2(\mathcal{X}\times [n])$ is used as the codomain of the Fourier transform. To ensure recovery of the original signal, we need to re-group the signals in the spectral domain by taking a “weighted average” (more precisely, an integral), and this leads to $\calF^{\dagger}_\mathcal{X}$.

Regarding $\calF_\mathcal{X}$ and $\calF^{\dagger}_\mathcal{X}$, the following observations can be verified directly from the definitions and we omit the proof.

* $\calF_{\mathcal{X}}$ and $\calF^{\dagger}_\mathcal{X}$ are both well-defined.

* (Left inverse) $\calF^{\dagger}_{\mathcal{X}} \circ \calF_\mathcal{X}$ is the identity map on $L^2(V)$.

* (Parseval's identity) $\norm{f}_{L^2(V)} = \norm{\calF_{\mathcal{X}}(f)}_{L^2(\mathcal{X}\times[n])}$ for each graph signal $f\in L^2(V)$.

Another interpretation of Parseval's identity is as follows. The set of vectors $\mathfrak{F}_{\mathcal{X}} = \{u_{\bX,i} : \bX\in \mathcal{X}, i\in [n] \} \subset L^2(V)$ forms a (continuous) tight frame of the finite dimensional vector space $L^2(V)$ (cf. <cit.>). With this point of view, the operator $\calF_\mathcal{X}$ can be interpreted as an analysis operator, while the operator $\calF^{\dagger}_\mathcal{X}$ is a synthesis operator.

We want to further analyze the left inverse $\calF^{\dagger}_{\mathcal{X}}$ for use in subsequent sections. We observe that $\calF^{\dagger}_{\mathcal{X}}$ can be decomposed as follows:

\begin{align}
\calF^{\dagger}_{\mathcal{X}}: L^2(\mathcal{X}\times [n]) \xrightarrow{\alpha_{\mathcal{X}}} L^2(\mathcal{X}\times V) \xrightarrow{\beta_{\mathcal{X}}} L^2(V), \label{eq:invdecomposition}
\end{align}

where here for each $\bX\in \mathcal{X}$, $v\in V$, and $g\in L^2(\mathcal{X}\times [n])$,

\begin{align}
\alpha_{\mathcal{X}}(g)(\bX,v) &= \sum_{i=1}^n g(\bX,i) u_{\bX,i}(v), \label{alphaX}
\intertext{$u_{\bX,i}(v)$ is the $v$-component of the eigenvector $u_{\bX,i}$, and for $q\in L^2(\mathcal{X}\times V)$,}
\beta_{\mathcal{X}}(q)(v) &= \int_{\mathcal{X}} q(\bX,v) \ud\mu_{\mathcal{X}}(\bX). \label{betaX}
\end{align}

Suppose we let $q(\bX,v)=\alpha_{\mathcal{X}}(g)(\bX,v)$. Then $q(\bX,v)$ is the $v$-component of the inverse GFT wrt the operator $\bX$, and further applying $\beta_\mathcal{X}$ to $q$ yields $\beta_{\mathcal{X}}(q)(\cdot) = \E_{\mu_{\mathcal{X}}}[q(\bX,\cdot)]$, where the expectation is over $\bX\in \mathcal{X}$. The map $\alpha_\mathcal{X}$ is invertible with its inverse being the fiberwise GFT for each $\bX\in \mathcal{X}$.

Note that if $g= \calF_{\mathcal{X}}(f)$ for some $f\in L^2(V)$, then in the above discussion, $q(\bX,v)=f(v)$ for all $\bX\in \mathcal{X}$. Hence, $\beta_{\mathcal{X}}(q)(v)$ gives back $f(v)$. The function $\beta_\mathcal{X}$ and hence the above decomposition of $\calF^{\dagger}_\mathcal{X}$ will only output a signal different from $f$ (the one we started with), when we introduce two filter families in subsequent sections:

* the convolution filter family that “inserts” a transformation on $L^2(\mathcal{X}\times [n])$ before applying $\alpha_\mathcal{X}$; and

* a kind of base change family that “inserts”, between $\alpha_\mathcal{X}$ and $\beta_\mathcal{X}$, a transformation $L^2(\mathcal{X}\times V)\mapsto L^2(\mathcal{Z}\times V)$, with $\mathcal{Z}$ also a probability space.

We end this section with some examples.

* Suppose $\mu_{\mathcal{X}} = \delta_{\bX_0}$ is a distribution concentrated at a single $\bX_0\in \mathcal{X}$. Then $\calF_\mathcal{X}$ is simply the GFT wrt the shift operator $\bX_0$ in classical GSP.

We consider a special case of <cit.>: Graph Slicing, as an example of heterogeneous graphs. Suppose the graph $G$ on the vertex set $V$ is fixed. In addition, we have Laplacians $\bL_0$ and $\bL_1$ associated with two subgraphs $G_1, G_2$ of $G$ such that $G_1$ and $G_2$ have disjoint edge sets and $G=G_1\cup G_2$. We form $\mathcal{X}$ parametrized by the unit interval $T = [0,1]$: for each $t \in [0,1]$, let $\bL_t = (1-t)\bL_0+t\bL_1$ be an element of $\calX$. On $T$, $\scA$ is the Borel $\sigma$-algebra and $\mu_{\calX}$ is a probability measure on $T$. We shall use the following setup in subsequent sections. Let $G$ be a $2$D-lattice (e.g., corresponding to the graph of an image). Then $\bL_0$, $\bL_1$ are the Laplacians of the subgraph containing horizontal and vertical edges only, respectively. Intuitively, we are thinking of vertical and horizontal edges having different importance for signal processing purposes.

§ CONVOLUTION FILTERS

In this section, we introduce the family of convolution filters and some important subfamilies. We show that a convolution filter in our framework is an expectation (in the sense of a Bochner integral) of convolution filters in the classical GSP theory. We then show that under some technical conditions, every convolution filter is a bi-polynomial filter, which is a polynomial filter whose coefficients are themselves polynomials of another parameter.

Similar to classical GSP, convolutions are defined using multiplication of functions on the frequency domain. Given $\Gamma \in L^2(\mathcal{X}\times [n])$, multiplication by $\Gamma$ induces a mapping $L^2(\mathcal{X}\times [n]) \mapsto L^2(\mathcal{X}\times [n])$, which we also denote by $\Gamma$. For any $g\in L^2(\mathcal{X}\times[n])$, we define $\Gamma g(\bX,i)=\Gamma(\bX,i)g(\bX,i)$ for all $(\bX,i)\in \mathcal{X}\times[n]$.

Let $\Gamma \in L^2(\mathcal{X}\times [n])$. A convolution filter $\bigstar_{\Gamma}: L^2(V) \mapsto L^2(V)$ is defined by the composition $\calF^{\dagger}_{\mathcal{X}}\circ \Gamma\circ\calF_\mathcal{X}$, i.e., for $f\in L^2(V)$ and any $v\in V$,

\begin{align}\label{def:conv}
\bigstar_{\Gamma}(f)(v) = \int_{\mathcal{X}} \sum_{i=1}^n \Gamma(\bX,i)\ip{f}{u_{\bX,i}} u_{\bX,i}(v) \ud\mu_{\mathcal{X}}(\bX).
\end{align}

To further understand the expression, for each $\bX\in \mathcal{X}$, we use $\Gamma_{\bX}$ to denote the function in $L^2([n])$ via the formula

\begin{align}\label{eq:gg}
\Gamma_{\bX}(i) = \Gamma(\bX,i).
\end{align}

Each $\bigstar_{\Gamma_{\bX}}$ is then a convolution filter on $L^2(V)$ in classical GSP theory. We call this a fiberwise convolution. This filter is nothing but the composition $\alpha_{\mathcal{X}}\circ \Gamma \circ \calF_\mathcal{X}$ evaluated at $(\bX,\cdot)$ (cf. <ref>).

We next show that a convolution filter can be expressed as an expectation of fiberwise convolutions.

For any $\Gamma \in L^2(\mathcal{X}\times [n])$, $\bigstar_{\Gamma} = \E_{\mu_{\mathcal{X}}}[\bigstar_{\Gamma_{\bX}}]$ and it is a bounded linear operator.

Note that the expectation of operators here is in the sense of a Bochner integral <cit.>. Notation-wise, an expectation of operators remains an operator.

In the following examples, we discuss special cases of convolution filters, to highlight some differences between our framework and classical GSP.

*  For any signal $g \in L^2(V)$, its Fourier transform $\hat{g}=\calF_{\mathcal{X}}(g)$ belongs to $L^2(\mathcal{X}\times [n])$. It induces a convolution filter $\bigstar_{\hat{g}} : L^2(V) \mapsto L^2(V)$ that maps $f$ to $h = \calF^{\dagger}_{\mathcal{X}}(\calF_{\mathcal{X}}(g)\cdot\calF_{\mathcal{X}}(f))$, i.e.,

\begin{align*}
h(v) = \int_{\mathcal{X}} \sum_{i=1}^n \ip{g}{u_{\bX,i}}\ip{f}{u_{\bX,i}} u_{\bX,i}(v) \ud\mu_{\mathcal{X}}(\bX),
\end{align*}

for each $v\in V$. A general convolution filter is associated with an element in $L^2(\mathcal{X}\times [n])$. And not every convolution filter is obtained from an element of $L^2(V)$ as in this example. On the contrary, in classical GSP where $\mathcal{X}$ is a singleton, every convolution operator corresponds to a graph signal, as described in this example.

*  If there is a uniform upper bound on the operator norm of $\bX\in \mathcal{X}$, then $\Lambda: (\bX,i) \mapsto \lambda_{\bX,i}$ belongs to $L^2(\mathcal{X}\times [n])$. From <ref>, we obtain

\begin{align*}
\bigstar_{\Lambda}(f) & = \int_{\mathcal{X}} \sum_{i=1}^n \lambda_{\bX,i} \ip{f}{u_{\bX,i}} u_{\bX,i} \ud\mu_{\mathcal{X}} \\
& = \int_{\mathcal{X}} \bX(f) \ud\mu_{\mathcal{X}} = \E_{\mu_{\mathcal{X}}}[\bX](f).
\end{align*}

Consequently, $\bigstar_{\Lambda} = \E_{\mu_{\mathcal{X}}}[\bX]$, the expectation of operators in $\mathcal{X}$. More generally, $\bigstar_{\Lambda^k} = \E_{\mu_{\mathcal{X}}}[\bX^k]$ for positive integers $k$. As $\calF_{\mathcal{X}}\circ \calF^{\dagger}_\mathcal{X}$ is not the identity map, ${\bigstar_{\Lambda^k}} \neq {(\bigstar_{\Lambda})^k}$. This is different from classical GSP where $\mu_{\mathcal{X}}$ is concentrated at a single shift operator $\bX_0$.

More concretely, following <ref><ref>, suppose parameters $\bX$ follows the uniform distribution on $[0,1]$. It is straight forward to compute that $({\bigstar_{\Lambda}})^2 = (\bL_0^2+\bL_1^2 + \bL_0\bL_1+\bL_1\bL_0)/4$ and ${\bigstar_{\Lambda^2}} = (2\bL_0^2+2\bL_1^2+\bL_0\bL_1+\bL_1\bL_0)/6$. In general, not only $({\bigstar_{\Lambda}})^2\neq \bigstar_{\Lambda^2}$ are distinct, they do not have common eigenbasis, i.e., they are not shift invariant with each other.

Through <ref>, we see differences between our framework and classical GSP on convolution filters. However, there are also common phenomena between them. In classical GSP, under favorable conditions, a convolution filter is always a polynomial of the shift operator. This observation is important as it facilitates computation with convolution filters. We next discuss an analogy in our case.

If almost surely every $\bX\in \mathcal{X}$ does not have repeated eigenvalues, any convolution filter has the form $\E_{\mu_{\calX}}[R(\bX)]$, where $R$ is a mapping $\mathcal{X} \mapsto M_n(\mathbb{R})$ such that $R(\bX)$ is a degree $n-1$ polynomial in $\bX$.

Under the conditions of <ref>, for almost surely every $\bX\in \mathcal{X}$, the fiberwise convolution $\bigstar_{\Gamma_{\bX}}$ of $\bigstar_{\Gamma}$ is nothing but a polynomial $R(\bX)= \sum_{0\leq i\leq n-1}a_i(\bX)\bX^i$ of degree at most $n-1$ in $\bX$. On the other hand, for each fixed degree $0\leq i\leq n-1$, we may look at the coefficient $a_i(\bX)$ of the $i$-th monomial for each $R(\bX)$, which gives rise to a function on $\bX\in \mathcal{X}$. Motivated by the above discussions, we now consider the following subspace of convolution filters, called bi-polynomial filters. Each bi-polynomial filter is a convolution filter, and we are also interested in conditions that ensure a convolution filter is bi-polynomial.

Suppose $\mathcal{X}$ is parametrized by $T\subset \mathbb{R}$ via a homeomorphism $t \in T\mapsto \bX_t\in \mathcal{X}$. Let the measure induced by $\mu_{\mathcal{X}}$ on $T$ be $\mu_T$. A convolution filter $\bigstar_{\Gamma}$ with $\Gamma \in L^2(\mathcal{X}\times [n])$ is called a bi-polynomial filter on $\mathcal{X}$ if for each $\bX_t\in \mathcal{X}$, there are polynomials $a_i(t)$, $0\leq i \leq k$, with degrees (in $t$) bounded by $d$ such that $\bigstar_{\Gamma_{\bX_t}} = \sum_{0\leq i\leq k} a_i(t) \bX_t^i$. We say that its bi-degree is bounded by $(d,k)$.

For the rest of this section, we assume the existence of such a parameter space $T$ as in <ref>. To give some simple examples, in the $k$-NN construction, the parameter space $T$ for the space $\mathcal{X}$ of graph shifts associated with different values of $k$ can be chosen as the discrete set $\{1,\ldots, n-1\}$, with $n$ being the number of nodes. In <ref><ref>, the parameter space is $T=[0,1]$ as we are taking convex combinations of two given graph shift operators. By <ref> below, every convolution filter is a bi-polynomial filter. If such a filter $\bF_t$ has its bi-degree bounded by $(1,1)$, then it takes the explicit form
\bF_t = (a_0 + a_1t)\big(t\bL_1+(1-t)\bL_0)\big) + (b_0 + b_1t)\bI, t\in [0,1],
where $\bL_0$ and $\bL_1$ are defined in <ref><ref>, $\bI$ is the identity transform and $a_0,a_1,b_0,b_1$ are real coefficients.

Suppose $\calX$ is parametrized by $T\subset\Real$, a finite set or bounded interval. If almost surely every $\bX_t \in \mathcal{X}, t\in T$ has no repeated eigenvalues and has uniformly bounded operator norm, then every convolution filter is a bi-polynomial filter.

§ BAND-PASS FILTERS AND SAMPLING

In this section, we develop the notion of band-pass filters, which is a special family of convolution filters. We then introduce the concept of bandlimited signals in our framework and discuss their sampling results.

We start by introducing band-pass filters together with the notion of bandlimited signals. Suppose $\mathcal{Y} \subset \mathcal{X}\times [n]$ is a measurable subset. Recall that the indicator function ${\bf 1}_{\mathcal{Y}}$ on $\mathcal{Y}$ is defined as ${\bf 1}_{\mathcal{Y}}(\bY) = 1$ if $\bY\in \mathcal{Y}$ and $0$ otherwise.

For a measurable subset $\mathcal{Y} \subset X\times [n]$, the band-pass filter $\bB_{\mathcal{Y}}$ w.r.t. $\mathcal{Y}$ is defined as the convolution filter associated with ${\bf 1}_{\mathcal{Y}} \in L^2(\mathcal{X}\times [n])$.[Note that $\bB_{\calY}=\bigstar_{{\bf 1}_{\calY}}$. To simplify notations, we use $\bB_{\calY}$ here instead.] For $\epsilon\geq 0$, the set of $(\mathcal{Y},\epsilon)$-bandlimited signals consists of graph signals $f\in L^2(V)$ such that $\norm{\bB_{\mathcal{Y}}(f)-f}\leq \epsilon$.

It is important to note that the filter $\bB_{\mathcal{Y}}$ is not a projection in general as it is the expectation of fiberwise projections (cf. <ref>). This means that $\bB_{\mathcal{Y}}$ may not even have non-zero fixed points $f= \bB_{\mathcal{Y}}(f)$. Therefore, we are not able to define bandlimited signals as the space of fixed points of a band-pass filter as in classical GSP. As a consequence of <ref>, the set of $(\mathcal{Y},\epsilon)$-bandlimited signals may not be a vector space. However, if $\mu_{\mathcal{X}}$ concentrates on a single $\bX_0 \in \mathcal{X}$ and $\epsilon = 0$ in <ref>, we recover the theory of band-pass filters and bandlimited signals in classical GSP.

Having introduced band-pass filters and bandlimited signals, we now discuss sampling theory that studies the recoverability of bandlimited signals from partial signal observations. To start, we have the following basic observation.

The set of $(\mathcal{Y},\epsilon)$-bandlimited signals is convex. Moreover, if $\epsilon>0$, then the signal $f(v) \equiv 0$, for all $v\in V$, is an interior point.

As a consequence of <ref>, if $\epsilon>0$ and $V'$ is a proper subset of $V$, then the signal values at $V'$ of a $(\mathcal{Y},\epsilon)$-bandlimited signal $f$ do not uniquely determine $f$. To see this, <ref> implies that the intersection of $(\mathcal{Y},\epsilon)$-bandlimited signals and the set of vectors with their $V'$ components fixed is an open subset of the latter. In particular, such an intersection is either empty or contains more than one element. Therefore, for signal recovery from sub-samples in the case $\epsilon>0$, we can only aim for approximations instead of exact recovery.

All the eigenvalues of $\bB_{\mathcal{Y}}$ are contained within the closed interval $[0,1]$.

Let $0\leq \lambda_1 \leq \ldots \leq \lambda_n\leq 1$ be the eigenvalues of $\bB_{\mathcal{Y}}$ and $u_1,\ldots, u_n$ be the associated eigenvectors chosen to form an orthonormal eigenbasis, since $\bB_{\mathcal{Y}}$ can be represented as an $n\times n$ symmetric matrix. In general, $\lambda_j$, $1\leq j\leq n$ can be distinct, while they can only be either $0$ or $1$ when $\mu_{\mathcal{X}}$ concentrates on a singleton.

Suppose $f = \sum_{1\leq i\leq n}a_iu_i$ is a $(\mathcal{Y},\epsilon)$-bandlimited signal. Then, for $1\leq j \leq n$ such that $\lambda_j\neq 1$, we have $\sum_{1\leq i\leq j}a_i^2\leq \epsilon^2/(1-\lambda_j)^2$.

As a consequence, if $\lambda_j$ is close to $0$ for some $1\leq j\leq n$, then the components of a signal spanned by $u_1, \ldots, u_j$ have a small contribution. Therefore, if one is restricted to sampling from a subset of $n-j$ vertices, one should choose vertices whose signals are spanned by $u_{j+1},\ldots, u_n$.

Consider the subspace of signals spanned by $u_{j+1},\ldots, u_n$. Let $\bU_{>j}=[u_{j+1},\ldots, u_n]$ be the matrix with columns $u_{j+1},\ldots, u_n$. Each of its rows corresponds to a node in $V$. A subset $V_j \subset V$ of size $n-j$ is called a uniqueness set <cit.> if the submatrix of $\bU_{>j}$ consisting of the rows of $\bU_{>j}$, corresponding to vertices in $V_j$, is invertible. Let this submatrix be $\bG_{V_j}$, the recovery matrix.  If $f$ belongs to the span of $u_{j+1},\ldots, u_n$ (cf. bandlimitedness in classical GSP), then $\bU_{>j}\bG_{V_j}^{-1}(f_{V_j})$ recovers the graph signal $f$ perfectly. Denote the operator norm of $\bG_{V_j}^{-1}$ by $\sigma_{V_j}$.

Suppose observation of a $(\mathcal{Y},\epsilon)$-bandlimited signal $f\in L^2(V)$ is made at a uniqueness set $V_j$, denoted by $f_{V_j}$. Let $f'$ be the linear combination of $u_{j+1},\ldots,u_n$ with coefficients the entries of $\bG_{V_j}^{-1}(f_{V_j})$, i.e., $\bU_{>j}\bG_{V_j}^{-1}(f_{V_j})$. Then:

*  $\norm{f'-f} \leq \epsilon(1+\sigma_{V_j})/(1-\lambda_j)$.

* $f'$ is $(\mathcal{Y},\epsilon')$-bandlimited with

\begin{align*}
\epsilon'= \epsilon\parens*{1+2\frac{1+\sigma_{V_j}}{1-\lambda_j}}.
\end{align*}

From the above discussions, we see that $\lambda_j$ is an important quantity that controls how well we can recover a bandlimited signal, with a smaller $\lambda_j$ leading to better recovery.

§ BASE CHANGE

So far, we have been dealing with convolution filters exclusively. <ref> summarizes an important observation regarding convolution filters: each is an expectation of a “random variable of fiberwise convolutions” on $\mathcal{X}$. On the other hand, once we have such random variables, we may pass them from one probability space to another. This yields the base change filters, which we discuss in this section. We first conceptualize base change in an abstract fashion, and then in special cases, describe explicit formulas with illustrated examples. We start with the basic setup for base change.

Suppose $(\mathcal{Z},\mu_{\calZ})$ and $(\mathcal{X},\mu_{\mathcal{X}})$ are probability spaces of operators.
In some applications, it is more natural to define a filter family and its distribution on another space $\mathcal{Z}$ while the graph signal is associated with the base space $\mathcal{X}$, or vice versa. A measurable function
h: \mathcal{Z} \mapsto \mathcal{X}
allows pushforward to $\mathcal{X}$ or pullback to $\mathcal{Z}$, on which further signal processing is then performed.

We motivate the need for a base change framework with the following examples.

Consider infection propagation over a network <cit.> in which a piece of information or disease from a source $s\in V$, propagates on a network represented as a graph. We assume that once $v\in V$ is infected by a neighbor, it remains in this status, and the infection spreads from one node to a neighbor according to a known stochastic model $\zeta$. An infection propagation then forms a random tree rooted at a source node. The graph signal here is the infection status of every node in the graph at some observation time. The adjacency matrices of such infection trees form the space of shift operators $\calZ$ that we can utilize to perform inference on graph signals generated by the model $\zeta$. Now suppose the infection model changes to $\zeta'$ and the tree adjacency matrices under this model form the space $\calX$. Suppose we can model the relationship between the two models using a mapping $h : \calZ \mapsto \calX$. For example, the spreading rate across some edges may have changed from the models $\zeta$ to $\zeta'$. For a concrete example, see <ref>. Now, we are given a training set of graph signals generated by model $\zeta$ from which we can learn $\mu_{\calZ}$ (cf. <ref>). Are we able to say anything about a corresponding measure $\mu_{\calX}$ for $\calX$? This will involve a pushforward of the learned measure $\mu_{\calZ}$ to the measurable space $\calX$.

A fisheye camera has an ultra wide-angle lens, but the image captured by it is distorted. Suppose we have a family of convolution filters (polynomials of shift operators in a space $\calZ$, where each shift operator captures certain relationships between pixels in the 2D image) with an associated distribution designed for normal 2D images. We would like to apply these filters to a fisheye image using the framework in <ref>. While rectification techniques <cit.> can recover a normal 2D image from a fisheye image, artifacts may occur. Furthermore, this pre-processing step leads to additional computational costs. Another approach is to pushforward the family of convolution filters to the space $\calX$ that consists of shifts that capture relationships between pixels in a fisheye image. For ease of illustration, we provide a toy example in <ref><ref> below involving the horizontal stretching of an image.

Let $F_{\mathcal{Z}}: \mathcal{Z} \mapsto M_n(\mathbb{R})$ (resp. $F_{\mathcal{X}}: \mathcal{X} \mapsto M_n(\mathbb{R})$) denote a mapping of each $\bZ\in \mathcal{Z}$ (resp. $\bX\in \mathcal{X}$) to a linear transformation ($n\times n$ matrix). We thus view $F_{\calZ}$ as a family of filters. For example, in the previous sections (cf. $\bigstar_{\Gamma_{\bX}}$ before <ref>), we have considered convolution filters in which each $F_{\mathcal{X}}(\bX)$ for $\bX \in \mathcal{X}$, is a convolution wrt the shift operator $\bX$.

We first introduce the pushforward of a measure. In some cases, we may start with a filter family $F_{\calX}$ and would like to have a corresponding family of filters on $\calZ$. This is the pullback of a filter family. These notions are defined formally as follows:

* Pushforward of measure: $h$ and $\mu_{\mathcal{Z}}$ induce a measure $h_*(\mu_{\mathcal{Z}})$ on $\mathcal{X}$, defined as $h_*(\mu_{\mathcal{Z}})(\mathcal{X}') = \mu_{\mathcal{Z}}(h^{-1}(\mathcal{X}'))$ for any measurable subset $\mathcal{X}' \subset \mathcal{X}$.

* Pullback of filter family: $h$ and $F_\mathcal{X}$ induce a family of filters $h^*(F_{\mathcal{X}})$ on $\mathcal{Z}$, defined by $h^*(F_{\mathcal{X}})(\bZ) = F_{\mathcal{X}}(h(\bZ))$.

Similarly, we can define the pullback of measure and pushforward of filters. Here, we need an additional assumption that there is a fiberwise measure. More specifically, for each $\bX\in \mathcal{X}$, there is a probability measure $\mu_{h^{-1}(\bX)}$ on the set $h^{-1}(\bX)=\set{\bZ\in \mathcal{Z}\given h(\bZ)=\bX}$. For example, one may choose $\mu_{h^{-1}(\bX)}$ to be the uniform distribution on $h^{-1}(\bX)$. In addition, we assume the technical condition that $\mathcal{Z}$ is locally compact and Hausdorff <cit.>, which is the case for most of the spaces we are interested in.

*  Pullback of measure: $h$, $\mu_{\mathcal{X}}$ and $\{\mu_{h^{-1}(\bX)}\}_{\bX\in \mathcal{X}}$ induce a measure $h^*(\mu_{\mathcal{X}})$ on $\mathcal{Z}$ defined by the integration formula

\begin{align*}
\int_{\mathcal{Z}} g \ud h^*(\mu_{\mathcal{X}}) = \int_{\mathcal{X}} \int_{h^{-1}(\bX)} g \ud\mu_{h^{-1}(\bX)} \ud\mu_{\mathcal{X}},
\end{align*}

where $g$ is any compactly supported continuous function on $\mathcal{Z}$. The measure $h^*(\mu_{\mathcal{X}})$ is uniquely determined by the Riesz–Markov–Kakutani representation theorem <cit.>.

*  Pushfoward of filter family: $h$, $F_{\mathcal{Z}}$ and $\{\mu_{h^{-1}(\bX)}\}_{\bX\in \mathcal{X}}$ induce a family of filters $h_*(F_{\mathcal{Z}})$ on $\mathcal{X}$ defined by

\begin{align*}
h_*(F_{\mathcal{Z}})(\bX) = \int_{h^{-1}(\bX)} F_{\mathcal{Z}}(\bZ) \ud\mu_{h^{-1}(\bX)}(\bZ),\ \bX\in \mathcal{X}.
\end{align*}

We summarize the essential ingredients of the different types of based changes in fig:bchange.

A summary of the different types of base changes.

For the rest of this section, we discuss the base change of a convolution filter more concretely and provide explicit formulas. Recall from <ref> that we have the following decomposition of the identity transform on $\mathcal{X}$:

\begin{align*}
\bI = \calF^{\dagger}_{\mathcal{X}}\circ \calF_{\mathcal{X}} = \beta_{\mathcal{X}}\circ \alpha_{\mathcal{X}} \circ \calF_{\mathcal{X}}.
\end{align*}

We shall insert base changes in this decomposition.

As we have seen, a convolution filter is constructed from $\Gamma \in L^2(\mathcal{X}\times [n])$. The map $h$ pulls it back to a function $h^*(\Gamma) \in L^2(\mathcal{Z}\times [n])$ defined by $h^*(\Gamma)(\bZ,i) = \Gamma(h(\bZ),i)$.  Let $h^\#: L^2(\mathcal{X}\times V) \mapsto L^2(\mathcal{Z}\times V)$ be $h^\#(q)(\bZ,v) = q(h(\bZ),v)$. We have the following associated base change convolutions.

For $\Gamma \in L^2(\mathcal{X}\times [n])$, the filter $\bF_{h^*(\Gamma)}$ defined as the composition

\begin{align}
\bF_{h^*(\Gamma)} = \calF^{\dagger}_{\mathcal{Z}}\circ h^*(\Gamma)\circ \calF_{\mathcal{Z}}: L^2(V) \mapsto L^2(V),
\end{align}

is a pullback of the filter $\Gamma$. The filter $\bF_{h^\#,\Gamma}$ is defined as the composition

\begin{align}
\bF_{h^\#,\Gamma} = \beta_{\mathcal{Z}} \circ h^\# \circ \alpha_{\mathcal{X}}\circ\Gamma\circ \calF_{\mathcal{X}}: L^2(V) \mapsto L^2(V).
\end{align}

Recall the eigenbasis $\set{u_{\bX,i} \given i\in[n]}$ in <ref> for each $\bX\in \mathcal{X}$. Similarly, for each $\bZ\in \mathcal{Z}$, we assume that its eigenvectors form an orthonormal eigenbasis $\set{u_{\bZ,i} \given i\in[n]}$ of $\Real^n$. The filters $\bF_{h^*(\Gamma)}$ and $\bF_{h^\#,\Gamma}$ have the following forms: for $f\in L^2(V)$,

\begin{align}
\bF_{h^*(\Gamma)}(f) &= \int_{\mathcal{Z}} \sum_{i=1}^n \Gamma(h(\bZ),i)\ip{f}{u_{\bZ,i}} u_{\bZ,i} \ud\mu_{\mathcal{Z}}(\bZ),\\
\bF_{h^\#,\Gamma}(f) &= \int_{\mathcal{Z}} \sum_{i=1}^n \Gamma(h(\bZ),i)\ip{f}{u_{h(\bZ),i}} u_{h(\bZ),i} \ud\mu_{\mathcal{Z}}(\bZ).
\end{align}

If we examine the formulas, we get the intuition that $\bF_{h^\#,\Gamma}$ performs a fiberwise convolution and aggregates according to $\mu_{\mathcal{Z}}$. It can be viewed as a “re-arrangement” of “probability densities”. Effectively, it corresponds to the pushforward of measure. On the other hand, recall that $\Gamma\in L^2(\mathcal{X}\times [n])$ gives rise to a family of filters $\mathcal{X} \mapsto M_n(\mathbb{R})$ by $\bX\mapsto \bigstar_{\Gamma_{\bX}}$. Change of filter family essentially replaces “$\bX$” on the subscript by “$h(\bZ)$”. Hence, as we observe that $\bF_{h^*(\Gamma)}$ “re-arranges” the convolution kernel $\Gamma$, it is indeed the pullback of a filter family. Note that $\bF_{h^*(\Gamma)}$ is a convolution filter associated with $\mathcal{Z}$ as in <ref>. However, $\bF_{h^\#,\Gamma}$ may not be a convolution filter, as illustrated in the following example.

* In this example, we consider two cases where either $\mathcal{Z}$ or $\mathcal{X}$ is finite.

* Suppose $\mathcal{Z}$ is a finite subset of $\mathcal{X}$ and $h: \mathcal{Z} \mapsto \mathcal{X}$ is inclusion, i.e., $h(\bZ)=\bZ$ for all $\bZ\in\calZ$. Then the pushforward measure on $\mathcal{X}$ is a discrete measure supported on $\mathcal{Z} \subset \mathcal{X}$. The pullback of any filter, via $h$, is just the restriction to $\mathcal{Z}$. In this case, $\bF_{h^\#,\Gamma} = \bF_{h^*(\Gamma)}$ performs the following: apply a convolution, with the convolution kernel $\Gamma$ restricted to $\bZ$, at each $\bZ\in \mathcal{Z}$; and then take the expectation according to the discrete measure on $\mathcal{Z}$. The resulting filter is equivalent to the convolution filter on $\mathcal{Z}\times[n] \subset \mathcal{X}\times[n]$.

* Suppose $\mathcal{Z}$ is parametrized by $T= [0,1]$ with the Lebesgue measure $\mathfrak{m}$ and $\mathcal{X} = \{\bX_1,\ldots, \bX_k\}$ is parametrized by a finite subset $\{t_1,\ldots, t_k\}$ of $T$. Suppose $T = \bigcup_{1\leq i\leq k} T_i$ has a decomposition into disjoint intervals such that the length $\mathfrak{m}(T_i)>0$ and $t_i \in T_i$ for each $i=1,\ldots,k$. The map $h: \mathcal{Z} \mapsto \mathcal{X}$ is induced by sending the interval $T_i$ to $t_i$, for each $1\leq i\leq k$. The pushforward measure on $\mathcal{X}$ assigns $\mathfrak{m}(T_i)$ to $t_i$, as well as $\bX_i$. If $\Gamma = (\Gamma_{\bX_i} \in \mathbb{R}^n)_{1\leq i\leq k}$ is a function on $\mathcal{X}\times [n]$, then the filter $\bF_{h^\#,\Gamma}$ performs the following: apply a pointwise convolution, with convolution kernel $\Gamma_{\bX_i}$ at each $\bX_i$; and then take the expectation as the weighted sum with weights $\mathfrak{m}(T_i)$ for $1\leq i\leq k$. This is the coarsening procedure. Since $h$ is not injective, $\bF_{h^\#,\Gamma}$ is not a convolution filter and hence $\bF_{h^\#,\Gamma}\neq \bF_{h^*(\Gamma)}$.

*  Recall the setting of <ref><ref>. Suppose $\mathcal{X}, \mathcal{Z}$ are both parametrized by $T_{\mathcal{Z}} = T_{\mathcal{X}}=[0,1]$ equipped with the Lebesgue measure. Let $G$ be a square lattice and $\bL_0, \bL_1$ are the Laplacians of the subgraphs consisting of horizontal and vertical edges respectively. As in <ref><ref> (notice the notations are changed to deal with two spaces), for $x, z\in [0,1]$, we have the matrix $\bL_x = (1-x)\bL_0+x\bL_1$, and similarly $\bL_z = (1-z)\bL_0+z\bL_1$. For $\eta>0$, define $h: T_{\mathcal{Z}} \mapsto T_{\mathcal{X}}$ by the formula $h(z) = z\eta/(1-z+z\eta) \in [0,1]$. The map $h$ is invertible with inverse given by $x \mapsto x/(x+\eta-x\eta)$. It can be verified that if $\bH_x= x\bL_1 + (1-x)(\eta \bL_0)$ as a stretched version of $\bL_x$, then $\bH_x=\eta/(1-z+z\eta) \bL_z $ is a scalar multiple of $\bL_z$ for $x = h(z)$. In particular, $\bH_x$ and $\bL_z$ have the same eigenbasis.

As a consequence, suppose our knowledge on the graph distribution is on the unstretched version $\bL_z$ and the signal is stretched in the horizontal direction by a factor $\eta$. Then to match prior knowledge and the observed signal in a convolution process, one needs to use the filter with base change $\bF_{h^\#,\Gamma}$. An illustration is shown in fig:bc

The figures show the case $\eta=2$, i.e., the graph is stretched horizontally by a factor of $2$. Suppose $f$, e.g. the image below, is stretched and knowledge of the graph distribution is based on unstretched graphs. To correctly perform signal processing, one needs to apply base change.

§ LEARNING THE SHIFT DISTRIBUTION

The framework discussed in this paper relies on knowing a probability space $(\mathcal{X}, \mu_{\mathcal{X}})$ of shift operators. When such information is not directly available, we propose a Bayesian learning framework to estimate $(\mathcal{X},\mu_{\mathcal{X}})$ with sample data, which essentially involves the MCMC method. For this, we follow largely <cit.>, which we briefly discuss here for completeness.

Suppose we have a candidate space $\mathcal{X}$ and want to learn a discrete distribution $\mu_{\mathcal{X}}$ that approximates the true distribution. We need the following data:

* There is a set of training signals: $\scD = \set{f_i \in L^2(V)\given 1\leq i\leq m}$, possibly with labels $\scL = \set{z_i\in\Real \given 1\leq i\leq m}$.

* There is a loss function to be minimized: $\ell: M_n(\mathbb{R})\times L^2(V) \mapsto \mathbb{R}_+$ in the unlabeled case, and $\ell : M_n(\mathbb{R})\times L^2(V)\times \Real \mapsto \Real_+$ in the labeled case.

* There is a prior measure $\mu_0$ on $\mathcal{X}$. We usually choose the uninformative prior, i.e., the uniform measure.

For each $\bX\in \mathcal{X}$, we may now define the empirical risk as
$\theta(\bX) = \frac{1}{m} \sum_{1\leq i\leq m}\ell(\bX,f_i)$
in the unlabeled case, and
$\theta(\bX) = \frac{1}{m} \sum_{1\leq i\leq m}\ell(\bX,f_i,z_i)$
in the labeled case.

Suppose $\mu_0$ has density function $p_0$ (wrt some dominating measure like Lebesgue measure). For a fixed parameter $\gamma > 0$, discrete samples to approximate $\mu_{\mathcal{X}}$ are drawn proportional to $\exp\parens*{-\gamma \theta(\cdot)}p_0(\cdot)$, yielding the Gibbs posterior. The main insight from <cit.> is that if we treat $\ell$ as a prediction loss, then the learned distribution is an approximation of the actual distribution in the following sense: the expected “prediction” with the learned distribution has a good average performance. The exact statements are called the PAC-Bayesian inequalities <cit.>. To generate such samples, one may use the Metropolis-Hastings algorithm <cit.>.

For each application, it is important to choose an appropriate loss function $\ell$. As $\ell$ depends largely on the explicit situation, we shall describe case-by-case choices in <ref>.

We next discuss base changes in the unlabeled case, in parallel with <ref>. The labeled case can be treated similarly. As earlier, we assume that there is a measurable function between measure spaces $h: \mathcal{Z} \mapsto \mathcal{X}$. We can pushforward a measure from $\mathcal{Z}$ to $\mathcal{X}$ or pullback a measure from $\mathcal{X}$ to $\mathcal{Z}$ as in <ref>, depending on whether the training data $\scD$ is associated with $\mathcal{X}$ or $\mathcal{Z}$.

For the loss function $\ell$, the pullback version is natural to define without additional assumptions. More precisely, $h^*(\ell)$ is defined as $h^*(\ell)(\bZ,f) = \ell(h(\bZ),f)$. It induces pullback of the risk $h^*(\theta)$. As earlier, to define the pushforward $h_*$, we require fiberwise measures $\{\mu_{h^{-1}(\bX)}\}_{\bX\in \mathcal{X}}$, and $h_*(\ell)(\bX,f) = \int_{h^{-1}(\bX)}\ell(\cdot,f)\ud\mu_{h^{-1}}(\bX)$. It induces pushforward of the risk $h_*(\theta)$.

Therefore, depending on the situation, the learned $\mu_{\mathcal{X}}$, using the framework of <cit.> as described above, is proportional to one of the following:

* $\exp(-\gamma \theta(\cdot))p_0(\cdot)$,

* $\exp(-\gamma \theta(\cdot))h_*(p_0)(\cdot)$,

* $\exp(-\gamma h_*(\theta)(\cdot))p_0(\cdot)$,

* $\exp(-\gamma h_*(\theta)(\cdot) \big)h_*(p_0)(\cdot)$.

For $\mu_{\mathcal{Z}}$, we just replace $h_*$ by $h^*$ in the above expressions.

We study the MNIST dataset[http://yann.lecun.com/exdb/mnist/] using the setup of <ref><ref>. This means that we use a $2$D-lattice $G=(V,E)$ to model each image. Let $\bL_0$ and $\bL_1$ be the Laplacians of subgraphs of $G$ consisting of horizontal and vertical edges, respectively. $\mathcal{X}$ is parametrized by the unit interval $[0,1]$; and $t\in [0,1]$ gives rise to $\bL_t = (1-t)\bL_0+t\bL_1$. The size of each image is $28\times 28$. The pixel values can be viewed as a signal $f$ on $V$.  We find a distribution on $\mathcal{X}$ based on sparse encoding of $f$ with each $t\in [0,1]$. More precisely, let $\hat{f}(t,\cdot) : [784] \mapsto \mathbb{R}$ be the Fourier transform of $f$ w.r.t. $\bL_t$. The loss function is:

\begin{align*}
\ell(\bL_t,f)^2 = \frac{\sum_{400\leq i\leq 784} |\hat{f}(t,i)|^2}{\norm{f}^2}.
\end{align*}

For each digit $j=0,1,\ldots, 9$, the empirical distribution $\mu_{j,\mathcal{X}}$ of $\mathcal{X}$ is shown in fig:dsp5. It is interesting to observe that for several digits, e.g., digit $1$, there is an obvious shift of the distribution away from the center $t=0.5$.

Histogram of samples for each digit, with $t$ as the $x$-axis and the number of samples as the $y$-axis for each subplot.

For an $f$, we add i.i.d. Gaussian noise to each pixel and the resulting signal is denoted by $f'$. On each of the distributions $\mu_{j,\mathcal{X}}$, $j=0,\ldots,9$, we design a simple convolution filter $\bigstar_{\Gamma_j}$ induced by $\Gamma_j(t,i)=1$ if $i\leq 400, t\in [0,1]$ and $\Gamma_j(t,i)=0.1$ if $i> 400, t\in [0,1]$ (cf. <cit.>). We apply $\bigstar_{\Gamma_j}$ to the noisy image $f'$ and examples are shown in fig:dsp6. Alongside, we also show the image $f'$ and the image signal obtained by applying an ordinary convolution filter $\bigstar_{\Gamma_G}$ (constructed similarly as above) with $\bL_G = 2\bL_{0.5}$. We see that in general $\bigstar_{\Gamma_j}$ produces arguably sharper images of the digits, with the contrast between a digit and its surrounding region higher. More examples are shown in the supplementary material.

For each digit, there are three images, corresponding to the noisy image $f'$, the image processed with $\bigstar_{\Gamma_G}$ and the image processed with $\bigstar_{\Gamma_j}$, respectively from left to right.

§ NUMERICAL RESULTS

In this section, we present simulation results to showcase applications of the framework proposed in this paper. We first present an example where the graph construction can be based on different attributes. We then consider a sampling and recovery problem on a weather station dataset where the underlying graph can be constructed using a $k$-NN approach, and anomaly detection on a ECoG dataset where graphs are constructed by thresholding pairwise node correlations. Finally, we present a network infection spreading example to illustrate the use of base changes.

§.§ IMDB dataset: heterogeneous graphs

We associate movies in the IMDB dataset[<https://www.imdb.com/interfaces/>] with nodes $V$ of a graph. We form two graphs: an actor graph $G_a$ where two nodes are connected by an edge if they share a common actor (or actress), and a director graph $G_d$ where two nodes are connected if they share a common director. The graph $G_a$ is denser as each movie has multiple actors. Movies are categorized into a few classes, resulting in a signal $f$ on $V$ consisting of integer labels. We assume $f$ is corrupted by (integer) noise. More specifically, we add independent additive white Gaussian noise to each entry of $f$. Each entry in the resulting signal is then rounded to the nearest integer. Let the final corrupted signal be $f_c$. In our experiments, we add different amounts of noise to vary the SNR as $-5$ dB, $-3$ dB and $-1$ dB, so that $f_c$ has a significant number of wrong labels.

Let $\bL_{G_a}$ and $\bL_{G_d}$ be the graph Laplacians of $G_a$ and $G_d$, respectively.
We apply convolution filters to $f_c$ to recover $f$ by scaling down high-frequency components, under the following frameworks:

* Classical GSP with the graph shift operator being either $\bL_{G_a}$ or $\bL_{G_d}$. Let $r_1, r_2 \geq 0$ be scaling factors and $c \in [n]$ be a cutoff threshold. These hyperparameters are denoted as $\omega = \{r_1,r_2,c\}$. The convolution filter we use is to apply the mask $g_{\omega}$ in the graph frequency domain, where $g_{\omega}(i) = r_1$ if $i\leq c$ and $g_{\omega}(i) = r_2$ if $i>c$. The hyperparameters $\omega$ are tuned using $30$ samples.

* Classical GSP with the graph shift operator being a mixture of $\bL_{G_a}$ and $\bL_{G_d}$, i.e., we choose an operator from $\set{\bL_t= t\bL_{G_a} + (1-t)\bL_{G_d}\given t\in (0,1)}$. We apply the same convolution filter given by the frequency mask $g_{\omega}$ as above, with the hyperparameters $\omega$ tuned using $30$ samples. In our experiment, we perform an exhaustive search over $19$ uniformly spaced $t$s in $(0,1)$ and show the result for the operator $\bL_{t^*}$ with the best performance.

* Proposed framework with $\calX = \{\bL_{G_a}, \bL_{G_d}\}$ and $\mu_{\calX}$ uniform over $\calX$. We use the convolution filter $\bigstar_{\Gamma}$ where $\Gamma(\bL_G,\cdot) = g_{\omega_G}(\cdot)$, for $G \in \set{G_a, G_d}$. The hyperparameters $\omega_{G_a}$ and $\omega_{G_d}$ are tuned using $30$ samples.

Boxplots of recovery accuracy by applying convolutions in diffeent frameworks. The figures are based on noisy and hence accuracy levels of $f_c$.

For performance evaluation, we compute the label accuracy of the recovered signal. We show in <ref> the boxplots of the results. Each plot is based on $300$ $f_c$'s with the same noise level. We notice that using $\bL_{G_a}$ and $\bL_{t^*}$ have almost identical performance, while using $\bL_{G_d}$ has much lower accuracy due probably to $G_d$ being more sparse. On the other hand, our proposed approach using $\bigstar_{\Gamma}$ outperforms the other approaches. It can more effectively make use of both $G_a$ and $G_d$ than taking a mixture of $\bL_{G_a}$ and $\bL_{G_d}$.

§.§ Weather station network: sampling and recovery

In this experiment, we study a real weather station network in the United States with $n=194$ nodes.[http://www.ncdc.noaa.gov/data-access/] Though geographic distances between pairs of stations are available, there is no explicit graph connecting the stations. The signals are based on daily temperature over 2013. By preliminary inspection, we notice that the signals are smooth (details are provided in the supplementary material). We want to estimate temperature readings over the entire network based on those sampled at $5$ stations.

To adopt the framework of this paper, we parametrize $\mathcal{X}$ by $k=5,10,\ldots, 190$. Using known geographical distance, for each $k$, we associate it with the $k$-NN graph $G_k$ and obtain the Laplacian $\bL_k$ of $G_k$. For each signal $f$, let $\hat{f}(\bL_k,\cdot)$ be the usual GFT of $f$ w.r.t. $\bL_k$. To learn a distribution on $\mathcal{X}$, we define the loss function

\begin{align}\label{sloss}
\ell(\bL_k,f)^2 = \frac{\sum_{6\leq i\leq 194}|\hat{f}(\bL_k,i)|^2}{\norm{f}^2}.
\end{align}

Intuitively, the loss function $\ell$ computes the fraction, in norm, of the high frequency components of $f$ w.r.t. $\bL_k$. For sampling, we want to find $G_k$ that best compresses the signals, as we want to sample at only a few stations. Therefore, it is reasonable to choose this loss function.
To construct the empirical risk with $\ell$ in the Bayesian framework in <ref>, we use $30$ randomly chosen signals, less than $10\%$ of the total number of signals. The resulting learned distribution $\mu_{\mathcal{X}}$ over the parameter space $\set{5,\ldots,190}$ of $\mathcal{X}$ is shown in fig:dsp1. We see that local peaks occur at $k=10$ and $k=40$, with the weights dropping sharply after $k=90$. Based on the observation, we further restrict $\mathcal{X} = \{5,\ldots,40\}$.

Empirical distribution for $k=5,\ldots, 190$.

For the sampling task, we follow <ref>. We choose $\mathcal{Y}$ to be $\mathcal{X}\times [5]$ and obtain $\bB_{\mathcal{Y}}$ using the empirical distribution $\mu_{\mathcal{X}}$. We apply the sampling and recovery procedure described in <ref> by choosing $V_{189}$ (cf. <ref>) consisting of $5$ stations. For each $f$, let $f'$ be the recovered signal as in <ref> using the recovery matrix $\bG_{V_{189}}$. We evaluate the performance by computing mean error $\epsilon$ between $f$ and $f'$ over all stations. On the other hand, for $k=5,\ldots, 90$, we apply the same procedure with the delta distribution $\delta_k$ at $k$ on $\mathcal{X}$. It is nothing but the procedure of recovery of bandlimited signals with $\bL_k$ as in classical GSP.

The stations are sampled randomly. However, $\bG_{V_{189}}$ associated with $V_{189}$ can be close to being singular. We perform $200$ trials with non-singular $\bG_{V_{189}}$. For each trial, we compute the average of $\epsilon$ over the whole year. The same is done for $\delta_k, k=5,\ldots,90$. Boxplots of the results are shown in fig:dsp14. We see that working with $\mu_{\mathcal{X}}$ has the overall best performance as compared with any $\delta_k$.

Boxplots of average error $\epsilon$. The solid magenta curve shows the mean for each case. It has local minimums as compared to the peaks of fig:dsp1. The dashed magenta line shows the mean for $\mu_{\mathcal{X}}$.

For either $\delta_k$ or $\mu_{\mathcal{X}}$ and chosen $V_{189}$, the determinant $\det(\bG_{V_{189}})$ of $\bG_{V_{189}}$ is another indication of sampling quality, as almost singular $\bG_{V_{189}}$ does not permit good recovery (e.g., the default numerical precision of MATLAB is $16$ digits). We randomly sample $V_{189}$ and show, in fig:dsp13, boxplots of $\log\det(\bG_{V_{189}})$ for different distributions used. We observe that with $\mu_{\mathcal{X}}$, the recovery matrix is less likely to be singular.

Boxplots of $\log\det$ of the recovery matrices.

§.§ ECoG dataset: anomaly detection

We consider the ECoG dataset corresponding to two periods (so-called “pre-ictal” and “ictal”) of a seizure in an epilepsy patient.[<https://math.bu.edu/people/kolaczyk/datasets.html>] ECoG signals are measurements taken at each of the $76$ electrodes in the brain of the patient. We test the performance of our framework with the anomaly detection task, by treating “pre-ictal” signals as normal and “ictal” signals as abnormal. We preprocess by normalizing the signals.

For graph construction, we follow <cit.>. Briefly, there are $76$ nodes associated with $76$ electrodes. To construct a graph, one first computes signal correlations between pairs of nodes. For a chosen $\tau<1$, the graph $G_{\tau}$ is then obtained by thresholding pairwise correlations with $\tau$. We form a sample space of graphs as $\mathcal{X} = \set{G_{\tau} \given \tau=0.25, 0.3, \ldots, 0.8}$. Let $\bL_{\tau}$ be the Laplacian of $G_{\tau}$, and its Fourier basis is $\set{u_{\tau,i}\given i\leq 76}$.

For each $\tau$, by preliminary inspection, we notice that normal (pre-ictal) signals tend to have smaller high frequency components. This prompts us to compute for a signal $f$, the norm $e_{f,\tau}$ of a high-pass filter applied to $f$ as: $e_{f,\tau}^2 = \sum_{i=60}^{76}\langle u_{\tau,i},f\rangle^2$. For each $\tau$, we assume that there is a known $\epsilon_\tau$, and $f$ is declared to be abnormal if $e_{f,\tau}>\epsilon_\tau$. Here, $\epsilon_{\tau}$ is obtained by average $e_{f,\tau}$ for a small sample of both pre-ictal and ictal signals. By going through every $\tau$, we notice that the top $3$ parameters are $\tau = 0.35, 0.4, 0.55$ with accuracy $76.4\%, 76.0\%, 74.6\%$ respectively.

To apply our framework, we first estimate an empirical distribution of $\mu_{\mathcal{X}}$ following <ref>. We randomly choose a sample consisting of $\kappa$ fraction of all signals. The label $y_f$ for a signal $f$ is $1$ if $f$ is ictal and $0$ otherwise. We modify the $0$-$1$-loss (cf. <cit.> Section 2) for the loss function $\ell(\bL_{\tau},f, y_f) = |1(e_{f,\tau}>\epsilon_{\tau})-y_f|$. The empirical distribution $\mu_{\mathcal{X}}$ depends on both $\kappa$ and chosen samples.

For anomaly detection, given signal $f$, we aggregate the normalized difference associated with high-pass filter norms $e_{f,\tau}$ and $\epsilon_{\tau}$: $b_f = \mathbb{E}_{\mu_{\mathcal{X}}} [e_{f,\tau}-\epsilon_{\tau}]/\epsilon_{\tau}$. The signal $f$ is declared to be abnormal if $b_f>0$. We show the detection accuracy in fig:dsp11 for different $\kappa$. We see that the distributional approach generally outperforms using any single $\bL_{\tau}$. As $\kappa$ increases, the general trend is that the performance improves and the standard deviation decreases. However, both changes are very gradual, and in practice $\kappa \approx 17.5\%$ seems to be sufficient.

The figure shows boxplots of the anomaly detection accuracy. For comparison, horizontal lines show the performance for using single graph operators: $\bL_{0.35}$ (top, blue), $\bL_{0.4}$ (middle, red) and $\bL_{0.55}$ (bottom, yellow).

We have noticed that for a single operator, $\bL_{0.35}$ performs the best. We want to investigate its role in the distributional approach by computing its probability weight in each $\mu_{\mathcal{X}}$. The results are shown in fig:dsp12. We notice that as $\kappa$ increases, the standard deviation decreases as expected. On the other hand, the median stays approximately constant near $0.3$. This suggests that contributions from operators other than $\bL_{0.35}$ are also significant.

Boxplots of the probability weight for $\bL_{\tau}$, where $\tau=0.35$, in the distribution $\mu_{\mathcal{X}}$.

§.§ Network infection spreading: base change

In this subsection, we consider infection propagation on a network to illustrate base change (cf. <ref>). Suppose $G=(V, E)$ is a graph and an infection (e.g., a piece of information or disease) is being propagated from a source $s\in V$. Recall that any node receiving the infection is called infected, and an infected node remains infected. A propagation path is a spanning tree $\mathcal{T}$ of $G$ <cit.>. For each $v\in V$, there is a unique path from $s$ to $v$ in $\mathcal{T}$. Therefore, $\mathcal{T}$ describes precisely how the infection is passed from $s$ to the entire graph $G$. Given a set of candidate sources $S\subset V$, our probability spaces $\calZ$ or $\calX$ are identified with pairs $(s, \mathcal{T})$, where $s\in S$ is understood as the source and $\mathcal{T}$ is a propagation path. Specifically, $\calZ$ or $\calX$ consists of the adjacency matrices of propagation paths $\calT$, where we allow self-loops at each vertex.

Transmission rates between vertices can be different. We assume that transmission between certain vertices in the graph is fast. For example, two co-workers in the same office may receive a piece of information at almost the same time or may display symptoms of a disease at around the same time. The collection of such connections $F$ is a subset of edges. Suppose we are given samples $\scD$ of pairs $(s,\mathcal{T})$, where any sample in $\scD$ does not contain edges with the fast transmission. For example, the training data based on contact tracing during a disease lockdown may not have included infections between co-workers while the new observations after lifting the lockdown do, and we need to learn an updated distribution for the current observations. The approach illustrated here can be adapted and applied to other types of changes in the graph attributes. To account for this in our inference, we need to perform a base change.

For ease of analysis, assume that $F$ does not contain cycles. Let $\calZ$ be the space of adjacency matrices corresponding to propagation paths that may be without edges from $F$ and $\calX$ be the space containing adjacency matrices that have all edges in $F$. We define a map $h: \calZ \mapsto \mathcal{X}$ as follows (an example is given in fig:ct):

* Given a pair $(s,\mathcal{T})$, let $\{\rho_1,\ldots,\rho_k\}$ be the edges in $F$ sorted according to the edges' distances to $s$ (the smaller of the distances of $s$ to either end points of each edge). At step $i=0$, we let $\mathcal{T}_0=\mathcal{T}$.

* Suppose in the $(i-1)$-th step, we have a spanning tree $\mathcal{T}_{i-1}$. If $\rho_i \in F$, we do not make changes and set $\mathcal{T}_i=\mathcal{T}_{i-1}$. Otherwise, suppose the end points of $\rho_{i}$ are $v_i, w_i$. There is a unique path $P$ in $T_{i-1}$ connecting $v_i$ and $w_i$. We find an edge $\rho$ in $P\backslash F$ whose distance to $s$ is the median among all edges in $P\backslash F$. Let $\mathcal{T}_i = \mathcal{T}_{i-1}\cup \set{\rho_i} \backslash \set{\rho}$.

We give an example of $h: \calZ \mapsto \mathcal{X}$: the red edges belong to $F=\{\rho_1,\rho_2,\rho_3,\rho_4\}$ and dotted edges are not contained in the spanning tree. The blue node $s$ is the source. Starting from $\mathcal{T}_0$, we have $\mathcal{T}_1=\mathcal{T}_2=\mathcal{T}_0$ (not shown in the illustration) as both $\rho_1,\rho_2$ are contained in $\mathcal{T}_0$. To obtain $\mathcal{T}_3$, we want to include $\rho_3$. We first identify the path in $\mathcal{T}_2=\mathcal{T}_0$ that connects the end points of $\rho_3$, which is the big cycle removing $\rho_3$ with $9$ edges. Seven of them do not belong to $F$. Among them, we pick $\rho$ who has median distance to $s$. The tree $\mathcal{T}_3$ is obtained by adding $\rho_3$ and removing $\rho$. Similarly, from $\mathcal{T}_3$, we add $\rho_4$ and remove the edge next to it to obtain $\mathcal{T}_4$, and $h(s,\mathcal{T}_0)=(s,\mathcal{T}_4)$.

Suppose that a snapshot observation of the infection status of the graph vertices at a particular time is given. We want to infer the source $s$. The snapshot observation gives rise to a graph signal $f$ which is $1$ at infected nodes and $0$ otherwise. We may interpret the source identification problem as learning a distribution on $\mathcal{X}$ with a single $f$, and marginalize $\mathcal{T}$ to find $s$ with the largest likelihood. For a $(s,\mathcal{T}) \in \mathcal{X}$, we define the loss $\ell$ as follows: Write $A_{\mathcal{T}}$ for the sum of the identity matrix and the adjacency matrix of $\mathcal{T}$. Let $\tau$ be the function that sends any non-zero components of a vector to $1$, and $\delta_s$ be the signal of unit length that is supported at $s$. If $H_{\mathcal{T}}$ is the height of $\mathcal{T}$ rooted at $s$, then we take $\ell((s,\mathcal{T}), f) = \min_{0\leq i\leq H_{\mathcal{T}}}\norm{\tau \big(A_{\mathcal{T}}^i(\delta_s)\big)-f}$, where $A_{\mathcal{T}}^i$ is applying the shift $A_{\calT}$, $i$ times.

Let $\mu_{\calZ}$ be the empirical probability mass function on the training set $\scD$. The distribution $\mu_{\mathcal{X}}$ on $\mathcal{X}$ sought after is proportional to $\exp(-\gamma \ell(\cdot,f))h_*(\mu_{\calZ})$ (cf. <ref>).

We perform simulations on a $2$D-lattice and Enron email graph with $225$ and $300$ nodes respectively. The source $s$ is uniformly randomly chosen from a list of $20\%$ candidates, and the propagation path $\mathcal{T}$ is uniformly randomly chosen from the breadth-first search (BFS) trees rooted at $s$, which is a reasonable spreading assumption <cit.>. The snapshot observation is made when around $40\%$ of nodes are infected. We run simulations for $F$ with sizes $20\%$, $40\%$, $60\%$, and $80\%$ of $|E|$. If $s^*$ is the estimated source node, we evaluate the performance by measuring its distance $d$ to the actual source $s$. For comparison, we summarize, in <ref>, the percentage improvement in the error distance with the base change using the mapping $h$ described above, over that without any base change.

In general, we do benefit from the base change. This is more prominent when $|F|$ is large, as expected. For small $|F|$, the performance working without base change may have a slightly better performance.

Improvement in error distance with base change.

$|F|/|E|$     $20\%$     $40\%$     $60\%$     $80\%$

% improvement     $-5.6\%$     $15.2\%$     $77.6\%$     $68.4\%$

$|F|/|E|$     $20\%$     $40\%$     $60\%$     $80\%$

% improvement     $9.65\%$     $27.7\%$     $41.2\%$     $19.1\%$

Enron email graph

§ CONCLUSION

We have presented a new GSP framework over a probability space of shift operators. This is useful in applications where the underlying graph topology is uncertain or where we do not know a priori what is a shift operator consistent with the observations. We develop the concepts of Fourier transform, convolution, and band-pass filters. We discuss and develop methods to allow a change of the underlying probability space of shift operators, which we call a base change. Finally, the usefulness of our framework is demonstrated with numerical experiments on both synthetic and real datasets.

For future work, we will explore how to gain knowledge of the distribution with even less prior information available, so that the current framework can be applied. Another interesting future research direction is to explore the possibility of using the framework to develop new graph neural network methods.

§ PROOFS OF THEORETICAL RESULTS

[Proof of <ref>]
From <ref>, it is clear that $\bigstar_{\Gamma}$ is linear. To show that it is bounded, taking norm of <ref>, we have from the triangle and Cauchy-Schwarz inequalities,

\begin{align*}
\norm{\bigstar_{\Gamma}(f)}
&\leq \int_{\mathcal{X}} \sum_{i=1}^n |\Gamma(\bX,i)| \cdot \norm{f} \ud\mu_{\mathcal{X}}(\bX)  \\
&\leq C \norm{\Gamma}_{L^2(\mathcal{X}\times[n])} \cdot \norm{f},
\end{align*}

for some constant $C$. To show the expectation form, we note that for all $f\in L^2(V)$,

\begin{align*}
\bigstar_{\Gamma}(f)(v)
&= \int_{\mathcal{X}} \sum_{i=1}^n \Gamma_{\bX}(i)\ip{f}{u_{\bX,i}} u_{\bX,i}(v) \ud\mu_{\mathcal{X}}(\bX)\\
&= \int_{\mathcal{X}} \Gamma_{\bX} * f (v) \ud\mu_{\mathcal{X}}(\bX) = \E_{\mu_{\mathcal{X}}}[\bigstar_{\Gamma_{\bX}}](f)(v),
\end{align*}

where the last equality holds because $\bigstar_{\Gamma_{\bX}}$ is a linear operator and can hence be written as an $n\times n$ matrix whose entries are functions of $\bX$. The result then follows from an interchange of integral and finite sum.

[Proof of <ref>]
If $\bX\in \mathcal{X}$ does not have repeated eigenvalues, then the space of fiberwise convolutions $\bigstar_{\Gamma_{\bX}}$ is isomorphic to the space of degree $n-1$ polynomials in $\bX$ (cf. <cit.>). Moreover, we have seen that $\bigstar_{\Gamma}$ is the expectation of $\bigstar_{\Gamma_{\bX}}$ from <ref>, and the result follows.

[Proof of <ref>]
We remark that the space of convolution filters on $L^2(V)$ is finite-dimensional, as it is a subspace of the finite-dimensional space $M_n(\mathbb{R})$ (see <ref>). We want to show that for each $\epsilon>0$ and $\bigstar_{\Gamma}$, there is a bi-polynomial filter $\bF$ such that the $\norm{\bF-\bigstar_{\Gamma}}\leq \epsilon$, where $\norm{\cdot}$ is the operator norm. However, all bi-polynomial filters form a subspace of the space of convolution filters. The above approximation property cannot hold for these two finite-dimensional vector spaces unless they are the same.

Let $c>0$ be the upper bound on the operator norm of $\bX_t^i, t\in T$, $0\leq i\leq n-1$ for almost every $\bX\in \mathcal{X}$. For $\bX_t \in \mathcal{X}$ with no repeated eigenvalues, $\bigstar_{\Gamma_{\bX_t}} = \sum_{0\leq i\leq n-1}a_i(t)\bX_t^i$ for some $a_i(t) \in L^2(T)$. The Weierstrass approximation theorem <cit.> says that any continuous function on $T$ can be approximated arbitrarily closely by a polynomial on $T$ with uniform norm. Moreover, the space of continuous functions is dense in $L^2(T)$ <cit.>. As a consequence, we can find a polynomial $b_i(t)$ such that $\norm{b_i(t)-a_i(t)}_{L^2(T)}$ is as small as we wish, say bounded by $\epsilon/nc$. Let $\bigstar_{\Gamma'_{\bX_t}} = \sum_{0\leq i\leq n-1}b_i(t)\bX_t^i$. We have

\begin{align*}
& \norm*{\E_{\mu_T}[\bigstar_{\Gamma_{\bX_t}}]-\E_{\mu_T}[\bigstar_{\Gamma'_{\bX_t}}]}_{L^2(T)} \\
= & \norm*{\E_{\mu_T} \sum_{0\leq i\leq n-1}(a_i(t)-b_i(t))\bX_t^i}_{L^2(T)} \\
\leq & \sum_{0\leq i\leq n-1} \norm{b_i(t)-a_i(t)}_{L^2(T)} \cdot c
\leq n \frac{\epsilon}{nc} c = \epsilon.
\end{align*}

This proves the claim of the first paragraph and hence the theorem.

[Proof of <ref>]
For $(\mathcal{Y},\epsilon)$-bandlimited signals $f_1,f_2$ and $0\leq a\leq 1$, let $f = af_1+(1-a)f_2$. We have $\norm{\bB_{\mathcal{Y}}(f)-f} \leq a\norm{\bB_{\mathcal{Y}}(f_1)-f_1}+(1-a)\norm{\bB_{\mathcal{Y}}(f_2)-f_2} \leq \epsilon$. This shows that $(\mathcal{Y},\epsilon)$-bandlimited signals form a convex set.

Moreover, from <ref> and the fact that norms are continuous, $f\mapsto \norm{\bB_{\mathcal{Y}}(f)-f}$ defines a continuous function $\mathbb{R}^n \mapsto \mathbb{R}$. If $\epsilon>0$, the inverse image of $(-\epsilon,\epsilon)$ is an open subset of $\mathbb{R}^n$ containing $0$. Hence, $0$ is an interior point of the set of $(\mathcal{Y},\epsilon)$-bandlimited signals.

[Proof of <ref>]
From <ref>, we have $\bB_{\mathcal{Y}} = \E_{\mu_{\mathcal{X}}}[\bF_{\bX}]$. In the expression, $\bF_{\bX}$ is a graph band-pass filter (in classical GSP <cit.> Section V.A.) parametrized by $\bX$ and its eigenvalues belong to $[0,1]$.  We want to show that $\bB_{\mathcal{Y}}$ is positive semi-definite and its operator norm is bounded by $1$.

For any signal $f\in L^2(V)$, we have

\begin{align*}
\norm{\bB_{\mathcal{Y}}(f)} = \norm*{\E_{\mu_{\mathcal{X}}}[\bF_{\bX}](f)} \leq \E_{\mu_{\mathcal{X}}}\norm{\bF_{\bX}(f)} \leq \norm{f}.
\end{align*}

On the other hand, we also have the lower bound

\begin{align*}
\ip{\bB_{\mathcal{Y}}(f)}{f} = \int_{\mathcal{X}} \ip{\bF_{\bX}(f)}{f}\ud\mu_{\mathcal{X}}(\bX) \geq \int_{\mathcal{X}} 0 \ud\mu_{\mathcal{X}} = 0.
\end{align*}

[Proof of <ref>]
As $f = \sum_{1\leq i\leq n}a_ie_i$, $\bB_{\mathcal{Y}}(f)-f = \sum_{1\leq i\leq n} a_i(\lambda_i-1)e_i$. Using orthogonality of $e_i$, $1\leq i\leq n$, we have $\norm{\bB_{\mathcal{Y}}(f)-f}^2 = \sum_{1\leq i\leq n}(1-\lambda_i)^2a_i^2$. As $f$ is $(\mathcal{Y},\epsilon)$-bandlimited, $\sum_{1\leq i\leq n}(1-\lambda_i)^2a_i^2\leq \epsilon^2$. Therefore, $\sum_{1\leq i\leq j}a_i^2\leq \epsilon^2/(1-\lambda_j)^2.$

[Proof of <ref>]

* Suppose we express $f$ as a vector in the coordinates given by the basis $u_1,\ldots,u_n$. We decompose $f = f_1+f_2$ where $f_1$ belongs to the span of $u_1,\ldots, u_j$ and $f_2$ belongs to the span of $u_{j+1},\ldots, u_n$. Correspondingly, we express $f_{V_j}= f_{1,V_j}+f_{2,V_j}$ where $f_{i,V_j},i=1,2$ takes the $V_j$-components of $f_i$. Hence, $f' = \bU_{>j}\bG_{V_j}^{-1}(f_{1,V_j})+\bU_{>j}\bG_{V_j}^{-1}(f_{2,V_j})$.

As $\bG_{V_j}$ is the recovery matrix associated with the uniqueness set $V_j$, we have $f_2 = \bU_{>j}\bG_{V_j}^{-1}(f_{2,V_j})$. Therefore, $\norm{f'-f} = \norm{f_1-\bU_{>j}\bG_{V_j}^{-1}(f_{1,V_j})}$.
On the other hand, since $f$ is $(\mathcal{Y},\epsilon)$-bandlimited, $\norm{f_1} \leq \epsilon/(1-\lambda_j)$ from <ref>. Hence, we have the following estimation:

\begin{align*}
& \norm{f_1-\bU_{>j}\bG_{V_j}^{-1}(f_{1,V_j})}  \leq \norm{f_1} +  \norm{\bU_{>j}\bG_{V_j}^{-1}(f_{1,V_j})} \\
= & \norm{f_1} +  \norm{\bG_{V_j}^{-1}(f_{1,V_j})}\leq \norm{f_1} + \sigma_{V_j}\norm{f_{1,V_j}} \\
\leq & \norm{f_1}(1+ \sigma_{V_j})
\leq  \epsilon\frac{1+\sigma_{V_j}}{1-\lambda_j}.
\end{align*}

* Using <ref> and the condition that $f$ is $(\mathcal{Y},\epsilon)$-bandlimited

\begin{align*}
& \norm{\bB_{\mathcal{Y}}(f')-f'} & \\
= &\norm{\bB_{\mathcal{Y}}(f'-f) + (f'-f) + (\bB_{\mathcal{Y}}(f)-f)} \\
\leq &\norm{\bB_{\mathcal{Y}}(f'-f)} + \norm{f'-f} + \norm{\bB_{\mathcal{Y}}(f)-f}\\
\leq & 2\norm{f'-f} + \epsilon
\leq \epsilon \parens*{1+2\frac{1+\sigma_{V_j}}{1-\lambda_j}}.
\end{align*}

The proof is now complete.