Generalized Sampling on Graphs With Subspace and Smoothness Priors

TL;DR

This paper introduces a generalized sampling framework for graph signals that accommodates arbitrary signals and different sampling and reconstruction filters, utilizing priors like subspace and smoothness for improved recovery.

Contribution

It extends classical sampling theory to graph signals with novel correction filter design and two priors, enabling flexible and effective sampling and reconstruction methods.

Findings

Effective recovery with correction filters in graph domain

Comparison shows improved performance over existing methods

Numerical experiments validate the proposed framework

Abstract

We propose a framework for generalized sampling of graph signals that parallels sampling in shift-invariant (SI) subspaces. This framework allows for arbitrary input signals, which are not constrained to be bandlimited. Furthermore, the sampling and reconstruction filters may be different. We present design methods of the correction filter that compensate for these differences and lead to closed form expressions in the graph frequency domain. In this study, we consider two priors on graph signals: The first is a subspace prior, where the signal is assumed to lie in a periodic graph spectrum (PGS) subspace. The PGS subspace is proposed as a counterpart of the SI subspace used in standard sampling theory. The second is a smoothness prior that imposes a smoothness requirement on the graph signal. We suggest the use of recovery techniques for when the recovery filter can be optimized and under a setting in which a predefined filter must be used. Sampling is performed in the graph frequency domain, which is a counterpart of "sampling by modulation" used in SI subspaces. We compare our approach with existing sampling techniques on graph signal processing. The effectiveness of the proposed generalized sampling approach is validated numerically through several experiments.

Tables

Table 1. TABLE I: Correction and Reconstruction Filters for Shift-Invariant and Graph Spectral Filters where CF and RF are abbreviations of correction filter and reconstruction filter, respectively. DS, LS, and MX refer to direct-sum, least squares, and minimax solutions, respectively. Spectra R X ​ Y ​ ( ω ) subscript 𝑅 𝑋 𝑌 𝜔 R_{XY}(\omega) and R ~ X ​ Y ​ ( λ i ) subscript ~ 𝑅 𝑋 𝑌 subscript 𝜆 𝑖 \tilde{R}_{XY}(\lambda_{i}) are defined in ( 8 ) and ( 53 ).

	Shift-invariant subspace				Periodic graph spectrum subspace
	Unconstrained		Predefined ($W\left(\omega \right)$ is fixed)		Unconstrained		Predefined ($W\left({\lambda }_i\right)$ is fixed)
Filter	CF	RF	CF		CF	RF	CF
Subspace Prior	$\frac{1}{R_{SA}\left(\omega \right)}$	$A\left(\omega \right)$	$\frac{R_{WA}\left(\omega \right)}{R_{SA}\left(\omega \right)R_{WW}\left(\omega \right)}$	DS, MX	$\frac{1}{{\tilde{R}}_{SA}\left({\lambda }_i\right)}$	$A\left({\lambda }_i\right)$	$\frac{{\tilde{R}}_{WA}\left({\lambda }_i\right)}{{\tilde{R}}_{SA}\left({\lambda }_i\right){\tilde{R}}_{WW}\left({\lambda }_i\right)}$	DS, MX
			$\frac{1}{R_{SW}\left(\omega \right)}$	LS			$\frac{1}{{\tilde{R}}_{SW}\left({\lambda }_i\right)}$	LS
Smoothness Prior	$\frac{1}{R_{SW}\left(\omega \right)}$	$\frac{S\left(\omega \right)}{{\left|V\left(\omega \right)\right|}^2}$	$\frac{1}{R_{S\tilde{W}}\left(\omega \right)}$	LS	$\frac{1}{{\tilde{R}}_{SW}\left({\lambda }_i\right)}$	$\frac{S\left({\lambda }_i\right)}{V^2\left({\lambda }_i\right)}$	$\frac{1}{{\tilde{R}}_{S\tilde{W}}\left({\lambda }_i\right)}$	LS
			$\frac{R_{W\tilde{W}}\left(\omega \right)}{R_{S\tilde{W}}\left(\omega \right)R_{WW}\left(\omega \right)}$	MX			$\frac{{\tilde{R}}_{W\tilde{W}}\left({\lambda }_i\right)}{{\tilde{R}}_{S\tilde{W}}\left({\lambda }_i\right){\tilde{R}}_{WW}\left({\lambda }_i\right)}$	MX

Table 2. TABLE II: Average MSEs of Reconstructed Signals after 1000 Independent Runs (in Decibels). Columns with BL Refer to Bandlimited Sampling, and Those with Non-BL Refer to Non-Bandlimited Sampling.

		Generator function #1 (70)				Generator function #2 (71)
	Solution/	Noiseless signals		Noisy signals		Noiseless signals		Noisy signals
Prior	Strategy	BL	Non-BL	BL	Non-BL	BL	Non-BL	BL	Non-BL
Subspace	Unconstrained	-296.65	-302.56	-9.66	-10.54	-299.84	-303.80	-10.22	-10.35
	Predefined: DS and MX	-18.52	-18.52	-9.00	-9.92	-12.16	-12.16	-7.91	-8.29
	Predefined: LS*	-13.82	-18.45	-7.84	-9.92	-7.25	-12.10	-5.12	-8.26
Smoothness	Unconstrained	-3.99	-20.30	-3.14	-10.06	-8.80	-18.67	-6.58	-9.87
	Predefined: MX	-5.14	-18.49	-3.95	-9.92	-8.33	-12.15	-6.14	-8.28
BL sampling and reconstruction [40]		-3.99		-3.86		-8.80		-8.41
FastGSSS [10]		-2.31		-1.85		-7.18		-5.88
* Same as the predefined solution for smoothness prior with LS strategy

Equations

$$\hat{x}[i]=⟨{\mathbf{u}}_i,\mathbf{x}⟩=\sum _{n=0}^{N-1}{u}_i^{\ast }[n]x[n].$$

$${\mathbf{x}}_{\text{out}}=\mathbf{G}\mathbf{x}.$$

$$x_{\text{out}}[n]:=\sum _{i\in {\mathcal{N}}_n}[\mathbf{G}{]}_{ni}x[i],$$

$$x_{\text{out}}[n]:=\sum _{i=0}^{N-1}\hat{x}[i]G({\lambda }_i)u_i[n]$$

$${\mathbf{x}}_{\text{out}}=\mathbf{U}G(\mathbf{\Lambda }){\mathbf{U}}^{\ast }\mathbf{x},$$

$$c[n]=⟨s(t-nT),x(t)⟩={x(t)\ast s(-t)|}_{t=nT},$$

$$C(\omega )=R_{SX}(\omega ),$$

$$R_{SX}(\omega ):=\frac{1}{T}\sum _{k=-\infty }^{\infty }{S}^{\ast }\left(\frac{\omega -2\pi k}{T}\right)X\left(\frac{\omega -2\pi k}{T}\right)$$

$$\tilde{x}=WHc=WH(S^{\ast }x),$$

$$\tilde{x}(t)=\sum _{n\in \mathbb{Z}}(h[n]\ast c[n])w(t-nT),$$

$$x(t)=\sum _{n\in \mathbb{Z}}d[n]a(t-nT),$$

$$X(\omega )=D(e^{j\omega T})A(\omega ),$$

$$\mathcal{H}=\mathcal{A}\oplus {\mathcal{S}}^{\perp }.$$

$$\tilde{x}=A(S^{\ast }A{)}^{-1}{S}^{\ast }x=x.$$

$$H(\omega )=\frac{1}{R_{SA}(\omega )}.$$

$$\tilde{x}=\underset{x\in \mathcal{A},S^{\ast }x=c}{\mathrm{argmin}}\Vert x{\Vert }^2,$$

$$\tilde{x}=A(S^{\ast }A{)}^{†}{S}^{\ast }x.$$

$$H(\omega )=\{\begin{array}{ll}\frac{1}{R_{SA}(\omega )}& R_{SA}(\omega )\ne 0\\ 0& R_{SA}(\omega )=0.\end{array}$$

$$\tilde{x}=\underset{\tilde{x}}{\mathrm{argmin}}\underset{x\in \mathcal{A},S^{\ast }x=c}{\max }\Vert \tilde{x}-x{\Vert }^2.$$

$$H=(W^{\ast }W{)}^{-1}{W}^{\ast }A(S^{\ast }A{)}^{-1}.$$

$$H(\omega )=\frac{R_{WA}(\omega )}{R_{SA}(\omega )R_{WW}(\omega )}.$$

$$\tilde{x}=W(S^{\ast }W{)}^{†}{S}^{\ast }x.$$

$$\tilde{x}=W(W^{\ast }W{)}^{-1}{W}^{\ast }A(S^{\ast }A{)}^{†}{S}^{\ast }x,$$

$${\int }_{-\infty }^{\infty }|V(\omega )X(\omega ){|}^{2}d\omega \le {\rho }^2.$$

$$\tilde{x}=\underset{x\in \{x|S^{\ast }x=c\}}{\mathrm{argmin}}\Vert Vx{\Vert }^2.$$

$$\tilde{x}=\tilde{W}(S^{\ast }\tilde{W}{)}^{-1}{S}^{\ast }x$$

$$H(\omega )=\frac{1}{R_{S\tilde{W}}(\omega )}$$

$$\tilde{W}(\omega )=\frac{S(\omega )}{|V^2(\omega )|}.$$

$$\tilde{x}=\underset{x\in \{x|x\in \mathcal{W},S^{\ast }x=Pc\}}{\mathrm{argmin}}\Vert Vx{\Vert }^2,$$

$$\tilde{x}=\hat{W}(S^{\ast }\hat{W}{)}^{†}{S}^{\ast }x$$

Full text

I Introduction

Sampling theory for graph signals has been recently studied with the goal of building parallels of sampling results in standard signal processing [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Since the pioneering Shannon–Nyquist sampling theorem [11, 12], sampling theories that encompass more general signal spaces beyond that of bandlimited signals in shift-invariant (SI) spaces have been widely studied with many promising applications [13, 14, 15, 16, 17, 18, 19]. More relaxed priors have also been considered such as smoothness priors. These theories allow for sampling and recovery of signals in arbitrary subspaces using almost arbitrary sampling and recovery kernels. These results are particularly useful in the SI setting in which sampling and recovery reduce to simple filtering operations.

Graph signal processing (GSP) [20, 21] is a relatively new field of signal processing that studies discrete signals defined on a graph. Recent work on GSP ranges from theory to practical applications including wavelet/filter bank design [22, 23, 24, 25], learning graphs from observed data [26, 27, 28, 29], restoration of graph signals [30, 31], image/point cloud processing [32], and deep learning on graphs [33].

One of the topics of interest in GSP is graph sampling theory [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], which is aimed at recovering a graph signal from its sampled version. Most studies on sampling of graph signals have considered recovery of discrete graph signals from their sampled version [1, 2, 3, 4, 5, 6, 7, 10]. Current approaches generally rely on vertex sampling. The graph can be highly irregular, namely, the number of edges connected to a vertex may vary significantly. Hence, the “best” sampling set depends on the graph and assumed signal model; sampling set selection with different models of signals or features has been studied extensively in sensor networks and machine learning [34, 35, 36, 37]. Graph sampling theory typically assumes that the signal is smooth on the graph, in which the smoothness is frequently measured based on the number of nonzero coefficients in the graph Fourier spectrum [8, 1], allowing for perfect recovery, which will differ depending on the sampling set.

Vertex domain sampling parallels sampling of the (discrete) time domain signals in standard signal processing. However, the sampled spectrum of a graph signal does not preserve its original shape. Whereas time domain sampling has a corresponding frequency (i.e., DFT) domain representation that preserves the shape of the spectrum (up to possible aliasing) [38, 19, 39], vertex domain sampling does not have such a simple relationship. Instead, we propose to build an analog of standard sampling in the graph frequency domain. Graph frequency sampling has been proposed in [40]. Vertex and graph frequency domain sampling approaches coincide under certain conditions (see Section IV-C).

Here, we expand on [40] by building a generalized graph sampling framework that allows for (perfect) recovery of graph signals beyond bandlimited signals, and parallels SI sampling for time domain signals. In SI sampling, the input subspace has a particular SI structure. Sampling is modeled by uniformly sampling the output of the signal convolved with an arbitrary sampling filter. Under a mild condition on the sampling filter, recovery is obtained using a correction filter having an explicit closed-form frequency response. Herein, we demonstrate how one can extend these ideas to graphs by defining an appropriate input space of graph signals and sampling in the graph frequency domain [40]. In addition, our generalized sampling framework enables recovery of non-bandlimited graph signals from vertex domain sampling for bipartite graphs by applying a relationship between vertex and graph frequency domain sampling. In contrast to graph filter bank approaches [41, 23, 42, 43, 24, 44, 45], our proposed generalized framework only requires one branch (i.e., channel) to recover the full-band graph signals.

Our framework relies on graph sampling performed in the graph frequency domain [40] as a counterpart to “sampling by modulation” under the SI setting [19, 39]. This sampling method maintains the shape of the graph spectrum. Unlike in SI sampling, in which sampling in the time domain coincides with that in the frequency domain, under the graph setting, vertex domain sampling and frequency domain sampling are generally different. Sampling by modulation enables a generalized graph sampling framework that is analogous to SI sampling, exhibiting a symmetric structure in which the sampling and reconstruction steps contain similar building blocks as those in SI sampling. Our approach reduces to the standard SI results in the case of a graph representing the conventional time axis whose graph Fourier basis is the discrete Fourier transform (DFT). We consider two priors on graph signals:

Subspace prior, where the signal lies in a known subspace characterized by a given generator; and 2. 2.

Smoothness prior, where the signal is smooth on a given graph.

Both priors parallel those considered in SI sampling [13, 19].

For the subspace prior, we define the periodic graph spectrum (PGS) subspace that serves as a counterpart of SI subspaces. In particular, this subspace maintains the repeated graph frequency spectra of SI signals. In the smoothness prior, we assume that the quadratic form of the graph signal is small for a given smoothness function. Under this setting, perfect recovery is no longer possible. Nonetheless, following the work on general Hilbert space sampling, we show how to design graph filters that allow to best approximate the input signal under several different criteria [46, 47, 48, 19].

Generalized sampling for standard and graph sampling paradigms allows for the use of arbitrary sampling and reconstruction filters that are not necessarily ideal low-pass filters. It also allows for fixed recovery filters that may have advantages in terms of implementation. For all settings, and under all recovery criteria considered, we show that reconstruction is given by the spectral graph filters, the response of which has a closed form solution that depends on the generator function, smoothness, and sampling/reconstruction filters.

In the context of subspace sampling with a PGS prior, our results allow for a perfect recovery of graph signals beyond those that are bandlimited for almost all signal and sampling spaces. In particular, we require such subspaces to satisfy a direct-sum (DS) condition, as found in standard generalized sampling. When the DS condition does not hold, we design a correction filter that best approximates the input under both the least-squares (LS) and minimax (MX) criteria. These criteria have been studied in the context of standard sampling. We then introduce LS and MX strategies for recovery under a smoothness prior. In all cases, the graph filters have explicit graph frequency responses that parallel those in the SI setting.

Earlier work focusing on generalized sampling of graph signals can be found in [4]. This approach is based on generalized Hilbert space sampling [13, 19] and demonstrates the possibility of perfect recovery of graph signals that are not necessarily bandlimited. However, its framework does not parallel SI sampling in general, i.e., the reconstruction matrix does not have a simple diagonal graph frequency response. Likewise, most previous studies have considered vertex domain subsampling, including many graph sampling studies [1, 2, 3, 4, 5, 6, 7, 10], resulting in different building blocks in the sampling and reconstruction steps. Our framework, by contrast, leads to simple closed form recovery methods based on graph filters for both the sampling and recovery. We expand on the similarities and differences between our study and previous approaches in Section VII.

In our preliminary study [49], we considered generalized graph sampling with a subspace prior. In this study, the results are significantly expanded by introducing an integrated framework, applying different design criteria, and further considering a smoothness prior.

The remainder of this paper is organized as follows. The notations and basics of GSP are introduced in Section II. Section III reviews generalized sampling in Hilbert spaces and in the SI setting. A framework for generalized graph sampling is presented in Section IV. Section V proposes signal recovery methods assuming a PGS subspace prior. We describe the relaxation this prior into a smoothness prior in Section VI. Section VII describes the relationships between our work and existing methods. Numerical experiments are presented in Section VIII. Finally, Section IX provides some concluding remarks regarding this research.

II Spectral Graph Theory and Basics of GSP

We begin by reviewing graphs and their spectrum. We also introduce some basic GSP operators, like the graph Fourier transform (GFT) and filtering on graphs.

A graph $\mathcal{G}$ is represented as $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ and $\mathcal{E}$ denote sets of vertices and edges, respectively. The number of vertices is given by $N=|\mathcal{V}|$ unless otherwise specified. We define an adjacency matrix $\mathbf{A}$ with elements $a_{mn}$ that represents the weight of the edge between the $m$th and $n$th vertices; $a_{mn}=0$ for unconnected vertices. The degree matrix $\mathbf{D}$ is a diagonal matrix, with $m$th diagonal element $[\mathbf{D}]_{mm}=\sum_{n}a_{mn}$.

GSP uses different variation operators [20, 21] depending on the application and assumed signal and/or network models. Here, for concreteness, we use the graph Laplacian $\mathbf{L}:=\mathbf{D}-\mathbf{A}$ or its symmetrically normalized version $\underline{\mathbf{L}}:=\mathbf{D}^{-1/2}\mathbf{L}\mathbf{D}^{-1/2}$. The extension to other variation operators (e.g., adjacency matrix) is possible with a slight modification for properly ordering its eigenvalues as long as the graph is undirected without self-loops. Because $\mathbf{L}$ is a real symmetric matrix, it always possesses an eigendecomposition $\mathbf{L}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{*}$, where $\mathbf{U}=[{\bm{u}}_{0},\ldots,{\bm{u}}_{N-1}]$ is a unitary matrix containing the eigenvectors $\bm{u}_{i}$, and $\bm{\Lambda}=\text{diag}(\lambda_{0},\lambda_{1},\ldots,\lambda_{N-1})$ consists of the eigenvalues $\lambda_{i}$. We refer to $\lambda_{i}$ as the graph frequency.

A graph signal $x:\mathcal{V}\rightarrow\mathbb{C}$ is a signal that assigns a value to each vertex. It can be written as a vector $\bm{x}$ in which the $n$th element $x[n]$ represents the signal value at the $n$th vertex. The GFT is defined as

[TABLE]

Our generalized sampling can also use other GFT definitions, e.g., [50, 51, 52, 53], without changing the framework.

A (linear) graph filter is defined as $\mathbf{G}\in\mathbb{C}^{N\times N}$. The filtered signal is represented as

[TABLE]

Graph filtering may be defined in the vertex and frequency domains. Vertex domain filtering is defined as a linear combination of the neighborhood samples

[TABLE]

where $\mathcal{N}_{n}$ represents neighborhood vertex indices around the $n$th vertex. In graph frequency domain filtering, the output is defined as a generalized convolution [54]:

[TABLE]

where the filter response in the graph frequency domain is given by $G(\lambda_{i})\in\mathbb{R}$. This filtering is equivalently written as

[TABLE]

where $G(\bm{\Lambda}):=\text{diag}(G(\lambda_{0}),G(\lambda_{1}),\dots)$. Here, $\mathbf{G}=\mathbf{U}G(\bm{\Lambda})\mathbf{U}^{*}$. If $G(\lambda_{i})$ is a $P$th order polynomial, (4) coincides with vertex domain filtering (3) with a $P$-hop local neighborhood [20].

III Generalized Sampling in Hilbert Space

This section introduces prior results on generalized sampling in Hilbert spaces [14, 15, 19] and corresponding results in the SI setting, which are fundamental for our generalized graph sampling approach. This section briefly describes the generalized sampling framework. Detailed derivations, including error analysis which can be easily applied to graph sampling, may be found in [19]. Table I summarizes the main results of this section in the SI setting.

III-A Sampling and Recovery Framework

Fig. 1(a) illustrates the generalized sampling framework in Hilbert space. Its SI counterpart is shown in Fig. 1(b). Let $x$ be a vector in a Hilbert space $\mathcal{H}$ and $c[n]$ be its $n$th sample given by $c[n]=\langle s_{n},x\rangle$, where $\{s_{n}\}$ is a Riesz basis and $\langle\cdot,\cdot\rangle$ is an inner product. Denoting by $S$ the set transformation corresponding to $\{s_{n}\}$ we can write the samples as $c=S^{*}x$, where $\cdot^{*}$ represents the adjoint. The subspace generated by $\{s_{n}\}$ is denoted by $\mathcal{S}$.

In the SI setting, $s_{n}=s(t-nT)$ for a real function $s(t)$ and a given period $T$. The samples can then be expressed as

[TABLE]

where $\ast$ represents convolution. The continuous-time Fourier transform (CTFT) of the samples $c[n]$, denoted $C(\omega)$, can be written as

[TABLE]

where

[TABLE]

is the sampled cross correlation. Thus, we may view sampling in the Fourier domain as multiplying the input spectrum by the filter’s frequency response and subsequently aliasing the result with uniform intervals that depend on the sampling period. In bandlimited sampling, $s(-t)=\text{sinc}(t/T)$, where $\text{sinc}(t)=\sin(\pi t)/(\pi t)$. However, $s(t)$ can be chosen arbitrarily in the generalized sampling framework.

The recovery of the sampled signal $c$ is represented as

[TABLE]

where $W$ is a set transformation corresponding to a basis $\{w_{n}\}$ for the reconstruction space, which spans a closed subspace $\mathcal{W}$ of $\mathcal{H}$. The transform $H$ is called the correction transformation and operates on the samples $c$ prior to recovery.

In the SI setting, the recovery corresponding to (9) is given by

[TABLE]

where a discrete-time correction filter $h[n]$ is first applied to $c[n]$: The output $d[n]=h[n]\ast c[n]$ is interpolated by $w(t-nT)$, to produce the recovery $\tilde{x}(t)$.

Next, we describe known results on generalized sampling with subspace and smoothness priors.

III-B Subspace Prior

Suppose that $x$ lies in an arbitrary subspace $\mathcal{A}$ of $\mathcal{H}$ and assume that $\mathcal{A}$ is known. Hence, $x$ can be represented as $x=\sum d[n]a_{n}=Ad$, where $\{a_{n}\}$ is an orthonormal basis for $\mathcal{A}$ and $d[n]$ are the expansion coefficients of $x$. In the SI setting, $x(t)$ is expressed as

[TABLE]

for some sequence $d[n]$ where $a(t)$ is a real generator satisfying the Riesz condition. In the Fourier domain, (11) becomes

[TABLE]

where $A(\omega)$ is the CTFT of $a(t)$ and $D(e^{j\omega T})$ is the discrete-time Fourier transform (DTFT) of the sequence $d[n]$, and is $2\pi/T$ periodic.

III-B1 Unconstrained Case

We first consider the case in which the recovery is unconstrained, so that $W$ can be any transformation. In this setting, we may recover a signal in $\mathcal{A}$ by choosing $W=A$ in (9). If $S^{*}A$ is invertible, then perfect recovery of any $x\in\mathcal{A}$ is possible by using $H=(S^{*}A)^{-1}$. Invertibility can be ensured by the DS condition: $\mathcal{A}$ and $\mathcal{S}^{\bot}$ intersect only at the origin and span $\mathcal{H}$ jointly. This requirement is formally written as

[TABLE]

Under the DS condition, a unique recovery is obtained by an oblique projection operator onto $\mathcal{A}$ along $\mathcal{S}^{\bot}$ given by

[TABLE]

In the SI setting, the frequency response of the correction filter is

[TABLE]

If $\mathcal{A}$ and $\mathcal{S}^{\bot}$ intersect, then there is more than one signal in $\mathcal{A}$ that matches the sampled signal $c$. We may then consider several selection criteria to obtain an appropriate signal out of (infinitely) many candidates. Widely accepted strategies are the LS and MX approaches.

The LS recovery is the minimum energy solution obtained as

[TABLE]

and is given by

[TABLE]

Here, $H=(S^{*}A)^{\dagger}$ and $\cdot^{\dagger}$ represents the Moore-Penrose pseudo inverse. Its corresponding form in the SI setting is

[TABLE]

The MX criterion minimizes the worst-case error from the original signal:

[TABLE]

The solution with a subspace prior is the same as that in (17).

III-B2 Predefined Case

When the reconstruction transformation $W$ is predefined, perfect recovery is not possible in general. However, we can still design a correction transformation $H$ such that the solution is close to $x$ in some sense.

With the DS condition in (13), a minimal-error recovery can be obtained by the correction filter

[TABLE]

The recovered signal is $\tilde{x}=W(W^{*}W)^{-1}W^{*}A(S^{*}A)^{-1}S^{*}x$, which is the orthogonal projection of the unconstrained solution onto $\mathcal{W}$. In the SI setting,

[TABLE]

When the DS condition does not hold, the LS and MX strategies can be considered as in the unconstrained case. The LS solution is $H=(S^{*}W)^{\dagger}$, which results in the following reconstruction:

[TABLE]

This solution is the same as that in (17) by replacing $A$ with $W$. The MX solution is given by

[TABLE]

with $H=(W^{*}W)^{-1}W^{*}A(S^{*}A)^{\dagger}$. The corresponding SI solution is the same as that in (21) with $H(\omega)=0$ when the denominator is zero.

III-C Smoothness Prior

The smoothness prior is a less restrictive assumption than the subspace prior because the actual signal subspace $\mathcal{A}$ is not necessarily known. Instead, we assume the signal is smooth, which is formulated as $\|Vx\|\leq\rho$ for some invertible operator $V$: In the SI setting $V=V(\omega)$ is nonzero for all $\omega$. Smoothness is often measured by low energy in high frequency components:

[TABLE]

In general, with a smoothness prior, there are infinitely many solutions. Two approaches to select a solution are the LS and MX methods, which can be applied in both the unconstrained and constrained settings.

III-C1 Unconstrained Case

Suppose that $V^{*}V$ is a bounded operator. In the LS method, the objective function is formulated by choosing the smoothest signal among all the possibilities:

[TABLE]

The solution to (25) is given by

[TABLE]

where $\widetilde{W}=(V^{*}V)^{-1}S$. In the SI setting, the correction filter in (26) reduces to

[TABLE]

with

[TABLE]

The MX solution coincides with (26).

III-C2 Predefined Case

When the recovery space is predefined, the constraint on the feasible set is slightly different from that in (25). The LS objective for the predefined case may be formulated as

[TABLE]

where $P$ is the orthogonal projection onto the range space of $S^{*}W$. The solution can be shown to be given by

[TABLE]

where $\widehat{W}=W(W^{*}V^{*}VW)^{-1}W^{*}S$. In the SI setting, (29) reduces to the use of $H(\omega)=1/R_{SW}(\omega)$ prior to reconstruction with $W(\omega)$ [19, Section 7.2.1]. Therefore, constrained recovery under the LS objective is the same in the subspace and smoothness priors and the smoothness constraint is not included in the solution.

The MX criterion with a smoothness prior can be formulated as

[TABLE]

This solution is given by

[TABLE]

This is the orthogonal projection onto $\mathcal{W}$ of the unconstrained solution in (26): The correction transformation is $H=(W^{*}W)^{-1}W\widetilde{W}(S^{*}\widetilde{W})^{-1}$. In the SI setting, it reduces to

[TABLE]

IV Sampling and Recovery of Graph Signals

IV-A Sampling of Graph Signals

Two methods of sampling over graphs have been proposed in the literature: 1) sampling in the vertex domain [1, 2] and 2) sampling in the graph frequency domain [40].

IV-A1 Sampling in the Vertex Domain

For sampling in the vertex domain, samples on a predetermined vertex set $\mathcal{T}$ are selected. This corresponds to nonuniform subsampling in the time domain. In contrast to the SI setting, vertex domain sampling is conducted nonuniformly because vertex indices do not reflect the structure of the signal. Many approaches have been proposed to select the “best” sampling set from a given graph under different criteria [1, 5, 10, 6].

Let us define $\mathbf{I}_{\mathcal{T}}\in\{0,1\}^{K\times N}$ as a submatrix of the identity matrix $\mathbf{I}_{N}$, whose rows are determined by the sampling set $\mathcal{T}$ that identifies the vertices that remain after sampling, i.e., row indices in $\mathbf{I}_{N}$. Sampling in the vertex domain is defined as follows:

Definition 1** (Sampling of graph signals in the vertex domain [1, 2]).**

Let $\bm{x}\in\mathbb{C}^{N}$ be the original graph signal and $\mathbf{G}\in\mathbb{C}^{N\times N}$ be an arbitrary graph filter in (2). In addition, let $\mathbf{I}_{\mathcal{T}}$ be a submatrix of the identity matrix $\mathbf{I}_{N}$ extracting $K=|\mathcal{T}|$ rows corresponding to the sampling set $\mathcal{T}$. The sampled graph signal $\bm{c}\in\mathbb{C}^{K}$ is given as follows:

[TABLE]

The sampling matrix is therefore given by $\mathbf{S}^{*}=\mathbf{I}_{\mathcal{T}}\mathbf{G}$.

Aggregation sampling [3, 7] is a variant of vertex sampling that uses a specifically designed $\mathbf{G}$. For example, [3] defines

[TABLE]

where $[\bm{\Psi}]_{k\ell}=\lambda_{\ell}^{k}$, and [7] utilizes a random matrix to filter the signal, i.e.,

[TABLE]

where $\mathbf{\Xi}$ is a random matrix and $\circ$ represents an element-wise product. In general, $\mathbf{G}$ in (35) and (36) cannot be decomposed using $\mathbf{U}$ and therefore it does not have a diagonal graph frequency response.

The definitions above based on vertex domain operations result in nonuniform sampling in general. When the signal is bandlimited (in a graph frequency sense), perfect recovery is guaranteed if $\mathcal{T}$ is a uniqueness set [1, 8]. However, sampling and reconstruction are not symmetric in general: Recovery is not given through filtering and upsampling. This is a significant difference from SI sampling. Instead, we use frequency domain sampling to build a parallel of generalized SI sampling applied to the graph setting and enable recovery through frequency domain filtering.

IV-A2 Sampling in the Graph Frequency Domain

To define sampling over a graph, we extend sampling in SI subspaces expressed by (7) to the graph frequency domain [40]. In particular, the graph Fourier transformed input $\hat{\bm{x}}$ is first multiplied by a graph frequency filter $S(\bm{\Lambda})$; the product is subsequently aliased with period $K$. This results in the following definition:

Definition 2** (Sampling of graph signals within the graph frequency domain).**

Let $\hat{\bm{x}}\in\mathbb{C}^{N}$ be the original signal in the graph frequency domain, i.e., $\hat{\bm{x}}=\mathbf{U}^{*}\bm{x}$, and let $S(\bm{\Lambda})$ be an arbitrary sampling filter in the graph frequency domain. For any sampling ratio $M\in\mathbb{Z}$, the sampled graph signal in the graph frequency domain111$M$ is assumed to be a divisor of $N$ for simplicity. is given by $\hat{\bm{c}}\in\mathbb{C}^{K}$, where $K=N/M$, and

[TABLE]

In matrix form, the sampled graph signal can be represented as $\hat{\bm{c}}=\mathbf{D}_{\text{\emph{samp}}}S(\bm{\Lambda})\hat{\bm{x}}$ where $\mathbf{D}_{\text{\emph{samp}}}=\begin{bmatrix}\mathbf{I}_{K}&\mathbf{I}_{K}&\ldots\end{bmatrix}$.

Hereafter, we denote the sampling matrix $\mathbf{S}^{*}$ as follows.

[TABLE]

This graph frequency domain sampling “mixes” different frequency components obtained by the GFT. Different eigenvectors represent different variations on the graph [20, 54, 55]. The GFT coefficients of a graph signal provide a notion of a frequency content similar to the DFT; however, the GFT basis varies according to the underlying graph and variation operator used. The (weighted) sum of the two GFT coefficients in (37) can be seen as a counterpart of sampling in the Fourier domain in classical signal processing.

Suppose that $\mathbf{U}^{*}$ is the DFT matrix: For example, the DFT matrix diagonalizes the graph Laplacian $\mathbf{L}$ of the circular graph [56]. In this case, the GFT domain sampling in (37) coincides with that in the DFT domain [40], i.e., the sampled spectrum of (7) $C[i]=C(2\pi i/N)$ ($i=0,\dots,N-1$) yields the same output as in (37).

IV-B Recovery Framework for Generalized Graph Sampling

Our framework for generalized graph sampling is shown in Fig. 2. It parallels sampling in Hilbert and SI spaces, as illustrated in Fig. 1 [15, 13]. In this paper, we assume that sampling, filtering, and reconstruction are all performed in the graph frequency domain. This results in graph filters that can be interpreted as an analog of SI sampling. As in standard sampling theory, three filters are critical in the recovery problem: sampling, correction, and reconstruction filters.

To sample $\bm{x}$, we transform the input into the GFT domain, resulting in $\hat{\bm{x}}=\mathbf{U}^{*}\bm{x}$. The output is subsequently filtered using the sampling filter $S(\bm{\Lambda})$. The filtered signal is downsampled to yield a sampled signal $\hat{\bm{c}}=\mathbf{S}^{*}\bm{x}=\mathbf{D}_{\text{samp}}S(\bm{\Lambda})\hat{\bm{x}}$. In the reconstruction step, $\hat{\bm{c}}$ is filtered by the correction filter $\mathbf{H}=H(\bm{\Lambda}):=\text{diag}(H(\lambda_{0}),\dots,H(\lambda_{N-1}))$. Finally, $\tilde{\bm{d}}=H(\bm{\Lambda})\hat{\bm{c}}$ is upsampled to the original dimension by $\mathbf{D}_{\text{samp}}^{\top}$, and the reconstruction filter $W(\bm{\Lambda}):=\text{diag}(W(\lambda_{0}),\dots,W(\lambda_{N-1}))$ is applied to the upsampled signal. After performing an inverse GFT, we obtain the recovered signal $\tilde{\bm{x}}$. This can be written as $\tilde{\bm{x}}=\mathbf{U}W(\bm{\Lambda})\mathbf{D}_{\text{samp}}^{\top}\mathbf{H}\hat{\bm{c}}$, where the reconstruction matrix is given by $\mathbf{W}:=\mathbf{U}W(\bm{\Lambda})\mathbf{D}_{\text{samp}}^{\top}$.

The primary objective in this framework is to consider the design method of the correction and reconstruction filters, $\mathbf{H}$ and $\mathbf{W}$, that recover the original signal as accurately as possible with a given prior and constraint. We follow the same strategies as that of generalized sampling in Hilbert spaces introduced in Section III: DS, LS, and MX. The solutions with subspace and smoothness priors are presented in Sections III-B and III-C, respectively.

IV-C Sampling and Recovery on Bipartite Graphs

Before describing the filter design methods in the following sections, we introduce an interesting and special case, namely, sampling on bipartite graphs. In this scenario, vertex and spectral domain sampling coincide.

Studies on sampling for bipartite graphs are not only interesting as an interconnection between two sampling paradigms, but also of practical importance. Bipartite graph sampling has been well studied in the context of perfect reconstruction graph filter banks [41, 23, 42, 43, 24, 44], where the original signal is decomposed through low- and high-pass channels. Maximally decimated coefficients are obtained by applying sampling for each channel where the two disjoint vertex sets correspond to the transformed coefficients in each channel. To realize perfect recovery for an arbitrary graph, one needs to decompose the original (non-bipartite) graph into several bipartite graphs with disjoint edge subsets222An arbitrary $\mathcal{X}$-colorable graph can always be decomposed into $\lceil\log_{2}\mathcal{X}\rceil$ bipartite subgraphs with disjoint edge subsets [41, 57].. Multiband decompositions are then recursively conducted on these bipartite graphs. Such a filter bank has multiple channels to realize recovery of full-band graph signals. Instead, our graph sampling framework only has one branch, as shown in Fig. 2, while still allowing for perfect recovery under appropriate conditions.

Suppose that a graph is bipartite having two equal-sized vertex sets. Formally, let $\mathcal{B}=(\mathcal{V}_{1},\mathcal{V}_{2},\mathcal{E})$ be a bipartite graph that contains two disjoint vertex sets $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$, where $|\mathcal{V}_{1}|=|\mathcal{V}_{2}|=N/2$, i.e., edges only exist between $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$. Without loss of generality, we assume that its first $N/2$ vertices correspond to $\mathcal{V}_{1}$ and its last vertices correspond to $\mathcal{V}_{2}$. We also assume that the GFT matrix is the eigenvector matrix of the symmetric normalized graph Laplacian.

Fig. 3 illustrates the vertex domain representation of our generalized sampling framework of Fig. 2, where the sampled signal is transformed back into the vertex domain. Suppose that the reduced-size graph $\mathcal{G}_{\text{reduced}}$ of size $N/2$ is obtained by reconnecting edges within $\mathcal{V}_{1}$ with Kron reduction [58, 59]. The symmetric normalized graph Laplacian of $\mathcal{G}_{\text{reduced}}$ can be obtained from $\underline{\mathbf{L}}$ of $\mathcal{B}$ through the following:

[TABLE]

where $\underline{\mathbf{L}}_{\mathcal{X}\mathcal{Y}}$ is a submatrix of $\underline{\mathbf{L}}$, whose extracted rows and columns from $\underline{\mathbf{L}}$ are specified by $\mathcal{X}$ and $\mathcal{Y}$, respectively.

Considering the simplest direct sampling, i.e., there is no sampling filter. The relationship between sampling in the vertex domain (34) and the vertex domain representation of sampling in the graph frequency domain (37) is given in the following theorem, taken from [25]:

Theorem 1**.**

Suppose that the GFT matrix $\mathbf{U}_{\mathcal{B}}$ is the eigenvector matrix of the symmetric normalized graph Laplacian of $\mathcal{B}$, and $\mathbf{U}_{\emph{reduced}}$ is the eigenvector matrix of $\underline{\mathbf{L}}_{\emph{reduced}}$ in (39). The following relationship then holds:

[TABLE]

This relationship indicates that vertex domain sampling (RHS in (40)) coincides with graph frequency domain sampling (LHS in (40)) under the special case in which the graph is bipartite and the graph filter is the identity operator, i.e., $\mathbf{G}=\mathbf{I}$.

Theorem 1 can be generalized to allow sampling with a sampling filter $\mathbf{G}$ not necessarily being the identity:

Corollary 1**.**

Suppose the same setup as in Theorem 1 where a sampling filter $\mathbf{G}$ is conducted before subsampling, as in (34). If the sampling filter $\mathbf{G}$ is diagonalizable by $\mathbf{U}_{\mathcal{B}}$, then $\mathbf{U}_{\emph{reduced}}\mathbf{S}^{*}$, i.e., the vertex domain representation of graph frequency domain sampling, is identical to $\mathbf{I}_{\mathcal{V}_{1}}\mathbf{G}$.

Proof.

By assumption, $\mathbf{G}$ can be represented as follows:

[TABLE]

Therefore, the sampling matrix $\mathbf{S}^{*}$ in (38) is equal to

[TABLE]

Using (40),

[TABLE]

completing the proof. ∎

Corollary 1 addresses the sampling phase. Similarly, during reconstruction, the correction-then-reconstruction steps in the graph frequency domain, $\mathbf{W}\mathbf{H}\mathbf{U}^{*}_{\text{reduced}}$, can be jointly represented using vertex domain upsampling:

[TABLE]

where $W^{\prime}(\lambda_{i}):=W(\lambda_{i})H(\lambda_{i\text{ mod }N/2})$.

Interestingly, from (43) and (44), the sampling-then-reconstruction step is represented as follows:

[TABLE]

This is the same as vertex domain sampling and reconstruction because $\mathbf{W}^{\prime}$ and $\mathbf{G}$ are both graph filters with diagonal graph frequency responses, as in (41) and (44), which is illustrated in Fig. 4. The sampling-then-reconstruction (45) is regarded as one branch in a two-channel graph filter bank for a bipartite graph [41, 23, 43]. The filter bank approach requires at least two such branches to guarantee perfect recovery. In contrast, our generalized graph sampling realizes perfect reconstruction with one branch as long as the DS condition holds (presented in the next section). This may lead to an implementation advantage in practical applications.

V Graph Signal Recovery with Subspace Prior

In this section, we assume that the signal lies in a known PGS subspace that depends on the given graph. Subsequently, we present two possible correction filters. One is an unconstrained solution that guarantees perfect recovery of the graph signal with an arbitrary choice of the sampling filter. The other is a predefined solution in which a given filter must be used for reconstruction.

V-A PGS Subspace

We first consider a graph signal subspace that parallels the generation process in SI subspaces shown in (11) and (12). As discussed in the previous section, vertex domain sampling is in general a nonuniform sampling operator in contrast to the uniform SI sampling of (11). Hence, we utilize graph frequency domain sampling that mimics “sampling by modulation” in (7).

In (12), the $2\pi/T$-periodic spectrum of the expansion coefficients $D(e^{j\omega T})$ is multiplied by the (non-periodic) generator $A(\omega)$ to obtain the signal spectrum $X(\omega)$. We reflect this characteristic in the signal subspace for the graph setting.

The spectrum of the graph considered herein is finite and discrete. Suppose that we have a length $K$ spectrum $\hat{d}(\lambda_{i})$ ($i=0,\dots,K-1$, $K\leq N$) as the expansion coefficients. Although its original length is finite, we assume that $\hat{\bm{d}}$ is periodic beyond $i\geq K$, as in (12), i.e.,

[TABLE]

Under this assumption, we can naturally define the signal subspace for graph signals as a counterpart of the SI subspace, as follows:

Definition 3**.**

A PGS subspace of a given graph $\mathcal{G}$ is a space of graph signals that can be expressed as a GFT spectrum filtered by a given generator:

[TABLE]

where $A(\lambda_{i})$ is the graph frequency domain response of the generator and $\hat{d}(\lambda_{i})$ is an expansion coefficient.

This signal can be represented in the following matrix form:

[TABLE]

where $\hat{\bm{d}}:=[\hat{d}(\lambda_{0}),\dots,\hat{d}(\lambda_{K-1})]^{\top}$.

Bandlimited graph signals are a special case of signals in a PGS subspace. Suppose that $A(\lambda_{i})$ is a bandlimiting low-pass filter $G_{\text{BL},K}(\lambda_{i})$, i.e.,

[TABLE]

The graph signal $\bm{x}$ generated by (48) completely maintains $\hat{\bm{d}}$, i.e., is $K$-bandlimited under the GFT basis $\mathbf{U}$. However, the graph signal generated by (47) with an arbitrary $A(\lambda_{i})$ is in general not necessarily bandlimited; thus, our generalized sampling introduced in the following sections allows for reconstruction of non-bandlimited graph signals.

A similar assumption as (48) without periodicity is given in [26], where $\hat{\bm{d}}$ is assumed to be a random variable. Our PGS assumption includes this definition: In particular, if $K=N$ and $\hat{\bm{d}}$ is a random variable following a zero-mean multivariate Gaussian distribution, then the subspace defined by Definition 3 coincides with that used in [26]. However, note that we impose no constraint on the generator function or the expansion coefficients. Thus, Definition 3 treats a broader class of graph signals than the above.

Suppose that $T$ in (12) is a positive integer, i.e., the spectra $D(e^{j\omega T})$ are repeated $T$ times within $\omega\in[0,2\pi]$, and $A(\omega)$ in (12) has support $\omega\in[0,2\pi]$. In this case, a sequence $X[i]=\left.D(e^{j\omega T})A(\omega)\right|_{\omega=2\pi i/N}$ ($i=0,\dots,N-1$) corresponds to the DFT spectrum of length $N$. Therefore, this $X[i]$ can be regarded as a graph signal spectrum in a PGS subspace when $\mathbf{U}^{*}$ is the DFT matrix, e.g., the graph $\mathcal{G}$ is a circular graph.

V-B Unconstrained Case

Our solutions for generalized graph sampling can be defined following the general Hilbert space results of Section III. Owing to the definition of the PGS subspace and sampling in the graph frequency domain, the sampling, correction, and reconstruction filters can all be implemented in the graph frequency domain.

V-B1 Recovery Filters

For the unconstrained solution, we use a reconstruction filter $W(\lambda_{i})=A(\lambda_{i})$ in (47). Suppose that the DS condition (13) is satisfied for the signal and sampling subspaces. Following the expression in (14), the signal recovery is given as follows:

[TABLE]

where the correction filter is

[TABLE]

Its graph frequency response is

[TABLE]

where

[TABLE]

The inverse of $\tilde{R}_{SA}(\lambda_{i})$ is well defined under the DS condition. Note the similarity with (15).

The solution for the LS and MX strategies when $\mathcal{A}$ and $\mathcal{S}$ intersect can be derived from (17) as follows:

[TABLE]

The correction filter in this case is $\mathbf{H}=(\mathbf{S}^{*}\mathbf{A})^{\dagger}$, and it has the same graph frequency response as (52) but with $H(\lambda_{i})=0$ for $\lambda_{i}$ with $\tilde{R}_{SA}(\lambda_{i})=0$.

V-B2 Special Cases

Suppose that both the generator and sampling filters are $A(\bm{\Lambda})=S(\bm{\Lambda})=G_{\text{BL},K}(\bm{\Lambda})$ in (49). Subsequently, $H(\lambda_{i})=1$ and no correction filter is required. This is equivalent to the perfect recovery condition for bandlimited graph signals using graph frequency domain sampling [40].

Another interesting case is the bipartite graph introduced in Section IV-C. For example, suppose that $S(\lambda_{i})$ in (41) is $G_{\text{BL},N/2}(\lambda_{i})$ and the generator is $A(\lambda_{i})=G_{\text{IR}}(\lambda_{i})$ with

[TABLE]

where the correction filter again becomes $H(\lambda_{i})=1$; therefore, $W(\lambda_{i})=A(\lambda_{i})=G_{\text{IR}}(\lambda_{i})$. This implies that a non-bandlimited graph signal can be perfectly reconstructed from bandlimited measurements by applying the same filtering as in the generation process without an explicit correction filter.

In addition, as mentioned in (45), our sampling and recovery can be represented by using vertex domain sampling for bipartite graphs. Existing graph filter banks for bipartite graphs, e.g., [41, 23, 43], may require length-$N$ coefficients (for maximally decimated transforms) for reconstruction on the synthesis side. Instead, our framework only needs one channel which requires a length-$K$ spectrum for recovery, as demonstrated in Section VIII-B. If the generator function $A(\lambda_{i})$ of a given signal is losslessly encoded and sent along with the spectrum, we can reconstruct the original signal. This may be regarded as a one-branch compression of a graph signal.

V-C Predefined Case

Suppose that the reconstruction filter $W(\lambda_{i})$ is predefined. The reconstructed signal $\tilde{\bm{x}}$ will in general be different from $\bm{x}$ in this case. As in the unconstrained setting introduced in the previous subsection, the correction transforms in our framework are given through graph spectral filters.

If $\mathcal{A}$ and $\mathcal{S}$ satisfy the DS condition in (13), the solution in (20) reduces to the following:

[TABLE]

The corresponding graph filter is

[TABLE]

If $W(\lambda_{i})=A(\lambda_{i})$, the response above is identical to that of the unconstrained case shown in (52).

Without the DS condition, we can apply the LS and MX strategies. The LS solution is

[TABLE]

where the correction filter $\mathbf{H}=(\mathbf{S}^{*}\mathbf{W})^{\dagger}$ has spectral response

[TABLE]

The MX solution becomes

[TABLE]

with

[TABLE]

The spectral response of the filter now is the same as that in (57) but $H(\lambda_{i})=0$ if the denominator is zero.

The graph correction filters are summarized in Table I. The table demonstrates nicely the similarities with SI sampling.

VI Graph Signal Recovery with Smoothness Prior

The subspace prior introduced in the previous section enables the input graph signal to be recovered perfectly; however, it requires full knowledge of the given graph and generator. In this section, we consider a less restrictive assumption. We still assume that the GFT basis $\mathbf{U}$ is given; however, the generator function $A(\lambda_{i})$ is unknown.

We assume that the graph signal is smooth on the given graph where smoothness is measured by the signal energy in the high graph-frequency components as in the SI setting (24). Although several possible operators exist for measuring signal smoothness on a graph [20], we consider a simple quadratic form of $\bm{x}$:

[TABLE]

where $\mathbf{V}:=\mathbf{U}V(\bm{\Lambda})\mathbf{U}^{*}$ is an arbitrary graph filter with spectral response $V(\lambda_{i})$. The smoothness condition is given by $\|\mathbf{V}\bm{x}\|^{2}_{2}\leq\rho^{2}$ for some constant $\rho$. This can be seen as a generalization of the bandlimitedness of graph signals, which has been widely studied [1, 2, 5, 6, 8, 9], because a bandlimited graph signal corresponds to $\rho=0$ for a high-pass filter $\mathbf{V}=\mathbf{U}(\mathbf{I}-G_{K,\text{BL}}(\bm{\Lambda}))\mathbf{U}^{*}$. In addition, if we assume $\mathbf{V}=\mathbf{L}^{1/2}$, then $\|\mathbf{V}\bm{x}\|^{2}_{2}=\bm{x}^{*}\mathbf{L}\bm{x}$, which is a Laplacian quadratic form also used extensively in the literature. Hereinafter, for simplicity, we assume $V(\lambda_{i})\neq 0$ for all $i$.

For the unconstrained case, the LS recovery is given from (26) as follows:

[TABLE]

where $\widetilde{\mathbf{W}}=(\mathbf{V}^{*}\mathbf{V})^{-1}\mathbf{S}=\mathbf{U}V^{2}(\bm{\Lambda})\mathbf{U}^{*}$ and $\mathbf{S}^{*}\mathbf{V}^{-2}\mathbf{S}=\mathbf{D}_{\text{samp}}S^{2}(\bm{\Lambda})V^{-2}(\bm{\Lambda})\mathbf{D}_{\text{samp}}^{\top}$ is invertible if $\tilde{R}_{SS}(\lambda_{i})\neq 0$ for all $i$. This results in $\mathbf{H}=(\mathbf{S}^{*}\widetilde{\mathbf{W}})^{-1}$, where the spectral response is

[TABLE]

The MX solution coincides with (64) as in the SI solution.

We next consider the predefined case. For the LS criterion, the solution in Hilbert space (29) reduces to the constrained LS solution with a subspace prior (59): This does not depend on $V(\lambda_{i})$, i.e., the smoothness prior does not affect the solution.

The MX solution can be obtained from (32):

[TABLE]

This leads to

[TABLE]

wtih spectral response

[TABLE]

The smoothness prior $V(\lambda_{i})$ is incorporated appropriately in the correction filter, in contrast to the LS solution.

These correction filters are summarized in Table I.

VII Comparison with Existing Studies

VII-A Computational Complexity

Here, we compare the computational complexities of vertex and graph frequency domain sampling required for sampling, recovery, and preprocessing.

VII-A1 Sampling

For spectral domain sampling, its complexity is $\mathcal{O}(NK)$ because the sampling matrix $\mathbf{S}^{*}$ has a size of $K\times N$. However, its complexity is reduced to $\mathcal{O}(N)$ if we already have GFT coefficients $\hat{\bm{x}}$ because the sampling filter response $S(\bm{\Lambda})$ is diagonal. Typically, the GFT is necessary to perform only once even when the sampling ratio or filter is changed.

For vertex domain sampling, the subsampling itself in (34) only picks up elements specified by the sampling set $\mathcal{T}$: Its computation cost is negligible. In contrast, a sampling filter in the vertex domain requires $\mathcal{O}(N^{2})$ complexity in general. This can be reduced using a localized filter: $P$-hop filtering requires $\mathcal{O}(P|\mathcal{E}|)$ complexity [22].

VII-A2 Recovery

As mentioned in the above two sections, all correction and reconstruction filters in our generalized sampling framework have diagonal graph frequency responses, the computational complexity of which is $\mathcal{O}(K)$. When we change the sampling filter, the response of the correction filter can be immediately calculated because it is diagonal. An additional complexity is required for the inverse GFT (where its complexity depends on the GFT used) if reconstructed vertex domain signals are required.

For signal recovery using vertex domain sampling, a matrix of size $N\times K$ is multiplied by the sampled coefficients, the complexity of which is $\mathcal{O}(NK)$. Note that the reconstruction matrix depends on $\mathcal{T}$: If the sampling set $\mathcal{T}$ or sampling rate $|\mathcal{T}|$ is changed, we have to re-calculate the entire reconstruction matrix even if the graph is the same. This calculation typically requires the inversion of a $K\times K$ matrix, e.g., [2, 1, 10], and may lead to $\mathcal{O}(K^{3})$ complexity.

VII-A3 Preprocessing

Graph frequency domain sampling considered in this paper requires the GFT matrix, and hence we need to compute an eigendecomposition of the variation operator. This typically requires $\mathcal{O}(N^{3})$ complexity (whereas we can use several fast computation methods of GFT or spectral decomposition of the graph variation operator such as in [52, 53, 60, 61, 62]). It is important to note that the same graph is often used numerous times. In such a case, we only need to calculate the GFT basis once, and can reuse it even when we change the sampling rate $M=N/K$ or the sampling filter.

Vertex domain sampling always requires computing the best $\mathcal{T}$ from a given graph. The computation complexity highly depends on the sampling set selection methods. Major methods are compared in [10]. Typically, the complexity depends on $N$, the (assumed) cutoff graph frequency, and the edge density.

VII-B Literature Review

In [4], a generalized sampling method for graph signal processing was studied. As the results did not assume any particular structure on the input signals and sampling filters, the recovery procedures were in general given by matrix inversions. Here, we focus on a special case of [4] that extends SI sampling to the graph setting and enables explicit expressions for the recovery filter in the graph Fourier domain.

Our solution represented in (54) allows for a broad choice of $S(\lambda_{i})$ and $A(\lambda_{i})$. In particular, $A(\lambda_{i})$ is not restricted to a bandlimiting operator. If we have $S(\bm{\Lambda})=A(\bm{\Lambda})=G_{\text{BL},K}(\bm{\Lambda})$, our solution reduces to that of [4], which is equivalent to the sampling theory with graph frequency domain sampling studied in [40].

For the smoothness prior, if the smoothness is measured by $\mathbf{V}^{2}=\mathbf{L}+G_{\text{BL},K}(\bm{\Lambda})$, i.e.,

[TABLE]

and $S(\bm{\Lambda})=G_{\text{BL},K}(\bm{\Lambda})$, our solution also reduces to that introduced in [4]: The reconstruction and correction filters as shown in (63) exhibit the following form:

[TABLE]

for $\lambda_{i}\in[\lambda_{0},\lambda_{K-1}]$, and $W(\lambda_{i})=H(\lambda_{i})=0$ otherwise. This is a special case of [4] where the correction and reconstruction operators can be represented as spectral filters.

As mentioned in Section V-A, many studies on graph sampling theory such as [1, 2, 3] implicitly assume that the graph signal lies in the PGS subspace with a typical generator function $A(\lambda_{i})=G_{\text{BL},K}(\lambda_{i})$. While their subspace is a special case of the PGS assumption, the sampling matrices are different from that in (38). As described in Definition 1, the simple subsampling $\mathbf{S}^{*}=\mathbf{I}_{\mathcal{T}}$ has been used in many studies on graph sampling [1, 2]. In [3], aggregation sampling was used. However, its sampling matrix in (35) does not in general have a corresponding sampling expression in the graph frequency domain as that in (37). This results in the requirement of matrix inversion even for recovering the bandlimited graph signal although the signal lies in a PGS subspace.

In summary, most studies on graph sampling theory require inversion of the sampling operator for their reconstruction framework. Moreover, they focused on the design problem for the nonuniform sampling matrix $\mathbf{I}_{\mathcal{T}}$ that, for example, maximizes the bandwidth with perfect recovery. In contrast, frequency domain sampling is utilized in this study as a counterpart of “sampling by modulation” in SI spaces, thus resulting in a symmetric structure, i.e., both the sampling and reconstruction steps can be represented as similar sampling and filtering operations. We also allow for a broader set of input signals and design criteria.

VIII Graph Signal Recovery Experiments

In this section, we validate the proposed generalized sampling through signal recovery experiments. First, we demonstrate that the correction and recovery filters described in Sections V and VI reduce the reconstruction error of non-bandlimited graph signals compared to the bandlimited sampling in the graph frequency domain [40] and FastGSSS [10], which is a state-of-the-art vertex domain sampling method. This reveals the MSE improvements of our generalized sampling for non-bandlimited graph signals over bandlimited or smoothness-based reconstruction. Sampling under the bipartite case presented in Section IV-C is then conducted in which full-band graph signals are almost perfectly recovered with one branch of sampling and reconstruction, even without calculating the GFT matrix.

VIII-A Recovery Experiments for Bandlimited and Non-bandlimited Sampling

We first conduct signal recovery experiments following our generalized sampling framework shown in Fig. 2. Because we have two choices of priors, i.e., subspace and smoothness priors; three strategies, i.e., DS, LS, and MX; and two possible reconstruction filters, i.e., unconstrained or predefined, we compare these settings throughly in this simulation. In addition, two sampling filters are considered: bandlimited and non-bandlimited. The bandlimited sampling corresponds to a graph sampling theory described in [40], whereas our generalized sampling recovers the original full-band graph signals after bandlimited sampling. As mentioned previously, we allow for non-bandlimited sampling filters to obtain the sampled signal $\hat{\bm{c}}$ while still guaranteeing perfect recovery under certain conditions. Here, we also demonstrate this property in this subsection.

The graph used is a random sensor graph with $N=256$. We downsampled the input signal by a factor of two such that $K=32$. We used the following functions:

• •

Generator function.

[TABLE]

• •

Sampling functions.

[TABLE]

• •

Reconstruction function (used only for the predefined solutions).

[TABLE]

• •

Smoothness function (used only for the smoothness prior).

[TABLE]

We set $\epsilon=0.1$. All functions are visualized in Fig. 5. It is worth noting that both $A(\lambda_{i})$ are not bandlimited; therefore, the original signal retains its full band. The generator function in (71) is smoother than (70): Here, (71) decays more rapidly than (70) when $\lambda$ increases. Each element in the expansion coefficients $\hat{\bm{d}}$ is a random variable drawn from $\mathcal{N}(1,1)$. Examples of $\bm{x}$ generated by $A(\lambda_{i})$ in (70) are shown in Fig. 6(a). For the proposed sampling, we perform two samplings to highlight the difference between the proposed sampling and the bandlimited sampling, as in (72).

For comparison, we applied bandlimited signal recovery with GFT domain sampling [40]: $S(\lambda_{i})=W(\lambda_{i})=G_{\text{BL},K}(\lambda_{i})$ with no correction filter $H(\lambda_{i})=1$. In addition, signal recovery using FastGSSS [10] is also conducted. FastGSSS assumes that $\bm{x}$ is smooth on a graph. It determines the sampling set $\mathcal{T}$ and applies a recovery based on a polynomial graph filter with a given kernel. FastGSSS perfectly recovers the original signal when the polynomial order approaches infinity if the signal is bandlimited.

We applied $1,000$ independent runs and calculated the average MSEs. We then repeated the experiments using zero-mean Gaussian noise with variance $\sigma^{2}=0.1$ added to $\bm{x}$.

Table II summarizes the average MSEs for various sampling approaches. The visualization of the reconstructed signals are shown in Fig. 6. Their corresponding spectra are also shown in Fig. 7.

**Noiseless Signals: ** The unconstrained solution for the subspace prior perfectly recovers the original signal with machine precision. The predefined solutions for both the subspace and smoothness priors contain some reconstruction errors; however, in most cases, they are much smaller than those in the bandlimited reconstruction, especially for the generator function #1 in (70). Because the signals generated by (71) are smoother than those of (70), the recovery with the smoothness prior is comparable to the bandlimited reconstruction if bandlimited sampling is applied. By contrast, using the non-bandlimited sampling function, all predefined solutions are superior to the bandlimited sampling/reconstruction. The unconstrained solution with a smoothness prior (63) yields the same results as those in bandlimited sampling and reconstruction when using $G_{\text{BL},K}(\lambda_{i})$ as the sampling filter.

FastGSSS uses vertex domain sampling and assumes smoothness (or bandlimitedness) of the signals. Its MSE is slightly worse than the predefined solutions and comparable to the bandlimited sampling/reconstruction using spectral domain sampling. This is because the signals used in the experiment are not fully bandlimited. As shown in Fig. 6(c), its reconstructed signal is over-smoothed compared to the original signal in Fig. 6(a). In fact, the sampling set selection strategy of FastGSSS (and any other reconstruction method utilizing vertex domain-based sampling) is based on the assumption that the signal is sufficiently smooth, i.e., a special case of the smoothness prior described in Section VI. This leads to existing recovery methods based on graph sampling theory interpolating missing values with a smooth graph filter, which thus have difficulty using the recovery of non-bandlimited signals even when $M$ is increased. As mentioned in Section VII-A, recovery based on vertex domain sampling generally requires a matrix inversion whose computational complexity is typically $O(K^{3})$, and the entire reconstruction matrix needs to be re-calculated even when the sampling set $\mathcal{T}$ is slightly changed. In contrast, our reconstruction matrix only differs in a diagonal matrix $\mathbf{H}$ even when the sampling rate or sampling filter is changed. Therefore, $\mathbf{H}$ can be easily re-calculated, as long as sampling is applied to the same graph.

**Noisy Signals: ** All methods contain increased errors for noisy cases, as expected. The unconstrained solution for the subspace prior demonstrated a significantly worse performance than that of a noiseless case because it did not assume any smoothness of the reconstructed graph signals. The predefined filters, both with subspace and smoothness priors, have demonstrated a performance close to those of unconstrained solutions because their reconstruction filters yield a smooth signal.

Bandlimited reconstruction occasionally outperforms generalized sampling because it removes the noise in the high-graph-frequency band. Therefore, the bandlimited reconstruction for the smoother signals in (71) is comparable to generalized sampling with non-bandlimited sampling. However, for wider-band signals, as in (70), generalized sampling is much better than a bandlimited method even for noisy cases.

FastGSSS is inferior to spectral domain sampling-based approaches in most cases as in the experiment for noiseless signals. Because such a smoothing filter is naturally robust to noise, MSEs of FastGSSS are stabilized for noisy cases (but are still inferior to GFT domain sampling).

VIII-B Recovery Experiment on Bipartite Graphs with Vertex Domain Sampling

We next demonstrate the recovery of non-bandlimited graph signals from vertex domain sampling, as described in Sections IV-C and V-B. To the best of our knowledge, this example is the first attempt to recover full-band graph signals from vertex domain operations without utilizing a multi-band decomposition.

In the signal recovery of this experiment, we set $S(\bm{\Lambda})=G_{\text{BL},N/2}(\bm{\Lambda})$ and $A(\lambda_{i})=W(\lambda_{i})=G_{\text{IR}}(\lambda_{i})$, as described in Section V-B2. That is, the original signal is full-band, whereas the sampled signal is low-pass filtered, i.e., the high-graph-frequency components are discarded after sampling. We assume that the generator information is losslessly available for recovery. The sampling and recovery of this situation is formulated in (45). If graph filters $\mathbf{G}$ and $\mathbf{W}^{\prime}$ can be represented as vertex domain filters, all operations will be performed in the vertex domain, i.e., without applying GFT (or an eigendecomposition of the variation operator). The filters, i.e., $S(\bm{\Lambda})$ and $W^{\prime}(\bm{\Lambda})$ in (45), cannot be represented as vertex domain filters in general. Fortunately, they can be approximated as vertex domain operators by utilizing polynomial approximations of the spectral filter responses. In this experiment, we use a Chebyshev polynomial approximation (CPA) [22, 63]. As described in Section II, the $P$th order a CPA of an arbitrary graph spectral filter corresponds to a vertex domain filter with $P$-hop localization.

The original signal $\bm{x}$ is obtained as follows333Note that the generation process in (75) uses the non-polynomial $\mathbf{W}^{\prime}$.:

[TABLE]

where each element in $\bm{d}$ is a random variable drawn from a normal distribution $\mathcal{N}(1,1)$. Figure 8 shows the average MSEs of the reconstructed signals after $100$ independent runs according to the polynomial order. For comparison, we also plot the MSE of bandlimited reconstruction in which we use $\widetilde{S}(\lambda_{i})=\widetilde{G}_{\text{BL},N/2}(\lambda_{i})$ as the reconstruction filter, of which $\widetilde{\cdot}$ denotes the polynomial approximated filter. As shown in Fig. 8, the reconstruction error decreases monotonically as $P$ increases. The reconstructed signals are also shown in Fig. 9. The bandlimited reconstruction yields large errors, whereas the proposed reconstruction exhibits extremely similar signal values as the original values.

Future studies on this type of generalized sampling, particularly for the non-bipartite case, is an interesting topic for future research.

IX Conclusion

We proposed a framework for generalized sampling of graph signals. We assumed that graph signals lie in a PGS subspace, which extends the SI subspace in standard signal processing to the graph setting. Sampling is defined in the graph frequency domain. We considered two priors for the graph signals, subspace and smoothness priors, which are parallel to those studied for signals in SI subspaces. All filters used in our framework can be represented as graph spectral filters. Numerical experiments demonstrated that our proposed sampling can recover a class of sampled signals broader than that obtained through existing graph sampling theories. We also presented perfect recovery of non-bandlimited graph signals on bipartite graphs without explicit operations in the GFT domain.