Greedy Sampling of Graph Signals

TL;DR

This paper develops universal performance bounds for noisy graph signal sampling and demonstrates that greedy algorithms can effectively approximate optimal sampling sets despite the problem's combinatorial complexity.

Contribution

It introduces the concept of approximate submodularity to provide near-optimal guarantees for greedy sampling in graph signal processing.

Findings

Universal bounds for Bayesian estimation from noisy samples.

Greedy sampling achieves near-optimal performance with explicit bounds.

Application to reducing complexity in kernel PCA.

Abstract

Sampling is a fundamental topic in graph signal processing, having found applications in estimation, clustering, and video compression. In contrast to traditional signal processing, the irregularity of the signal domain makes selecting a sampling set non-trivial and hard to analyze. Indeed, though conditions for graph signal interpolation from noiseless samples exist, they do not lead to a unique sampling set. The presence of noise makes choosing among these sampling sets a hard combinatorial problem. Although greedy sampling schemes are commonly used in practice, they have no performance guarantee. This work takes a twofold approach to address this issue. First, universal performance bounds are derived for the Bayesian estimation of graph signals from noisy samples. In contrast to currently available bounds, they are not restricted to specific sampling schemes and hold for any sampling sets. Second, this paper provides near-optimal guarantees for greedy sampling by introducing the concept of approximate submodularity and updating the classical greedy bound. It then provides explicit bounds on the approximate supermodularity of the interpolation mean-square error showing that it can be optimized with worst-case guarantees using greedy search even though it is not supermodular. Simulations illustrate the derived bound for different graph models and show an application of graph signal sampling to reduce the complexity of kernel principal component analysis.

Figures

Equations

$$\bm{x}=[\begin{array}[]{ccc}\sigma(u_{1})&\cdots&\sigma(u_{n})\end{array}]^{T}\text{,}\quad u_{i}\in\mathcal{V}\text{.}$$

$$\bar{\mathbf{x}}={\mathbf{V}}^H\mathbf{x}\text{.}$$

$$\mathbf{x}={\mathbf{V}}_{\mathcal{K}}{\bar{\mathbf{x}}}_{\mathcal{K}}\text{.}$$

$$\mathbf{z}=\mathbf{H}\mathbf{x}\text{,}$$

$$\mathbf{y}=\mathbf{x}+\mathbf{w}\text{,}$$

$$\hat{\mathbf{z}}(\mathcal{S})=\mathbf{L}(\mathcal{S}){\mathbf{y}}_{\mathcal{S}}\text{,}$$

$$\underset{\mathbf{L}}{\mathrm{minimize}}$$

$$\hat{\mathbf{z}}(\mathcal{S})=\mathbf{L}{\mathbf{y}}_{\mathcal{S}}$$

$${\mathbf{L}}^{\star }\mathbf{C}\left(\mathbf{\Sigma }+{\mathbf{\Lambda }}_w\right){\mathbf{C}}^T=\mathbf{H}\mathbf{\Sigma }{\mathbf{C}}^T\text{.}$$

$${\mathbf{K}}^{\star }(\mathcal{S})=\mathbf{H}{\mathbf{V}}_{\mathcal{K}}{\left({\mathbf{\Lambda }}^{-1}+\sum _{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\mathbf{v}}_i{\mathbf{v}}_i^H\right)}^{-1}{\mathbf{V}}_{\mathcal{K}}^H{\mathbf{H}}^H\text{,}$$

$$\mathbf{K}(\mathbf{L}{\mathbf{y}}_{\mathcal{S}})=E\left[(\mathbf{H}\mathbf{x}-\mathbf{L}\mathbf{C}\mathbf{y})(\mathbf{H}\mathbf{x}-\mathbf{L}\mathbf{C}\mathbf{y}{)}^H|\mathbf{x},\mathbf{w}\right]\text{.}$$

$$\mathbf{K}(\mathbf{L}{\mathbf{y}}_{\mathcal{S}})$$

$$+\mathbf{L}\mathbf{C}(\mathbf{\Sigma }+{\mathbf{\Lambda }}_w){\mathbf{C}}^T{\mathbf{L}}^H\text{.}$$

$$J(\mathbf{L})={\mathbf{b}}^H\mathbf{K}(\mathbf{L}{\mathbf{y}}_{\mathcal{S}})\mathbf{b}$$

$$\frac{\partial J(\mathbf{L})}{\partial {\mathbf{b}}^H\mathbf{L}}=0\iff \mathbf{C}\left(\mathbf{\Sigma }+{\mathbf{\Lambda }}_w\right){\mathbf{C}}^T{\mathbf{L}}^H\mathbf{b}=\mathbf{C}\mathbf{\Sigma }{\mathbf{H}}^H\mathbf{b}\text{.}$$

$$\begin{gathered} {\mathbf{K}}^{\star }({\mathbf{L}}^{\star }{\mathbf{y}}_{\mathcal{S}})=\mathbf{H}{\mathbf{V}}_{\mathcal{K}}\{\mathbf{\Lambda }-\mathbf{\Lambda }{\mathbf{V}}_{\mathcal{K}}^H{\mathbf{C}}^T\times \\ {\left[\mathbf{C}\left({\mathbf{V}}_{\mathcal{K}}\mathbf{\Lambda }{\mathbf{V}}_{\mathcal{K}}^H+{\mathbf{\Lambda }}_w\right){\mathbf{C}}^T\right]}^{-1}\mathbf{C}{\mathbf{V}}_{\mathcal{K}}\mathbf{\Lambda }\}{\mathbf{V}}_{\mathcal{K}}^H{\mathbf{H}}^H\text{.} \end{gathered}$$

$${\mathbf{K}}^{\star }(\mathcal{S})=\mathbf{H}{\mathbf{V}}_{\mathcal{K}}{\left({\mathbf{\Lambda }}^{-1}+{\mathbf{V}}_{\mathcal{K}}^H{\mathbf{C}}^T\mathbf{C}{\mathbf{\Lambda }}_w^{-1}{\mathbf{C}}^T\mathbf{C}{\mathbf{V}}_{\mathcal{K}}\right)}^{-1}{\mathbf{V}}_{\mathcal{K}}^H{\mathbf{H}}^H\text{.}$$

$${\mathbf{V}}_{\mathcal{K}}^H{\mathbf{C}}^T\mathbf{C}{\mathbf{\Lambda }}_w^{-1}{\mathbf{C}}^T\mathbf{C}{\mathbf{V}}_{\mathcal{K}}=\sum _{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\mathbf{v}}_i{\mathbf{v}}_i^H\text{,}$$

$$\underset{\mathcal{S}\subseteq \mathcal{V}}{\mathrm{minimize}}$$

$$\left|\mathcal{S}\right|\le k$$

$${\mathbf{K}}^{\star }(\mathcal{S})=\mathbf{H}{\mathbf{V}}_{\mathcal{K}}\bar{\mathbf{K}}(\mathcal{S}){\mathbf{V}}_{\mathcal{K}}^H{\mathbf{H}}^H\text{,}$$

$$\mathbf{R}(\mathcal{S})=\sum _{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\mathbf{v}}_i{\mathbf{v}}_i^H\text{.}$$

$$\frac{{\left|\mathcal{K}\right|}^2}{\mathrm{Tr}\left[(\mathbf{W}\mathbf{\Lambda }{)}^{-1}\right]+{\bar{\ell }}_{|\mathcal{S}|}}\le \text{MSE}(\mathcal{S})\le \mathrm{Tr}\left(\mathbf{W}\mathbf{\Lambda }\right)\text{,}$$

$$\text{MSE}(\mathcal{S})=\mathrm{Tr}\left[\mathbf{W}{\left({\mathbf{\Lambda }}^{-1}+\sum _{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\mathbf{v}}_i{\mathbf{v}}_i^H\right)}^{-1}\right]\text{,}$$

$$\mathrm{Tr}\left(\mathbf{X}\right)\ge \frac{N^2}{\mathrm{Tr}\left({\mathbf{X}}^{-1}\right)}\text{,}$$

$$\text{MSE}(\mathcal{S})\ge \frac{{\left|\mathcal{K}\right|}^2}{\mathrm{Tr}\left[(\mathbf{W}\mathbf{\Lambda }{)}^{-1}\right]+\mathrm{Tr}\left[{\mathbf{W}}^{-1}\left({\sum }_{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\mathbf{v}}_i{\mathbf{v}}_i^H\right)\right]}\text{,}$$

$$\text{MSE}(\mathcal{S})\ge \frac{{\left|\mathcal{K}\right|}^2}{\mathrm{Tr}\left[(\mathbf{W}\mathbf{\Lambda }{)}^{-1}\right]+{\sum }_{i\in \mathcal{S}}{\lambda }_{w,i}^{-1}{\Vert {\mathbf{v}}_i\Vert }_{{\mathbf{W}}^{-1}}^2}\text{,}$$

$$\text{MSE}(\mathcal{S})\ge \frac{\left|\mathcal{K}\right|}{{\lambda }_{\text{min}}(\mathbf{W}\mathbf{\Lambda }{)}^{-1}+{\ell }_{\text{max}}}\text{,}$$

$${\bar{\ell }}_{|\mathcal{S}|}\ge \frac{{\left|\mathcal{K}\right|}^2-\eta \mathrm{Tr}\left[(\mathbf{W}\mathbf{\Lambda }{)}^{-1}\right]}{\eta }\text{.}$$

$$\left|\mathcal{S}\right|\ge \frac{{\left|\mathcal{K}\right|}^2-\eta \mathrm{Tr}\left[(\mathbf{W}\mathbf{\Lambda }{)}^{-1}\right]}{\eta {\ell }_{\text{max}}}\text{.}$$

$$f\left(\mathcal{A}\cup \{u\}\right)-f\left(\mathcal{A}\right)\le \alpha \left[f\left(\mathcal{B}\cup \{u\}\right)-f\left(\mathcal{B}\right)\right]\text{,}$$

$$\alpha =\underset{\begin{array}{c}\mathcal{A}\subseteq \mathcal{B}\subseteq \mathcal{V}\\ u\notin \mathcal{B}\end{array}}{\min }\frac{f(\mathcal{A}\cup \{u\})-f(\mathcal{A})}{f(\mathcal{B}\cup \{u\})-f(\mathcal{B})}\text{.}$$

Full text

I Introduction

Graph signal processing (GSP) is an emerging field that studies signals supported on irregular domains [1, 2]. It extends traditional signal processing techniques to more intricate data structures, finding applications in sensor networks, image processing, and clustering, to name a few [3, 4, 5]. Extensions of sampling, in particular, have attracted considerable interest from the GSP community [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. This is not surprising given the fundamental role of sampling in signal processing. Sampling methods in GSP are broadly divided into two categories: selection sampling, in which the graph signal is observed at a subset of nodes [10], and aggregation sampling, in which the signal is observed at a single node for many applications of the graph shift [14]. This work focuses on the former.

As in classical signal processing, samples are only useful inasmuch as they represent the original signal. Conditions under which it is possible to reconstruct a graph signal from noiseless samples can be found in [6, 7, 8, 9, 10, 11, 12]. These, however, do not necessarily lead to unique sampling sets. In fact, for finite graphs with bounded weights, it can be shown that almost every sampling set larger than the bandwidth of the signal guarantees perfect reconstruction [9, 10, 11, 12]. In the presence of noise, however, it is not straightforward which of these exponentially many sets performs best, an issue that becomes more severe as the measurement signal-to-noise ratio (SNR) decreases. In general, selecting an optimal sampling set is NP-hard [16, 17, 18, 19].

In [10, 13], this issue was addressed using randomized sampling schemes, for which optimal sampling distributions and performance bounds were derived for different types of graphs and graph signals. For high SNR, it was shown that sampling proportionally to the leverage score (or its square-root) approximates the sampling distribution that minimizes the reconstruction mean-square error (MSE). Alternatively, a convex relaxation approach was adopted in [20], where the sampling set selection problem was cast as a binary semi-definite program (SDP) and solved by relaxing the binary constraint and thresholding the solution. Rounding and truncation can also be used to approximate the solution of this binary problem. Nevertheless, greedy sampling remains pervasive and has proven successful in many applications [10, 11, 9, 21, 12, 20], though performance analyses are available only for surrogate figures of merit of the MSE, such as the log-determinant [15].

To be sure, this success is warranted by the attractive features of greedy algorithms for large-scale problems. First, their complexity is polynomial in the deterministic case and randomized versions exist that are linear in the size of the ground set, which in this case is the number of nodes in the graph [22]. Also, since they build the solution sequentially, they can be interrupted at any time if, for instance, a desired performance level is reached. More importantly, there is an upper bound on the suboptimality of the greedy solution to monotonic supermodular function minimization problems. This is indeed why greedy algorithms are often used in sensor selection, experimental design, and machine learning [19, 16, 18, 17, 23]. However, the main performance measure in GSP, namely the MSE, is not supermodular in general [24].

In this work, we study the reconstruction (interpolation) performance of greedy sampling schemes in GSP and set out to reconcile the empirical success of greedy MSE minimization with the fact that it is not supermodular. First, in contrast to [10, 13], we adopt a Bayesian approach to graph signal estimation and consider sampling to be deterministic (Section II). Then, we derive bounds on the interpolation MSE that are universal in the sense that they hold for all sampling sets and any sampling method (Section III). These universal bounds are explicit, tractable, and provide practical means of benchmarking the MSE performance of any sampling scheme. Numerical analyses show that the bounds are tight when signal and noise are homeoscedastic. Finally, we develop the concept of approximate supermodularity introduced in [24] and provide near-optimal guarantees for the greedy minimization of the interpolation MSE (Section IV). This result justifies the use of greedy sampling set selection in GSP and explains its success.

To illustrate the practical value of these results, we recall that the concept of sampling is also at the core of statistical methods, such as data subsetting and variable selection, that are crucial for big data applications [25, 26]. Kernel methods, in particular, are prone to complexity issues in large data sets. For instance, performing kernel principal component analysis (kPCA) on a data set of size $n$ requires $n^{2}$ kernel evaluations (KEs) and $\Theta(n^{3})$ operations, while extracting projections for new data takes $n$ KEs and $\Theta(np)$ operations, where $p$ is the number of principal components (PCs) retained [27, 28]. We show that this problem can be cast in the context of GSP and that greedy sampling can be used to reduce the complexity of projections by over $95\%$ at a small performance cost (Section V-B).

Notation: Lowercase boldface letters represent vectors ($\bm{x}$), uppercase boldface letters are matrices ($\bm{X}$), and calligraphic letters denote sets ($\mathcal{A}$). We write ${\left|\mathcal{A}\right|}$ for the cardinality of $\mathcal{A}$ and denote the empty set by $\{\}$. Set subscripts refer either to the vector obtained by keeping only the elements with indices in the set ($\bm{x}_{\mathcal{A}}$) or to the submatrix whose columns have indices in the set ($\bm{X}_{\mathcal{A}}$). To say $\bm{X}$ is a positive semi-definite (PSD) matrix we write $\bm{X}\succeq 0$, so that for $\bm{X},\bm{Y}\in\mathbb{C}^{n\times n}$, $\bm{X}\preceq\bm{Y}\Leftrightarrow\bm{b}^{H}\bm{X}\bm{b}\leq\bm{b}^{H}\bm{Y}\bm{b}$, for all $\bm{b}\in\mathbb{C}^{n}$. The set of PSD matrices is denoted $\mathbb{S}_{+}$ and the set of non-negative real numbers is denote $\mathbb{R}_{+}$. Finally, we take the derivative of a function $f$ with respect to an $n\times 1$ vector $\bm{x}$ to yield the $1\times n$ gradient vector, i.e., $\partial f/\partial\bm{x}=\left[\ \partial f/\partial x_{1}\ \cdots\ \partial f/\partial x_{n}\ \right]$ [29].

II Sampling and Interpolation of Graph Signals

A graph-supported signal, or graph signal for short, is an assignment of values to the nodes of a graph. Formally, let $\mathbb{G}$ be a weighted graph with node set $\mathcal{V}$, having cardinality ${\left|\mathcal{V}\right|}=n$, and define a graph signal to be an injective mapping $\sigma:\mathcal{V}\to\mathbb{C}$. For an ordering of the nodes in $\mathcal{V}$, this signal can be represented as an $n\times 1$ vector that captures its values at each node:

[TABLE]

In what follows, we assume that the node ordering is fixed, so that we can index $\bm{x}$ using elements of $\mathcal{V}$. For instance, we write $\bm{x}_{\{u_{i},u_{j},u_{k}\}}=\left[\ \sigma(u_{i})\ \sigma(u_{j})\ \sigma(u_{k})\ \right]^{T}$.

Of interest to GSP is the spectral representation of the signal $\sigma$ (or $\bm{x}$), which depends on the graph on which it is supported. Indeed, let $\bm{A}\in\mathbb{C}^{n\times n}$ be a matrix representation of $\mathbb{G}$. Usual choices include the adjacency matrix or one of the discrete Laplacians [1, 2]. Assume that $\bm{A}$ is consistent with the signal vector (1) in the sense that they employ the same ordering of the nodes in $\mathcal{V}$. Furthermore, assume that $\bm{A}$ is normal, i.e., that there exist $\bm{V}\in\mathbb{C}^{n\times n}$ unitary and $\bm{D}\in\mathbb{R}^{n\times n}$ diagonal such that $\bm{A}=\bm{V}\bm{D}\bm{V}^{H}$, where $\cdot^{H}$ is the Hermitian (conjugate transpose) operator [30]. Then, the graph Fourier transform of $\bm{x}$ is given by [1, 2]

[TABLE]

Observe that if $\bm{A}$ is normal we obtain a spectral energy conservation property analog to Parseval’s theorem in classical signal processing: it is ready to see that $\left\|\bar{\bm{x}}\right\|_{2}=\left\|\bm{x}\right\|_{2}$ if and only if $\bm{V}$ in (2) is unitary, which holds if and only if $\bm{A}$ is normal [30].

Similar to traditional signal processing, a graph signal $\bm{x}$ is said to be spectrally sparse (ssparse) when its spectral representation is sparse. Explicitly, $\bm{x}$ is $\mathcal{K}$-ssparse if $\bar{\bm{x}}$ in (2) is such that $\bar{\bm{x}}_{\mathcal{V}\setminus\mathcal{K}}$ is a zero vector. Then,

[TABLE]

Note that spectrally sparse signals are a superset of bandlimited (“low-pass”) signals. Hence, all results in this work apply to bandlimited signals regardless of the graph frequency order adopted [13, 11, 10].

The interest in $\mathcal{K}$-ssparse or bandlimited graph signals is motivated similarly to traditional signal processing: these signals can be sampled and interpolated without loss of information. Indeed, take sampling to be the operation of observing the value of a graph signal on $\mathcal{S}\subseteq\mathcal{V}$, the sampling set. Then, there exists a set $\mathcal{S}$ of size ${\left|\mathcal{K}\right|}$ such that $\bm{x}$ can be recovered exactly from $\bm{x}_{\mathcal{S}}$ [9, 10, 11, 12]. If, however, only a corrupted version of $\bm{x}_{\mathcal{S}}$ is available, then $\bm{x}$ can only be approximated. To do so, the next section poses noisy interpolation as a Bayesian estimation problem, from which the minimum MSE interpolation operator can be derived. This then allows us to provide universal bounds on the reconstruction error and give near-optimal guarantees for greedy sampling strategies.

II-A Graph signal interpolation

We study graph signal interpolation as a Bayesian estimation problem. Formally, let $\bm{x}\in\mathbb{C}^{n}$ be a graph signal and $\mathcal{S}\subseteq\mathcal{V}$ be a sampling set. We wish to estimate

[TABLE]

for some matrix $\bm{H}\in\mathbb{C}^{m\times n}$ based on the samples $\bm{y}_{\mathcal{S}}$ taken from

[TABLE]

where $\bm{w}\in\mathbb{C}^{n}$ is a circular zero-mean noise vector. By circular we mean that its relation matrix vanishes, i.e., that $\operatorname{\mathbb{E}}\bm{w}\bm{w}^{T}=\bm{0}$ [31]. Note that (4) accounts for scenarios in which we are not interested in the graph signal itself but on a post-processed value, such as the output of a linear classifier or estimator (e.g., Section V-B). The usual graph signal interpolation problem from [24, 32, 9, 10, 11, 12, 13] is recovered by taking $\bm{H}=\bm{I}$.

The prior distribution of $\bm{x}$ reflects the fact that the graph signal is $\mathcal{K}$-ssparse by assuming it is a circular zero-mean distribution with covariance matrix $\bm{\Sigma}=\bm{x}\bm{x}^{H}=\bm{V}_{\mathcal{K}}\bm{\Lambda}\bm{V}_{\mathcal{K}}^{H}$ for $\bm{\Lambda}=\operatorname{diag}(\lambda_{i})$, $\lambda_{i}\in\mathbb{R}_{+}$. We assume without loss of generality that $\bm{\Lambda}$ is full-rank. Otherwise, remove from $\mathcal{K}$ any element $i$ for which $\lambda_{i}=0$. Note that this is equivalent to placing a zero-mean uncorrelated prior on $\bar{\bm{x}}$ in (2). Hence, this model can also be interpreted as the generative model for a wide-sense stationary random process on $\mathbb{G}$ [33, 34, 35]. The noise prior is taken as a zero-mean circular distribution with covariance matrix $\bm{\Lambda}_{w}=\operatorname{diag}(\lambda_{w,i})$, $\lambda_{w,i}\in\mathbb{R}_{+}$ and $\lambda_{w,i}>0$.

We consider estimates of $\bm{z}$ of the form

[TABLE]

for some $\bm{L}(\mathcal{S})\in\mathbb{C}^{n\times{\left|\mathcal{S}\right|}}$. Because $\bm{L}$ recovers (approximates) $\bm{z}$ from the samples $\bm{y}_{\mathcal{S}}$, it is referred to as a linear interpolation operator [10, 13, 11]. An optimal interpolation operator can be found for each $\mathcal{S}$ by minimizing the interpolation error covariance matrix as in

[TABLE]

where $\bm{K}\left[\hat{\bm{z}}(\mathcal{S})\right]=\operatorname{\mathbb{E}}\left[\left(\bm{z}-\hat{\bm{z}}(\mathcal{S})\right)\left(\bm{z}-\hat{\bm{z}}(\mathcal{S})\right)^{H}\mid\bm{x},\bm{w}\right]$ and the minimum is taken with respect to the partial ordering of the PSD cone (see Remark 1). We have omitted the dependence of $\bm{L}$ on $\mathcal{S}$ for clarity. Our interest in solving (7) instead of minimizing the MSE directly is that it is more general. In particular, a solution of (7) also minimizes any spectral function of $\bm{K}$, including the MSE and the log-determinant. The following proposition gives an explicit solution to this problem. To clarify the derivations, we define the selection matrix $\bm{C}\in\{0,1\}^{{\left|\mathcal{S}\right|}\times N}$ composed of the identity matrix rows with indices in $\mathcal{S}$, so that the samples of (5) can be written as $\bm{y}_{\mathcal{S}}=\bm{C}\bm{y}$.

Proposition 1**.**

Let $\bm{y}=\bm{x}+\bm{w}$ be noisy observations of a graph signal $\bm{x}$. Let the priors on $\bm{x}$ and $\bm{w}$ be zero-mean circular distributions with covariances $\bm{\Sigma}=\bm{V}_{\mathcal{K}}\bm{\Lambda}\bm{V}_{\mathcal{K}}^{H}$, $\bm{\Lambda}=\operatorname{diag}(\lambda_{i})$, and $\bm{\Lambda}_{w}=\operatorname{diag}(\lambda_{w,i})$ respectively. Given a sampling set $\mathcal{S}$, the optimal Bayesian linear interpolator $\bm{L}^{\star}$ that solves problem (7) is obtained as a solution of

[TABLE]

The error covariance matrix of the optimal interpolation $\hat{\bm{z}}^{\star}=\bm{L}^{\star}\bm{y}_{\mathcal{S}}$ is given by

[TABLE]

where \bm{V}_{\mathcal{K}}=[\begin{array}[]{ccc}\bm{v}_{1}&\cdots&\bm{v}_{N}\end{array}]^{H}.

Proof.

Start by substituting (4) and (6) into the definition of $\bm{K}$ to get

[TABLE]

Note that we used the fact that $\bm{y}_{\mathcal{S}}=\bm{C}\bm{y}$. Then, using the priors on $\bm{x}$ and $\bm{w}$, $\bm{K}$ expands to

[TABLE]

From the partial ordering of the PSD cone, $\bm{L}^{\star}$ can be obtained by minimizing the scalar cost function

[TABLE]

simultaneously for all $\bm{b}\in\mathbb{C}^{n}$ [29]. Substituting (10) into (11) and setting its gradient with respect to $\bm{b}^{H}\bm{L}$ to zero gives

[TABLE]

Since this must hold for all $\bm{b}$ simultaneously, we obtain (8).

To determine the error covariance matrix $\bm{K}^{\star}$ of the optimal interpolator, replace any $\bm{L}^{\star}$ satisfying (8) into (10) and expand $\bm{\Sigma}=\bm{V}_{\mathcal{K}}\bm{\Lambda}\bm{V}_{\mathcal{K}}^{H}$ to get

[TABLE]

Note that (12) does not depend on $\bm{L}^{\star}$ or $\bm{y}_{\mathcal{S}}$, only on the sampling set $\mathcal{S}$ through the selection matrix $\bm{C}$. Moreover, since $\bm{\Lambda}_{w}$ is diagonal and full rank, $(\bm{C}\bm{\Lambda}_{w}\bm{C}^{T})^{-1}=\bm{C}\bm{\Lambda}_{w}^{-1}\bm{C}^{T}$, so that the inverse in (12) always exists. Therefore, using the matrix inversion lemma [30] gives

[TABLE]

Given that $\bm{C}^{T}\bm{C}$ is a diagonal matrix with ones on the indices in $\mathcal{S}$ and zeros everywhere else, we obtain (9) by noting that

[TABLE]

for \bm{V}_{\mathcal{K}}=[\begin{array}[]{ccc}\bm{v}_{1}&\cdots&\bm{v}_{N}\end{array}]^{H}. ∎

Given prior distributions for the graph signal and noise, Proposition 1 determines the optimal linear interpolator from the samples in $\mathcal{S}$. If the priors on $\bm{x}$ and $\bm{w}$ are moreover Gaussian, then $\hat{\bm{z}}^{\star}=\bm{L}^{\star}\bm{y}_{\mathcal{S}}$ is also the maximum likelihood estimate of $\bm{z}$ [29]. An important consequence of the Bayesian statement of Proposition 1 is that $\bm{\Lambda}$ and $\bm{\Lambda}_{w}$ are taken from prior distributions on the signal and noise. Thus, their actual values need not be known exactly, as illustrated in Section V-B. Note that the optimal error covariance matrix $\bm{K}^{\star}$ now depends only on the sampling set $\mathcal{S}$, since it measures the error of the optimal estimator $\bm{L}^{\star}$. Moreover, although we assume that the interpolation is performed as a single step projection, iterative procedures can also be used [36, 37].

Despite our assumption that $\bm{\Lambda}_{w}$ is full-rank, (8) also holds in the noiseless case ($\bm{\Lambda}_{w}=\bm{0}$). Its solution, however, may no longer be unique. In particular, this happens if the sampling set is not sufficient to determine $\bm{z}$, i.e., if $\bm{C}\bm{V}_{\mathcal{K}}$ is rank-deficient [9, 10, 11, 12]. In contrast, when $\bm{\Lambda}_{w}\succ\bm{0}$, the matrix on the left-hand side of (8) is always invertible and $\bm{L}^{\star}$ is unique for each $\mathcal{S}$. This is similar to the well-known regularization effect of noise in Kalman filtering [29]. The interpolation performance given in (9), however, is not the same for all sampling sets.

Remark 1**.**

Problem (7) is a PSD matrix minimization problem that searches for the optimal interpolator $\bm{L}^{\star}$ that minimizes the error covariance matrix $\bm{K}$. In general, optimization problems in the PSD cone need not have a solution. Since the ordering of PSD matrices is only partial, the existence of a matrix that is smaller than all other matrices is not guaranteed [38]. As shown in Proposition 1, this is not the case here. Problem (7) admits a dominant solution $\bm{L}^{\star}$ in the PSD cone, i.e., it holds that $\bm{K}(\bm{L}^{\star}\bm{y}_{\mathcal{S}})\preceq\bm{K}(\bm{L}\bm{y}_{\mathcal{S}})$ for all $\bm{L}\in\mathbb{C}^{n\times{\left|\mathcal{S}\right|}}$. This means that $\bm{L}^{\star}$ minimizes all the eigenvalues of $\bm{K}$ simultaneously. Equivalently, it implies that $\bm{L}^{\star}$ is a solution to the minimization of any spectral function of $\bm{K}$. In particular, it follows that $\bm{L}^{\star}$ minimizes the MSE, since $\textup{MSE}(\hat{\bm{z}}):=\operatorname{\mathbb{E}}\left\|\bm{z}-\hat{\bm{z}}\right\|_{2}^{2}=\operatorname{Tr}\left[\bm{K}(\hat{\bm{z}})\right]$, and the $\log\det\left[\bm{K}(\hat{\bm{z}})\right]$.

II-B Sampling set selection

Proposition 1 allows us to evaluate the optimal interpolator $\bm{L}^{\star}$ that minimizes the estimation error covariance matrix for a given sampling set. This does not guarantee, however, that there is no other sampling set of the same size for which the interpolation error is smaller. To address this issue, we investigate the sampling set selection problem which sets out to find the sampling set that minimizes the interpolation MSE over all sampling sets. Explicitly, we wish to solve

[TABLE]

where $\textup{MSE}(\mathcal{S})=\operatorname{Tr}[\bm{K}^{\star}(\mathcal{S})]$.

An important fact about (13) is that increasing $\mathcal{S}$ always decreases MSE. This has two important consequences. First, the unconstrained version of (13) is trivial, i.e., its solution is $\mathcal{S}=\mathcal{V}$. Second, it implies that the constraint in (13) is tight, i.e., it can be replaced by the equality constraint ${\left|\mathcal{S}\right|}=k$ without changing the problem solution. This property is a direct corollary of the following lemma and the monotonicity of the trace operator [39]:

Lemma 1**.**

The matrix-valued set function $\bm{K}^{\star}(\mathcal{S})$ in (9) is monotonically decreasing with respect to the PSD cone, i.e., $\bm{K}^{\star}(\mathcal{A})\succeq\bm{K}^{\star}(\mathcal{B})$ whenever $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{V}$.

Proof.

Start by noting that $\bm{K}^{\star}$ in (9) can be written as

[TABLE]

with $\bm{{\bar{K}}}(\mathcal{S})=\left[\bm{\Lambda}^{-1}+\bm{R}(\mathcal{S})\right]^{-1}$ and

[TABLE]

Since $\bm{K}^{\star}$ and $\bm{{\bar{K}}}$ are congruent, it suffices to show that $\bm{{\bar{K}}}$ is a monotonically decreasing set function [30].

To do so, note that $\bm{{\bar{K}}}$ only depends on $\mathcal{S}$ through $\bm{R}(\mathcal{S})$ and that $\bm{R}$ is additive, i.e., $\bm{R}(\mathcal{A}\cup\mathcal{B})=\bm{R}(\mathcal{A})+\bm{R}(\mathcal{B})$. Then, since $\lambda_{w,i}>0$, $\bm{R}$ is a sum of PSD matrices, which implies that $\mathcal{A}\subseteq\mathcal{B}\Rightarrow\bm{R}(\mathcal{A})\preceq\bm{R}(\mathcal{B})$, i.e., $\bm{R}$ is monotonically increasing. From the antitonicity of the matrix inverse [39], it follows that $\bm{{\bar{K}}}$ is monotonically decreasing. ∎

Although Lemma 1 reduces the searching space to sampling sets of size $k$, (13) remains a combinatorial optimization problem: $\binom{n}{k}$ sampling sets must still be checked, which is impractical even for moderately small $n$. In fact, due to the irregularity of the domain of graph signals, sampling set selection is NP-hard in general. It is straightforward to see that it is equivalent to the sensor placement or forward regression problems in [16, 17, 18, 19], so that the typical reduction from set cover applies [40].

In the following sections, we address this issue in two ways. First, we derive universal performance bounds that hold for all sampling sets (Section III). These bounds can therefore be used to evaluate the quality of a sampling set or selection heuristic a posteriori. Second, we study the greedy sampling algorithm and provide near-optimal a priori guarantees based on the concept of approximate submodularity (Section IV). Special cases of these results that considered real-valued homeoscedastic signal and noise ($\bm{\Lambda}=\sigma_{x}^{2}\bm{I}$ and $\bm{\Lambda}_{w}=\sigma_{w}^{2}\bm{I}$) and no transformation of the graph signal ($\bm{H}=\bm{I}$) can be found in [32, 24].

III Universal Bounds on the Interpolation MSE

In this section, we derive interpolation performance bounds that hold for all $\mathcal{S}$. These universal bounds can be used to inform the sampling set selection by (i) describing how different factors influence the reconstruction performance and (ii) gauging the quality of sampling set instances. The main result of this section is presented below.

Theorem 1**.**

Let $\bm{x}$ be a $\mathcal{K}$-ssparse stationary graph signal and $\bm{y}=\bm{x}+\bm{w}$ be its noisy observations. Take $\hat{\bm{z}}^{\star}=\bm{L}^{\star}\bm{y}_{\mathcal{S}}$ to be the minimum MSE linear interpolation of $\bm{z}=\bm{H}\bm{x}$ based on a sampling set $\mathcal{S}$ given zero-mean circular priors on the signal and noise such that $\operatorname{\mathbb{E}}\bm{x}\bm{x}^{H}=\bm{V}_{\mathcal{K}}\bm{\Lambda}\bm{V}_{\mathcal{K}}^{H}$ and $\operatorname{\mathbb{E}}\bm{w}\bm{w}^{H}=\bm{\Lambda}_{w}$. For $\bm{W}=\bm{V}_{\mathcal{K}}^{H}\bm{H}^{H}\bm{H}\bm{V}_{\mathcal{K}}\succ 0$, the reconstruction error $\textup{MSE}(\mathcal{S})=\operatorname{\mathbb{E}}\left\|\bm{z}-\hat{\bm{z}}^{\star}\right\|^{2}=\operatorname{Tr}[\bm{K}^{\star}(\mathcal{S})]$ is bounded by

[TABLE]

where $\bar{\ell}_{m}$ is the sum of the $m$ largest weighted structural SNRs $\ell_{i}=\lambda_{w,i}^{-1}\left\|\bm{v}_{i}\right\|_{\bm{W}^{-1}}^{2}$, with \bm{V}_{\mathcal{K}}=[\begin{array}[]{ccc}\bm{v}_{1}&\cdots&\bm{v}_{N}\end{array}]^{H} and $\left\|\bm{x}\right\|_{\bm{A}}^{2}=\bm{x}^{H}\bm{A}\bm{x}$. Explicitly, $\bar{\ell}_{m}=\max_{\mathcal{X}:{\left|\mathcal{X}\right|}=m}\sum_{j\in\mathcal{X}}\ell_{j}$.

Proof.

Start with the upper bound that is achieved for an empty sampling set, i.e., for $\mathcal{S}=\{\}$. Indeed, recall from Lemma 1 that $\bm{K}^{\star}$ is a monotone decreasing set function, i.e., it achieves its maximum for the empty set. Thus, it holds that $\bm{K}(\hat{\bm{x}}^{\star})\preceq\bm{H}\bm{V}_{\mathcal{K}}\bm{\Lambda}\bm{V}_{\mathcal{K}}^{H}\bm{H}^{H}$, from which the upper bound in (14) follows by the monotonicity of the trace operator [39].

To obtain the lower bound, start by using (9) to get

[TABLE]

where we used the circular commutation property of the trace. Then, since the trace of a matrix is the sum of its eigenvalues, the arithmetic/harmonic means inequality can be used to get, for any $N\times N$ positive-definite matrix $\bm{X}$,

[TABLE]

with equality if and only if $\bm{X}=\gamma\bm{I}$, $\gamma>0$ [30]. Since $\bm{W}\succ 0$, the matrix in (15) is positive-definite and we have

[TABLE]

which from the commutation property of the trace gives

[TABLE]

where $\left\|\bm{x}\right\|_{\bm{A}}^{2}=\bm{x}^{H}\bm{A}\bm{x}$ is a weighted norm. Finally, replacing the sum in the denominator by its maximum value $\bar{\ell}_{{\left|\mathcal{S}\right|}}$ gives the desired lower bound in (14).∎

The bounds in (14) were derived from a Bayesian perspective, so that the expectation in the MSE is taken over realizations of the signal and noise and the bounds hold for all sampling sets $\mathcal{S}\subseteq\mathcal{V}$. Also, it is worth noting that (14) depends only on statistics of the graph signal ($\bm{\Lambda}$, $\bm{\Lambda}_{w}$, and $\mathcal{K}$), the transform ($\bm{H}$), the structural properties of the underlying graph ($\bm{V}$), and the sampling set size (${\left|\mathcal{S}\right|}$). These are all quantities known a priori, i.e., before the sampling occurs.

As expected, (14) decreases with the sampling set size. The rate of decay, however, depends on the weighted structural SNRs $\{\ell_{i}\}$. These quantities represent the relation between the signal of interest and the noise at each node, taking into account the structure of the graph and the subspace of interest [$\operatorname*{colspan}(\bm{V}_{\mathcal{K}})$]. Moreover, they are related to statistical estimates such as the leverage score and the Mahalanobis distance in regression. In a sequential sampling scheme, their value can be used to inform whether a new sample is worth acquiring by bounding the MSE improvement. A bound on the decay rate can also be obtained using the fact that $\bar{\ell}_{m}\leq m\ell_{\text{max}}$ for $\ell_{\text{max}}=\max_{i}\ell_{i}$. Then, if the sampling set is chosen so as to uniquely determine the graph signal, i.e., ${\left|\mathcal{S}\right|}\geq{\left|\mathcal{K}\right|}$, (14) reduces to

[TABLE]

where we used the fact that $\operatorname{Tr}\left[(\bm{W}\bm{\Lambda})^{-1}\right]\leq{\left|\mathcal{K}\right|}\lambda_{\text{min}}(\bm{W}\bm{\Lambda})^{-1}$. It is clear from (16) that the reconstruction error increases linearly with the bandwidth of the graph signal, which is a fundamental limitation for large dimensional signals. It also shows the importance of working with low bandwidth signals and, consequently, of appropriately identifying the signal’s underlying graph.

Although these observations give insights into graph signal interpolation, one of the main motivation behind Theorem 1 is addressing the issue of sampling set selection. Towards this end, we propose the following corollary:

Corollary 1**.**

For any graph signal and its interpolation as in Theorem 1, all sampling set $\mathcal{S}$ for which $\textup{MSE}(\mathcal{S})\leq\eta$ satisfy

[TABLE]

Since $\bar{\ell}_{{\left|\mathcal{S}\right|}}\leq{\left|\mathcal{S}\right|}\ell_{\text{max}}$, it also holds that

[TABLE]

Corollary 1 gives a lower bound on the number of samples needed to achieve a desired MSE. It is worth noting that this bound does not inform whether the specific MSE value $\eta$ is achievable with less than $n$ samples. Whenever it is not, Corollary 1 holds trivially. From (18), note that the minimum number of samples increases as the MSE decreases. Moreover, although (18) suggest that the sample set size required to achieve a certain MSE grows with $\mathcal{O}({\left|\mathcal{K}\right|}^{2})$, it is not necessarily the case. Indeed, recall that $\ell_{\text{max}}$ is a function of ${\left|\mathcal{K}\right|}$ (through $\left\|\bm{v}_{i}\right\|$ and $\bm{V}_{\mathcal{K}}$). Still, as in the noiseless case, the signal bandwidth is a dominating factor in the determination of the minimum sampling set size.

Although (18) characterizes the overall behavior of the sampling set size, it is not informative in practice because it largely underestimates ${\left|\mathcal{S}\right|}$. On the other hand, (17) yields a tighter bound which can be used, together with (14), to evaluate a sampling set or sampling technique for direct reconstruction of a graph signal ($\bm{H}=\bm{I}$). Indeed, Figures 1 and 2 compares (14) and (17) to the minimum interpolation MSE and optimal set size, found by exhaustive search, for three graph models ($n=20$): Erdős-Rényi, preferential attachment, and a random undirected graph with weights uniformly distributed in $[0,1]$ (see details of these models in Section V-A). The graph signal is assumed to be homeoscedastic with $\bm{\Lambda}=\bm{I}$ and $\bm{\Lambda}_{w}=\sigma_{w}^{2}\bm{I}$, $\sigma_{w}^{2}=10^{-2}$. Note that the bounds are conservative for ${\left|\mathcal{S}\right|}<{\left|\mathcal{K}\right|}$, but become tighter as ${\left|\mathcal{S}\right|}$ increases. This is because the inequality used to derive (14) becomes tighter as the eigenvalues of $\bm{K}^{\star}$ become more similar.

When $\bm{H}$ is arbitrary and variance of signal and noise can vary across nodes, the eigenvalues of $\bm{K}^{\star}$ can become different from each another and deteriorate the bound in (14). This is illustrated in Figures 3 and 4, where $\bm{H}$ was taken as a $30\times 20$ matrix whose entries are zero-mean unit variance Gaussian random variables, $\bm{\Lambda}=\bm{I}$, and the noise variance was uniformly distributed in $[10^{-3},10^{-1}]$.

Remark 2**.**

Bounds on the interpolation MSE of graph sampling techniques have also been derived in [10, 13]. These works consider randomized sampling set selection schemes, including uniform and leverage score sampling, and derive performance bounds on the optimal sampling distributions and interpolation error. The bounds in Theorem 1 and Corollary 1 differ from those in [10, 13] in that the latter take the spectrum of the graph signal to be deterministic and the sampling to be random. Thus, these bounds hold in expectation over different sampling realizations for a specific randomized strategy. The bounds in (14) hold in expectation over realizations of the signal and noise and give a worst-case performance bound for any—possibly randomized—sampling strategy.

IV Near-Optimal Sampling Set Selection

Although the bounds from Section III can be used to evaluate specific sampling set instances, they do not provide performance guarantees for any sampling strategy. To do that, this section studies a specific sampling scheme, namely greedy sampling set selection, and derives near-optimality results that hold for all problem instances.

Greedy sampling set selection is pervasive in GSP and has proven successful in many applications [10, 11, 9, 21, 12, 15]. This is illustrated in Figure 5 which uses the bounds derived in (17) to assess the quality of sampling sets obtained by greedily minimizing the MSE (see Algorithm 2) on larger instances ($n=1000$) of the three random graph models found in Figures 1 to 4. Note that the final greedy sampling set size remains within 10% of the lower bound in these realizations.

Despite strong empirical evidences, typical performance guarantees for greedy search do not hold for greedy sampling set selection. Indeed, the well-established result from [41] states that greedy minimization (Algorithm 1) yields guaranteed near-optimal results for monotonically decreasing and supermodular set functions. The MSE, however, is not supermodular in general. This can be seen from [18, Thm. 2.4] and the fact that $f(t)=t^{-2}$ is not operator antitone [39]. Thus, although greedily minimizing the MSE appears to work in practice, there has yet to be a theoretical justification for it. The following sections bridge this gap by expanding the notion of approximate supermodularity introduced in [24] and updating the performance bound from [41] for this class of functions. This novel framework then allows near-optimality bounds to be derived for the MSE.

IV-A Approximate supermodularity and greedy minimization

Supermodularity (and its dual submodularity) encodes the “diminishing returns” property of certain functions that leads to bounds on the suboptimality of their greedy minimization [41]. Well-known supermodular functions include the rank, $\log\det$, or Von-Neumann entropy of a matrix [23]. Still, supermodularity is a stringent condition. In particular, it does not hold for the MSE in (13). To provide suboptimality bounds for its greedy minimization, we therefore define the concept of approximate supermodularity.

A set function $f:2^{\mathcal{V}}\to\mathbb{R}$ is $\alpha$-supermodular if for all sets $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{V}$ and all $u\notin\mathcal{B}$ it holds that

[TABLE]

for $\alpha\geq 0$. We say $f$ is $\alpha$-submodular if $-f$ is $\alpha$-supermodular. For $\alpha\geq 1$, (19) is equivalent to the traditional definition of supermodularity, in which case we refer to the function simply as supermodular/submodular [23]. For $\alpha\in[0,1)$, however, $f$ is said to be approximately supermodular/submodular. Notice that (19) always holds for $\alpha=0$ if $f$ is monotone decreasing. Indeed, $f\left(\mathcal{A}\cup\{u\}\right)-f\left(\mathcal{A}\right)\leq 0$ in this case. Thus, $\alpha$-supermodularity is only of interest when $\alpha$ takes the largest value for which (19) holds, i.e.,

[TABLE]

Before proceeding, it is worth noting that $\alpha$ is related to the submodularity ratio from [17]. It is, however, more amenable to give explicit bounds on its value (see Section IV-B). The submodularity ratio bounds derived in [17] depend on the minimum sparse eigenvalue of a matrix, which cannot be evaluated efficiently. Due to the relation between $\alpha$ and the submodularity ratio, it is not surprising that similar near-optimal bounds hold for $\alpha$-supermodular functions.

Theorem 2**.**

Let $f^{\star}=f(\mathcal{S}^{\star})$ be the optimal value of the problem

[TABLE]

and $\mathcal{G}_{\ell}$ be the set obtained by applying Algorithm 1. If $f$ is (i) monotone decreasing and (ii) $\alpha$-supermodular, then

[TABLE]

where $f(\{\})$ is the value of function for the empty set. If $f$ is normalized, i.e., $f(\{\})=0$, (22) reduces to

[TABLE]

Proof.

Using the fact that $f$ is monotone decreasing, it holds for every set $\mathcal{G}_{j}$ that

[TABLE]

Using a telescopic sum then gives

[TABLE]

where $\mathcal{T}_{i}=\mathcal{G}_{j}\cup\{s_{1}^{\star},\dots,s_{i}^{\star}\}$ and $s_{i}^{\star}$ is the $i$-th element of $\mathcal{S}^{\star}$. Since $f$ is $\alpha$-supermodular and $\mathcal{G}_{j}\subseteq\mathcal{T}_{i}$ for all $i$, the incremental gains in the summation in (23) can be bounded using (19) to get

[TABLE]

Finally, given that $\mathcal{G}_{j+1}=\mathcal{G}_{j}\cup\{u\}$ is chosen to minimize $f(\mathcal{G}_{j+1})$ (see Algorithm 1),

[TABLE]

To obtain a recursion, let $\delta_{j}=f(\mathcal{G}_{j})-f(\mathcal{S}^{\star})$ so that (24) becomes

[TABLE]

Noting that $\delta_{0}=f(\{\})-f(\mathcal{S}^{\star})$, we can solve this recursion to get

[TABLE]

Using the fact that $1-x\leq e^{-x}$ yields (22). ∎

Theorem 2 bounds the relative suboptimality of the greedy solution to problem (21) when $f$ is decreasing and $\alpha$-supermodular. Under these conditions, it guarantees a minimum improvement of the greedy solution over the empty set. What is more, it quantifies the effect of relaxing the supermodularity hypothesis in (19). Indeed, when $f$ is supermodular ($\alpha=1$) and the greedy search in Algorithm 1 is repeated $k$ times ($\ell=k$), we recover the $e^{-1}\approx 0.37$ guarantee from [41]. On the other hand, if $f$ is not supermodular ($\alpha<1$), (22) shows that the same $37\%$ guarantee can be obtained by greedily selecting a set of size $\alpha^{-1}k$. Thus, $\alpha$ not only quantifies how much $f$ violates supermodularity, but also gives a factor by which a solution set must increase to maintain supermodular near-optimality. In other words, it measures the constraint violation needed to recover the $37\%$ guarantee. It is worth noting that, as with the original bound in [41], (22) is not tight and that better results are common in practice (see Section V-A).

In the sequel, we show that $\textup{MSE}(\mathcal{S})$ is a monotone decreasing and $\alpha$-supermodular function of the sampling set $\mathcal{S}$. We also provide an explicit lower bound on $\alpha$ as a function of the SNR. This result simultaneously provide near-optimal performance guarantees based on Theorem 2 and sheds light on why greedy algorithms have been so successful in GSP applications.

IV-B Near-optimality of greedy sampling set selection

The main result of this section is:

Theorem 3**.**

Let $\mathcal{S}^{\star}$ be the solution of (13) and $\mathcal{G}_{\ell}$ be the result of applying Algorithm 1 for $f(\mathcal{S})=\textup{MSE}(\mathcal{S})$. Then,

[TABLE]

where

[TABLE]

for $\mu_{\textup{min}}\leq\lambda_{\textup{min}}\left[\bm{\Lambda}^{-1}\right]$, $\mu_{\textup{max}}\geq\lambda_{\textup{max}}\left[\bm{\Lambda}^{-1}+\bm{V}_{\mathcal{K}}^{H}\bm{\Lambda}_{w}^{-1}\bm{V}_{\mathcal{K}}\right]$, and $\kappa_{2}(\bm{W})$ is the 2-norm condition number of $\bm{W}$. Assuming $\bm{\Lambda}=\sigma_{x}^{2}\bm{I}$ and $\bm{\Lambda}_{w}=\sigma_{w}^{2}\bm{I}$, (25) reduces to

[TABLE]

Theorem 3 establishes that a near-optimal solution to the sampling set selection problem in (13) can be obtained efficiently using greedy search. Though strong empirical evidence exists that greedily minimizing the MSE yields good results in contexts such as regression, dictionary learning, and graph signal processing [17, 10, 11, 9, 21, 12], this result is counter-intuitive given that the MSE is not supermodular in general. For instance, restrictive and often unrealistic conditions on data distribution are required to obtain supermodularity in the context of regression [17].

Theorem 3 therefore reconciles the empirical success of greedy sampling set selection and the non-supermodularity of the MSE by bounding the suboptimality of greedy sampling. In particular, (26) gives a simple bound on $\alpha$ in terms of the SNR and the condition number of $\bm{W}$ that gives clear insights into its behavior. Indeed, as $\gamma\to\infty$ and we approach the noiseless case, $\alpha\to 0$. This is expected as in the noiseless case almost every set of size ${\left|\mathcal{K}\right|}$ achieves perfect reconstruction, so that the choice of sampling nodes is irrelevant. On the other hand, $\alpha\to 1$ as $\gamma\to 0$, i.e., the MSE becomes closer to supermodular as the SNR decreases. Given that reconstruction errors are small for high SNR, Theorem 3 guarantees that greedy sampling performs well when it is most needed. Similar trends can be observed in the more general setting of (25). These observations are illustrated in Figure 6 that compares the bound in (26) to the true value of $\alpha$ for the MSE (found by exhaustive search) in $100$ realizations of random graphs (see Section V-A for details).

Theorem 3 stems directly from Theorem 2 and the following characterizations of the MSE function:

Lemma 2**.**

The scalar set functions $\textup{MSE}(\mathcal{S})=\operatorname{Tr}[\bm{K}^{\star}(\mathcal{S})]$ is (i) monotone decreasing and (ii) $\alpha$-supermodular with

[TABLE]

where $\mu_{\textup{min}}\leq\lambda_{\textup{min}}\left[\bm{\Lambda}^{-1}\right]$, $\mu_{\textup{max}}\geq\lambda_{\textup{max}}\left[\bm{\Lambda}^{-1}+\bm{V}_{\mathcal{K}}^{H}\bm{\Lambda}_{w}^{-1}\bm{V}_{\mathcal{K}}\right]$, and $\kappa_{2}(\bm{W})$ is the 2-norm condition number of $\bm{W}$ [30].

Lemma 3**.**

Assuming $\bm{\Lambda}=\sigma_{x}^{2}\bm{I}$ and $\bm{\Lambda}_{w}=\sigma_{w}^{2}\bm{I}$, the set functions $\textup{MSE}(\mathcal{S})=\operatorname{Tr}[\bm{K}^{\star}(\mathcal{S})]$ is (i) monotone decreasing and (ii) $\alpha$-supermodular with

[TABLE]

where $\kappa_{2}(\bm{W})$ is the 2-norm condition number of $\bm{W}$ [30].

The proof of Lemma 2 is deferred to Appendix A. Here, we proceed with the proof of Lemma 3 after stating a pertinent remark.

Remark 3**.**

Since the MSE is not supermodular, it is common to see surrogate supermodular figures of merit used instead, specially in statistics and experiment design [23, 17, 18, 16]. In particular, the log-determinant $\log\det[\bm{K}^{\star}(\mathcal{S})]$ is a common alternative to the objective $\textup{MSE}(\mathcal{S})=\operatorname{Tr}[\bm{K}^{\star}(\mathcal{S})]$ used in (13). This is justified because the $\log\det[\bm{K}^{\star}(\mathcal{S})]$ is proportional to the volume of the confidence ellipsoids of the estimate when the data is Gaussian [42, 18]. This choice of objective is also common in the sensor placement literature due to its relation to information theoretic measures, such as entropy and mutual information [16]. By replacing the trace operator in (13) by the $\log\det$, the problem becomes a supermodular function minimization that can be efficiently approximated using greedy search, as shown in [24, 15]. We remark that minimizing the $\log\det$ of the error covariance matrix and the MSE are not equivalent problems.

IV-C Proof of Lemma 3

Start by noticing that part (i) stems directly from Lemma 1. Indeed, the monotonicity of the trace implies that $\bm{X}\succeq\bm{Y}\Rightarrow\operatorname{Tr}(\bm{X})\geq\operatorname{Tr}(\bm{Y})$, for any PSD matrices $\bm{X}$ and $\bm{Y}$.

Then, to obtain part (ii), use the homeoscedasticity assumption to rewrite (9) as

[TABLE]

where $\bm{Z}(\mathcal{S})=\bm{I}+\gamma\sum_{i\in\mathcal{S}}\bm{v}_{i}\bm{v}_{i}^{H}$ and $\gamma=\sigma_{x}^{2}/\sigma_{w}^{2}$ is the SNR. Then, proceed to obtain a closed form expression for the increments in (20) by using (29) to get

[TABLE]

From the matrix inversion lemma [30], this expression reduces to

[TABLE]

which using the commutation property of the trace yields

[TABLE]

From (30), the expression for $\alpha$ in (20) becomes

[TABLE]

To bound (31), first notice that for any set $\mathcal{X}\subseteq\mathcal{V}$

[TABLE]

where $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ denote the minimum and maximum eigenvalues of a matrix. These bounds are achieved for the empty set and $\mathcal{V}$, respectively. Then, using the Rayleigh quotient inequalities [30]

[TABLE]

we get that (31) is bounded by

[TABLE]

To simplify this expression, let $\sigma_{i}(\bm{A})$ denote the $i$-th singular value of $\bm{A}$ and recall that for $\bm{A}\succeq 0$, $\sigma_{i}^{2}(\bm{A})=\lambda_{i}(\bm{A}\bm{A}^{H})$. Thus, we can write $\lambda_{i}\left[\bm{Z}(\mathcal{A})^{-1}\bm{W}\bm{Z}(\mathcal{A})^{-1}\right]=\sigma_{i}^{2}\left[\bm{Z}(\mathcal{A})^{-1}\bm{W}^{1/2}\right]$, which is well-defined since $\bm{W}\succeq 0$. Using the fact that $\sigma_{\text{max}}(\bm{A}\bm{B})\leq\sigma_{\text{max}}(\bm{A})\sigma_{\text{max}}(\bm{B})$ and $\sigma_{\text{min}}(\bm{A}\bm{B})\geq\sigma_{\text{min}}(\bm{A})\sigma_{\text{min}}(\bm{B})$ [43, Thm. 9.H.1, p. 338], we obtain

[TABLE]

where $\kappa_{2}(\bm{W})=\lambda_{\text{max}}(\bm{W})/\lambda_{\text{min}}(\bm{W})$ is the 2-norm condition number of $\bm{W}$ [30].

Finally, to obtain the expression in Lemma 3, notice that (33) is decreasing with respect to $\left\|\bm{v}_{u}\right\|_{2}^{2}$. Indeed, since $\kappa_{2}\geq 1$ and $\gamma\geq 0$,

[TABLE]

Given that $\bm{v}_{u}$ is a row of $\bm{V}_{\mathcal{K}}$, i.e., it is composed of a subset of elements from a unit vector, $\left\|\bm{v}_{u}\right\|_{2}^{2}\leq 1$ and we obtain the result in (28).∎

V Numerical Examples and Applications

Before proceeding with the simulations, the complexity issue of greedy sampling set selection must be addressed. The greedy search in Algorithm 1 requires $n\ell c_{f}$ operations, where $c_{f}$ is the cost of evaluating the objective $f$. As it is, problem (13) has $c_{f}=\mathcal{O}({\left|\mathcal{K}\right|}^{3})$. It can, however, be reduced using the matrix inversion lemma [30].

Indeed, start by noticing that the first step of the greedy approximation of problem (13) involves finding (see Algorithm 1)

[TABLE]

which, using the definition of $\bm{K}^{\star}$ in (9) and the circular commutation property of the trace, requires the evaluation of

[TABLE]

where once again $\bm{W}=\bm{V}_{\mathcal{K}}^{H}\bm{H}^{H}\bm{H}\bm{V}_{\mathcal{K}}$. Letting $\bm{K}^{\star}_{j}=\bm{K}^{\star}(\mathcal{G}_{j})$ and using the matrix inversion lemma, we can reduce the update of $\bm{K}^{\star}$ to

[TABLE]

From linearity, it is then straightforward to see that finding the minimum of the trace of (34) is equivalent to finding the maximum of

[TABLE]

The greedy sampling set selection procedure obtained by leveraging (34) and (35) is presented in Algorithm 2. This algorithm now requires only $\mathcal{O}(n\ell{\left|\mathcal{K}\right|}^{2})$ operations.

V-A Simulations

In this section, we start by evaluating the performance greedy sampling set selection (Algorithm 2). For comparison, we also display the results obtained by the uniform and leverage score randomized methods from [13] and a deterministic heuristic based on sampling nodes with the highest leverage score ($\left\|\bm{v}_{i}\right\|_{2}^{2}$). In the following examples, we use undirected graphs generated using the Erdős-Rényi model, in which an edge is placed between two nodes with probability $p=0.2$; the preferential attachment model [44], in which nodes are added one at a time and connected to a node already in the graph with probability proportional to its degree; and a random undirected graph, obtained by assigning a weight to all possible edges uniformly at random from $[0,1]$.

The figure of merit in the following simulations is the relative suboptimality from (22). Since it depends on the optimal sampling set which needs to be determined by exhaustive search, we focus on graphs with $n=20$ nodes. Bandlimited graph signals are generated by taking $\bm{V}_{\mathcal{K}}$ in (3) to be the eigenvectors of the graph adjacency matrix relative to the five eigenvalues with largest magnitude (${\left|\mathcal{K}\right|}=5$). The random vectors $\bar{\bm{x}}$ in (3) and $\bm{w}$ in (5) are realizations of zero-mean Gaussian random variables with covariance matrices $\bm{\Lambda}=\bm{I}$ and $\bm{\Lambda}_{w}=\sigma_{w}^{2}\bm{I}$, where $\sigma_{w}^{2}$ is varied to obtain different SNRs. The transform in (4) is taken to be the identity ($\bm{H}=\bm{I}$) and the sampling set size is chosen as $\ell={\left|\mathcal{K}\right|}=5$.

Figures 7 and 8 display histograms of the relative suboptimality for $1000$ realizations of graphs and graph signals with $\sigma_{w}^{2}=10^{2}$ ($\text{SNR}=-20$ dB) and $\sigma_{w}^{2}=10^{-2}$ ($\text{SNR}=20$ dB), respectively. As predicted by Theorem 3, greedy sampling set selection performs better in low SNR environments, where the optimal sampling set was obtained in more than $95$% of the realizations. Nevertheless, even in high SNRs, it found the optimal sampling set almost half of the time. In fact, note that Algorithm 2 typically performs much better than the bounds in Theorem 3 (see details in Fig. 8). For comparison, results for greedily optimizing $\log\det\left[\bm{K}^{\star}(\mathcal{S})\right]$, a supermodular function, are shown in Figure 9. Although the MSE is $\alpha$-supermodular with $\alpha<1$, the relative suboptimality obtained by using both cost functions is comparable.

It is worth noting that, although the deterministic leverage score ranking technique often yields good results, there are advantages to greedy sampling set selection, specially for higher SNR. The randomized sampling schemes, on the other hand, do not perform as well for single problem instances. To be fair, these methods are more appropriate when several sampling sets of the same graph signal are considered. Indeed, the performance measures in [13] hold in expectation over sampling realizations.

Evaluating the relative suboptimality for larger graphs is untractable. However, since these sampling set selection techniques build the sampling set sequentially, we can assess their performance in terms of the sampling set size required to obtain a given MSE reduction. Figure 10 displays the distribution of the sampling set size required to achieve a $90\%$ reduction in the MSE with respect to the empty set. The plots are obtained from $1000$ graphs and signals realizations with $n=100$ nodes, $\bm{V}_{\mathcal{K}}$ in (3) composed the eigenvectors relative to the seven eigenvalues with largest magnitude (${\left|\mathcal{K}\right|}=7$), and $\sigma_{w}^{2}=10^{-2}$. Although Theorem 3 estimates that Algorithm 2 requires sets considerably larger to recover the same near-optimal guarantees as supermodular functions, greedy sampling obtained a sampling set of size exactly ${\left|\mathcal{K}\right|}$ in more than $50\%$ of the realizations. Moreover, as noted in [13], we can now see that leverage score sampling has similar performance to uniform sampling for Erdős-Rényi graphs, but gives better results for the preferential attachment model.

V-B Application: Subsampled Kernel PCA

Kernel PCA is a nonlinear version of PCA [27] that also identifies data subspaces by truncating the eigenvalue decomposition (EVD) of a Gram matrix $\bm{\Phi}$. However, whereas PCA uses the empirical covariance matrix, kPCA constructs $\bm{\Phi}$ by evaluating inner products between data points in a higher dimensional space $\mathbb{F}$ known as the feature space. Since the map $\varphi:\mathbb{R}^{m}\to\mathbb{F}$ can be nonlinear and $\mathbb{F}$ typically has infinite dimensionality, kPCA results in richer subspaces than PCA [27, 28, 45].

Naturally, the dimensionality of $\mathbb{F}$ poses a challenge for constructing the Gram matrix. This problem is addressed using the so called kernel trick [27, 28, 45]. A kernel is a function $\kappa$ that allows the inner product in $\mathbb{F}$ to be evaluated directly from vectors in $\mathbb{R}^{m}$, i.e., $\kappa(\bm{r},\bm{s})=\langle\varphi(\bm{r}),\varphi(\bm{s})\rangle_{\mathbb{F}}$. We can use $\kappa$ to construct $\bm{\Phi}$ from a training set $\{\bm{u}_{i}\}_{i=1,\dots,n}$, $\bm{u}_{i}\in\mathbb{R}^{m}$, as in

[TABLE]

Kernel PCA identifies the data subspace as the span of the first $k$ eigenvectors of $\bm{\Phi}$, i.e., as $\operatorname*{colspan}(\bm{V}_{\mathcal{K}})$, where $\bm{\Phi}=\bm{V}\bm{\Lambda}\bm{V}^{H}$ is the EVD of $\bm{\Phi}$ with eigenvalues in decreasing order and $\mathcal{K}=1,\dots,k$. Using the representer’s theorem [45], any data point $\bm{y}$ can be projected onto this subspace by

[TABLE]

The projection in (37) requires $\Theta(kn)$ operations and $n$ KEs, making this method impractical for large data sets even if the dimension $k$ of the subspace of interest is small. Indeed, although the training phase in (36) is usually performed offline, (37) needs to be evaluated during the operation phase for every new data point. In [46], this issue was addressed by using a Gaussian generative model for $\bm{\Phi}$ and showing that its maximum likelihood estimate depends only on a subset of the $\bm{u}_{i}$. Another approach is to impose sparsity on $\bm{V}$ a priori so that it depends only on a reduced number of training points [28]. Alternatively, one can find a representative subset of the training data and apply kPCA to that subset [47]. The issue with the latter method is that finding a good data subset is known to be a hard problem [25, 26]. In fact, it is related to the problem of sampling set selection in GSP.

Indeed, since we used the same notation as in Section II, formulating kPCA in the context of GSP is straightforward. Let the graph $\mathbb{G}$ have adjacency matrix $\bm{A}=\bm{\Phi}$, which is symmetric and normal, so that (37) has the form of a (partial) graph Fourier transform (3). In other words, (37) can be interpreted as enforcing graph signals of the form $\tilde{\bm{y}}$ to be bandlimited on $\bm{\Phi}$. Thus, we can apply the sampling and interpolation theory from GSP to put forward a subsampled kPCA.

Based on the guarantees given in Section IV, we use greedy search to obtain a sampling set $\mathcal{S}$ and use the interpolation techniques from Section II-A to recover $\tilde{\bm{y}}$ from its samples as in

[TABLE]

Then, (37) and (38) yield

[TABLE]

Notice that $\bm{P}$ is now $k\times{\left|\mathcal{S}\right|}$, so that the projection in (39) only takes $\Theta(k{\left|\mathcal{S}\right|})$ operations and ${\left|\mathcal{S}\right|}$ KEs, leading to a considerable complexity reduction (${\left|\mathcal{S}\right|}/n$) over the direct projection in (37). Moreover, kPCA is typically used for dimensionality reduction prior to regression or classification, so that we are actually interested in a linear transformation of $\bar{\bm{y}}$. Subsampled kPCA can account for this case by properly choosing $\bm{H}$ in (4). It is worth noting that contrary to [47], the full dataset is used during the training stage to obtain $\bm{V}_{\mathcal{K}}$. However, once $\bm{P}$ is determined, only the subset $\mathcal{S}$ is required.

In the sequel, we illustrate this method in a face recognition application using the faces94 data set [48]. It contains $20$ pictures ($200\times 180$) of $152$ individuals which were converted to black and white and normalized so that the value of each pixel is in $[-1,1]$. A training set is obtained by randomly choosing $14$ images for each individual ($70$% of the data set) and the remaining pictures are used for testing. In this application, we use a polynomial kernel of degree $d=2$ [45] and a one-against-one multiclass support vector machine (SVM) classifier, in which an SVM is trained for each pair of class and the classification is obtained by majority voting (see [49] for details on this scheme). Finally, note that since images in both training and testing sets come from the same data set there is no observation noise $\bm{w}$. Still, $\sigma_{w}^{2}$ can be used to regularize the matrix inversions in (8) and (9) [29].

Figure 11 shows the misclassification percentage on the test set as a function of the number of components ($k$) for both PCA and kPCA. Note that kPCA can achieve the same performance as PCA with less components. The results of using greedy subsampled kPCA are shown in Figure 12 and clearly illustrate the trade-off between complexity and performance: as the sampling set size increases, the classification errors decrease. However, since misclassification is a nonlinear function of the MSE, it may be advantageous to use more components instead of increasing the sampling set. For instance, kPCA requires $k=7$ components to achieve a misclassification of $1$%, so that evaluating the direct projection in (37) takes $2128$ KEs and $29785$ operations. Greedy subsampled kPCA, on the other hand, can achieve the same performance with $k=12$ components and ${\left|\mathcal{S}\right|}=33$, i.e., $33$ KEs and $780$ operations, a complexity reduction of more than $97$%. Nevertheless, using $7$ components, greedy subsampled kPCA would require $\mathcal{S}$ to be almost the full training set.

Naturally, the method in (36)–(39) is not restricted to sampling sets obtained greedily. Thus, in Figures 13 and 14 we compare greedy sampling to the other methods from Section V-A now based on their misclassification performance for $12$ and $15$ components. We omit the deterministic leverage score results because it performed consistently worst than the other methods. The average misclassification rates for the randomized schemes are from $1$ to $15\%$ higher than those of greedy sampling. Although some realizations yield good classification, their performance varies a lot, especially for smaller sampling sets. Comparing Figures 13 and 14, it is ready that in this application the difference between uniform and leverage score sampling becomes less significant as the number of components increases.

Remark 4**.**

Although this section discussed kPCA, the same argument applies to the classical PCA. It is therefore straightforward to derive an analog subsampled PCA technique using (36)–(39).

VI Conclusion

This work provided a solution to graph signal sampling problems by addressing the issue of sampling set selection in two ways. First, it derived universal bounds on the interpolation MSE (Theorem 1) which allow the quality of any sampling set or sampling heuristic to be evaluated by gauging how close their reconstruction performance is to the lower bound. Second, it provided near-optimality results for greedy MSE minimization (Theorem 3), demonstrating that greedy sampling set selection is an effective sampling scheme, justifying its empirical success in the literature. The strength of Theorem 3 is that it gives a worst-case result: there exists no graph or graph signal for which the relative suboptimality of greedy sampling is worst than $e^{-\alpha}$. In fact, greedy sampling typically performs much better and consistently across graph signal realizations. In contrast, although randomized sampling schemes can be effective, their performance can vary widely across realizations. We should note that given the generality of Theorems 2 and 3, it is likely that they can be applied beyond GSP for sensor placement, experimental design, and variable selection. Moreover, despite the MSE’s ubiquity in signal processing, other performance metrics are sometimes more appropriate and we foresee that the theory from this paper can be extended to these cases. In particular, we believe that the concept of approximate submodularity can be used to provide guarantees for the greedy minimization of other non-supermodular functions.

Appendix A Proof of Lemma 2

Proof.

Once again, part (i) is a corollary of Lemma 1. To prove part (ii), we proceed as for Lemma 3. However, we now let $\bm{Z}(\mathcal{A})=\bm{\Lambda}^{-1}+\sum_{i\in\mathcal{A}}\lambda_{w,i}^{-1}\bm{v}_{i}\bm{v}_{i}^{H}$ so that again the increment in (20) reads

[TABLE]

which using the matrix inversion simplifies to

[TABLE]

Using this expression, the expression for $\alpha$ in (20) becomes

[TABLE]

Note that, although similar, (40) is not the same as (31).

We now bound (40) by noticing that for any set $\mathcal{X}\subseteq\mathcal{V}$

[TABLE]

Thus, using the Rayleigh quotient inequalities leads to

[TABLE]

which can be simplified using the same singular value bounds as in Lemma 3 [43, Thm. 9.H.1, p. 338] to yield

[TABLE]

where again $\kappa_{2}(\bm{W})=\lambda_{\text{max}}(\bm{W})/\lambda_{\text{min}}(\bm{W})$ is the 2-norm condition number of $\bm{W}$. To obtain the expression in (27), notice that (41) is decreasing with respect to $\left\|\bm{v}_{u}\right\|_{2}^{2}$ and $\lambda_{w,u}^{-1}$. Indeed,

[TABLE]

are both non-positive because $0<\mu_{\text{min}}\leq\mu_{\text{max}}$ and $\kappa_{2}(\bm{W})\geq 1$ [30]. We then use the fact that $\left\|\bm{v}_{u}\right\|_{2}^{2}\leq 1$ to get (27). ∎