Low-complexity Graph Sampling with Noise and Signal Reconstruction via Neumann Series

TL;DR

This paper introduces a low-complexity graph sampling method using Neumann series and fast Fourier transforms, enabling efficient node selection and robust signal reconstruction in noisy and dynamic networks.

Contribution

It proposes a novel greedy sampling strategy based on an approximated A-optimality criterion via Neumann series and FGFT, improving efficiency and robustness in graph signal sampling.

Findings

Sampling strategy outperforms existing methods in experiments.

Reconstruction method is more noise-robust than least squares.

Efficient node selection without large matrix inversions.

Abstract

Graph sampling addresses the problem of selecting a node subset in a graph to collect samples, so that a K-bandlimited signal can be reconstructed in high fidelity. Assuming an independent and identically distributed (i.i.d.) noise model, minimizing the expected mean square error (MMSE) leads to the known A-optimality criterion for graph sampling, which is expensive to compute and difficult to optimize. In this paper, we propose an augmented objective based on Neumann series that well approximates the original criterion and is amenable to greedy optimization. Specifically, we show that a shifted A-optimal criterion can be equivalently written as a function of an ideal low-pass (LP) graph filter, which in turn can be approximated efficiently via fast graph Fourier transform (FGFT). Minimizing the new objective, we select nodes greedily without large matrix inversions using a matrix inverse lemma. Further, for the dynamic network case where node availability varies across time, we propose an extended sampling strategy that replaces offline samples one-by-one in the selected set. For signal reconstruction, we propose an accompanied biased signal recovery strategy that reuses the approximated filter from sampling. Experiments show that our reconstruction is more robust to large noise than the least square (LS) solution, and our sampling strategy far outperforms several existing schemes.

Tables

Table 1. TABLE I: Complexity Comparison of Different Sampling Strategies

	Preparation	Selection
Spectral Proxies	NONE	$𝒪\left(k\left|ℰ\right|M{T}_2\left(k\right)+NM\right)$
E-optimal	$𝒪\left(\left(\left|ℰ\right|M+R{M}^3\right)T_1\right)$	$𝒪\left(N{M}^4\right)$
MFN	$𝒪\left(\left(\left|ℰ\right|M+R{M}^3\right)T_1\right)$	$𝒪\left(N{M}^4\right)$
MIA	$𝒪\left(qN\left|ℰ\right|\right)$	$𝒪\left(NL{M}^{3.373}\right)$
Proposed GFS	$𝒪\left(N^2{\text{log}}^{2}N\right)$	$𝒪\left(N{M}^3\right)$

Table 2. TABLE II: Reconstruction MSE of The Proposed GFS Algorithm in Different Shift μ 𝜇 \mu at 0dB

Graph	$\mu$	Sample size
		100	110	120	130	140
G1	${10}^{-5}$	16.10	14.55	13.43	12.44	11.63
	$1/99$	16.07	14.59	13.43	12.46	11.64
	(47)	16.10	14.62	13.47	12.50	11.64
G2	${10}^{-5}$	20.77	18.68	17.09	15.77	14.63
	$1/99$	20.77	18.73	17.12	15.78	14.64
	(47)	21.31	19.36	17.59	16.33	15.08
G3	${10}^{-5}$	23.09	20.79	18.99	17.48	16.19
	$1/99$	22.98	20.77	18.99	17.49	16.21
	(47)	23.11	20.87	19.04	17.55	16.32

Equations

$$\begin{split}{{\mathbf{C}_{ij}}=\left\{\begin{array}[]{l}1,~{}~{}~{}{j={\mathcal{S}(i)}};\\ 0,~{}~{}~{}\textrm{otherwise},\end{array}\right.}\end{split}$$

$$\begin{array}{r}{\mathbf{C}}^{\ast }=\underset{\mathbf{C}\in {\mathbb{F}}^{M\times N}}{\mathrm{argmin}}\text{tr}\left([{(\mathbf{C}{\mathbf{V}}_K)}^{\top }\mathbf{C}{\mathbf{V}}_K{]}^{-1}\right).\end{array}$$

$$\begin{array}{r}0\le {\delta }_i<1\end{array}$$

$$\begin{array}{r}{\mathbf{C}}^{\ast }=\underset{\mathbf{C}\in {\mathbb{F}}^{M\times N}}{\mathrm{argmin}}\text{tr}\left([{(\mathbf{C}{\mathbf{V}}_K)}^{\top }\mathbf{C}{\mathbf{V}}_K+\mu \mathbf{I}{]}^{-1}\right)\end{array}$$

$$0\le {ϵ}_i\le 1$$

$${\mathbf{x}}^{\top }{\mathbf{T}}_{\mathcal{S}}\mathbf{x}={\left({\mathbf{C}}^{\top }\mathbf{x}\right)}^{\top }\left({\mathbf{V}}_K{\mathbf{V}}_K^{\top }\right)\left({\mathbf{C}}^{\top }\mathbf{x}\right)={\mathbf{y}}^{\top }\left({\mathbf{V}}_K{\mathbf{V}}_K^{\top }\right)\mathbf{y}$$

$$\begin{array}{r}0\le {\mathbf{x}}^{\top }{\mathbf{T}}_{\mathcal{S}}\mathbf{x}\le 1\end{array}$$

$$\begin{array}{r}{\mathcal{S}}^{\ast }={\mathrm{argmin}}_{\mathcal{S}:|\mathcal{S}|=M}\text{tr}{\left({\mathbf{T}}_{\mathcal{S}}+\mu \mathbf{I}\right)}^{-1},\end{array}$$

$$\begin{array}{r}[{(\mathbf{C}{\mathbf{V}}_K)}^{\top }\mathbf{C}{\mathbf{V}}_K+\mu \mathbf{I}{]}^{-1}=\sum _{l=0}^{\infty }(\mathbf{I}-\tilde{\mathbf{\Psi }}{)}^l\end{array}$$

$$\text{tr}(\mathbf{I}-\tilde{\mathbf{\Psi }}{)}^l=\text{tr}{\left[(1-\mu )\mathbf{I}-\mathbf{\Psi }\right]}^l$$

$$=\text{tr}\left[(1-\mu {)}^l\mathbf{I}+\sum _{d=1}^l(\genfrac{}{}{0ex}{}{l}{d})(1-\mu {)}^{(l-d)}(-\mathbf{\Psi }{)}^d\right]$$

$$\stackrel{(a)}{=}(1-\mu {)}^{l}K+\sum _{d=1}^l(\genfrac{}{}{0ex}{}{l}{d})(1-\mu {)}^{(l-d)}\text{tr}\left[(-{\mathbf{T}}_{\mathcal{S}}{)}^d\right]$$

$$=(1-\mu {)}^l(K-M)+\sum _{d=0}^l(\genfrac{}{}{0ex}{}{l}{d})(1-\mu {)}^{(l-d)}\text{tr}\left[(-{\mathbf{T}}_{\mathcal{S}}{)}^d\right]$$

$$=(1-\mu {)}^l(K-M)+\text{tr}{\left[(1-\mu )\mathbf{I}-{\mathbf{T}}_{\mathcal{S}}\right]}^l$$

$$\text{tr}\left[\sum _{l=0}^{\infty }(\mathbf{I}-\tilde{\mathbf{\Psi }}{)}^l\right]=\sum _{l=0}^{\infty }\text{tr}(\mathbf{I}-\tilde{\mathbf{\Psi }}{)}^l$$

$$=\sum _{l=0}^{\infty }(1-\mu {)}^l(K-M)+\text{tr}{\left[(1-\mu )\mathbf{I}-{\mathbf{T}}_{\mathcal{S}}\right]}^l$$

$$=(K-M)\frac{1-(1-\mu {)}^{\infty }}{\mu }+\text{tr}\sum _{l=0}^{\infty }{\left[(1-\mu )\mathbf{I}-{\mathbf{T}}_{\mathcal{S}}\right]}^l$$

$$\stackrel{(a)}{=}\frac{K-M}{\mu }+\text{tr}\sum _{l=0}^{\infty }{\left[\mathbf{I}-({\mathbf{T}}_{\mathcal{S}}+\mu \mathbf{I})\right]}^l$$

$$\begin{array}{r}\sum _{l=0}^{\infty }{\left[\mathbf{I}-({\mathbf{T}}_{\mathcal{S}}+\mu \mathbf{I})\right]}^l=({\mathbf{T}}_{\mathcal{S}}+\mu \mathbf{I}{)}^{-1}\end{array}$$

$$\text{tr}\left([{(\mathbf{C}{\mathbf{V}}_K)}^{\top }\mathbf{C}{\mathbf{V}}_K+\mu \mathbf{I}{]}^{-1}\right)=\frac{K-M}{\mu }+\text{tr}{\left({\mathbf{T}}_{\mathcal{S}}+\mu \mathbf{I}\right)}^{-1}$$

$${\min }_{\hat{\mathbf{\Lambda }},{\mathbf{S}}_1,...,{\mathbf{S}}_J}\Vert \mathbf{L}-{\mathbf{S}}_1...{\mathbf{S}}_J\hat{\mathbf{\Lambda }}{\mathbf{S}}_J^{\top }...{\mathbf{S}}_1^{\top }{\Vert }_F^2$$

$$\begin{array}{r}{\mathcal{S}}^{\ast }={\mathrm{argmin}}_{\mathcal{S}:|\mathcal{S}|=M}\text{tr}{\left({\mathbf{T}}_{\mathcal{S}}^{\text{FGFT}}+\mu \mathbf{I}\right)}^{-1}\end{array}$$

$$\kappa ({\mathbf{G}}_{\mathcal{S}})=\frac{{\lambda }_{\text{max}}({\mathbf{G}}_{\mathcal{S}})}{{\lambda }_{\text{min}}({\mathbf{G}}_{\mathcal{S}})}\le \frac{1+\mu }{\mu }\le {\kappa }_0$$

$$\displaystyle\hskip 8.5359pt\mathbf{M}^{-1}=\left[\begin{array}[]{l}\mathbf{A}~{}~{}~{}{\mathbf{U}}\\ {\mathbf{V}}~{}~{}~{}{\mathbf{C}}\\ \end{array}\right]^{-1}$$

$$\displaystyle=\left[\begin{array}[]{l}\mathbf{A}^{-1}+\mathbf{A^{-1}UH^{-1}VA^{-1}}~{}~{}~{}-{\mathbf{A^{-1}UH^{-1}}}\\ -{\mathbf{H^{-1}VA^{-1}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\kern 1.0pt}{\mathbf{H^{-1}}}\\ \end{array}\right]$$

$$\displaystyle\mathbf{G}_{\mathcal{S}\cup\{i\}}=\left[\begin{array}[]{l}\mathbf{G}_{\mathcal{S}}~{}~{}~{}~{}~{}~{}~{}{\mathbf{G}}_{\mathcal{S},\{i\}}\\ {\mathbf{G}}_{\{i\},\mathcal{S}}~{}~{}~{}{\mathbf{G}}_{ii}\\ \end{array}\right]\doteq\left[\begin{array}[]{l}\mathbf{G}_{\mathcal{S}}~{}~{}~{}\mathbf{g}_{i}\\ ~{}\mathbf{g}^{\top}_{i}~{}~{}\mathbf{G}_{ii}\\ \end{array}\right]$$

$$\displaystyle\mathbf{G}^{-1}_{\mathcal{S}\cup\{i\}}=\left[\begin{array}[]{l}\mathbf{G}^{-1}_{\mathcal{S}}+h^{-1}\mathbf{G}^{-1}_{\mathcal{S}}\mathbf{g}_{i}\mathbf{g}^{\top}_{i}\mathbf{G}^{-1}_{\mathcal{S}}~{}~{}~{}~{}-h^{-1}\mathbf{G}^{-1}_{\mathcal{S}}\mathbf{g}_{i}\\ ~{}~{}~{}~{}~{}~{}~{}~{}~{}-h^{-1}\mathbf{g}^{\top}_{i}\mathbf{G}^{-1}_{\mathcal{S}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}h^{-1}\\ \end{array}\right]$$

$$\alpha \ge \frac{{\mu }^3(\mu +2)}{(\mu +1{)}^4}$$

$$\frac{g(\mathcal{S})-g({\mathcal{S}}^{\ast })}{g(\{\})-g({\mathcal{S}}^{\ast })}\le e^{-\alpha }$$

$$\begin{array}{r}{\left(\mathbf{A}+\mathbf{u}{\mathbf{v}}^{\top }\right)}^{-1}={\mathbf{A}}^{-1}-\frac{{\mathbf{A}}^{-1}\mathbf{u}{\mathbf{v}}^{\top }{\mathbf{A}}^{-1}}{1+{\mathbf{v}}^{\top }{\mathbf{A}}^{-1}\mathbf{u}}\end{array}$$

$$\displaystyle\mathbf{G}_{\mathcal{S}}=\left[\begin{array}[]{l}{\mathbf{G}}_{\mathcal{S}\setminus\{j\}}~{}~{}~{}~{}~{}~{}~{}{\mathbf{G}}_{\mathcal{S}\setminus\{j\},\{j\}}\\ {\mathbf{G}}_{\{j\},\mathcal{S}\setminus\{j\}}~{}~{}~{}{\mathbf{G}}_{jj}\\ \end{array}\right]$$

$$\displaystyle\mathbf{G}_{\mathcal{\tilde{S}}}=\left[\begin{array}[]{l}{\mathbf{G}}_{\mathcal{S}\setminus\{j\}}~{}~{}~{}~{}~{}~{}~{}{\mathbf{G}}_{\mathcal{S}\setminus\{j\},\{k\}}\\ {\mathbf{G}}_{\{k\},\mathcal{S}\setminus\{j\}}~{}~{}~{}{\mathbf{G}}_{kk}\\ \end{array}\right]$$

$${\mathbf{G}}_{\tilde{\mathcal{S}}}={\mathbf{G}}_{\mathcal{S}}+{\mathbf{e}}_i{\mathbf{p}}^{\top }+\mathbf{q}{\mathbf{e}}_i^{\top }$$

$${\mathbf{F}}^{-1}={\mathbf{G}}_{\mathcal{S}}^{-1}-\frac{{\mathbf{G}}_{\mathcal{S}}^{-1}{\mathbf{e}}_i{\mathbf{p}}^{\top }{\mathbf{G}}_{\mathcal{S}}^{-1}}{1+{\mathbf{p}}^{\top }{\mathbf{G}}_{\mathcal{S}}^{-1}{\mathbf{e}}_i}$$

Full text

I Introduction

Graph signal processing (GSP) is the study of signals that reside on irregular data kernels described by graphs [1, 2, 3]. To analyze graph signals spectrally, a large amount of research strived to define frequencies on graphs, and subsequently design transforms and wavelets for spectral decomposition of signals [4, 5, 6]. Specifically, graph Fourier transform (GFT) of a graph signal $\mathbf{x}\in\mathbb{R}^{N}$ is defined by its projection on the eigenvector space of the graph Laplacian matrix [2] (or the adjacency matrix [7]) 111[8] proposed an alternative definition of generalized GFT that takes the irregularity of the graph into account. For direct graphs, [9] proposed the notion of graph direct variation. In this paper, we focus on GFT derived from the Laplacian operator, but our sampling method is applicable to other symmetrical graph operators.. A graph signal of bandwidth $K$ is a signal with non-zeros GFT coefficients for $K$ eigenvectors associated with the $K$ smallest eigenvalues. Sampling of bandlimited graph signals is an important basic problem in GSP, since sensing in practice is often expensive, e.g., SAR images from remote sensing [10] and label assignment of a medical examination [11]. Formally, graph sampling is the study of node selection in a graph to collect samples, such that the original signal can be reconstructed in high fidelity.

Previous works in graph sampling can be broadly divided into two categories: noiseless sampling and noisy sampling. Given that samples are noiselessly observed, [12] showed that random sampling can result in perfect signal reconstruction with high probability, assuming that the sample size is larger than the signal’s bandwidth. However, the condition number of the corresponding coefficient matrix can be large, resulting in an unstable signal reconstruction. In response, [13, 14] proved a necessary and sufficient condition of a uniqueness set [15], and then proposed a sampling strategy with simple matrix-vector computations based on the notion of spectral proxies (SP).

If the observed samples are corrupted with noise, one classical criterion is to optimize a minimum mean square error (MMSE) function, which leads to an A-optimality design criterion assuming an independent and identically distributed (i.i.d.) additive noise model [16]. Minimizing the A-optimality criterion for graph sampling is known to be difficult: evaluating the criterion for a fixed sample set requires computation of eigenvectors and matrix inverse, and optimizing the sample set for a fixed budget is combinatorial.

To avoid heavy computation, [17] used a greedy procedure to optimize the MMSE problem directly, called minimum Frobenius norm (MFN). Further, [18] proved that the performance of the greedy solution can be bounded via super-modularity analysis. Alternatively, [19] used the E-optimality criterion for graph sampling, which can be interpreted as the worst case MSE. Nonetheless, all proposals above required eigen-decomposition of the Laplacian operator and costly computations in each greedy step, such as matrix inverse for MFN and singular value decomposition (SVD) for E-optimal sampling.

To reduce the complexity of evaluating the objective, an eigen-decomposition-free sampling strategy was recently proposed in [20, 21] by constraining the coverage of a localization operator, implemented using Chebyshev polynomial approximation 222Our proposed method is a function of an ideal graph low-pass filter, which can be viewed as an ideal operator-based sampling defined in [21]. However, we do not implement this filter by a localized operator. . However, the performance of this method is competitive only when the sample size is smaller than signal’s bandwidth. Orthogonally, our previous work, called matrix inversion approximation (MIA) [22], developed a sampling strategy to optimize an approximate MSE function based on truncated Neumann series, but the number of terms in the truncated sum must be sufficiently large for the approximation error to be small, limiting its practicality.

Because noisy graph sampling problem is NP-hard in nature [18], all aforementioned criteria are optimized in a greedy manner. To lower complexity, random sample node selection was proposed to sample signals based on a specified probability [23, 24]. However, for the same sampling budget it cannot guarantee stable performance comparable to the deterministic counterparts [1, 21].

In this paper, we propose a novel low-complexity deterministic sampling strategy via Neumann series expansion. Unlike previous approaches [17, 19, 21], we minimize a variant of the A-optimality criterion directly but without eigen-decomposition or matrix inversion. Yet different from our work in [22], there is no explicit computation of the Neumann series sum in our augmented objective, and thus difficult truncation / approximation tradeoff is not necessary.

Specifically, we first propose an augmented objective to approximate the A-optimality criterion, and prove the inverse of its information matrix is equal to a matrix series based on the Neumann series theorem [25]. We then rewrite the objective as a function of an ideal low-pass (LP) graph filter with cutoff frequency $K$, which is efficiently approximated using fast Graph Fourier Transform (FGFT) [26]. To optimize the objective function, we select nodes greedily without any matrix inversions based on a matrix inverse lemma. Further, the greedy solution is proved to have upper-bounded performance based on super-modularity analysis [18].

All previous graph sampling works assumed a static node set, while in practice often the availability of nodes as sensors varies over time; e.g., one person goes offline in a social network, or one sensor runs out of battery in a sensor network [27]. We additionally address the sampling problem for dynamic time-varying node subsets, where the availability of each node varies as a function of time. Specifically, we develop an extended sampling strategy that replaces offline samples one-by-one in the selected set via Sherman-Morrison formula [28] to reduce complexity. *To the best of our knowledge, we are the first to formulate and tackle the sampling problem on dynamic node subsets in the literature. *

Finally, we develop an accompanied signal reconstruction strategy to reuse the LP graph filter from sampling. This reconstruction strategy provides a biased estimator with lower complexity, which is more robust to large noise than the conventional least square (LS) solution [14].

To summarize, our main contributions are as follows:

We propose a novel graph sampling objective based on Neumann series expansion, which well approximates the A-optimality criterion. 2. 2.

We design a fast sampling algorithm to minimize the proposed objective with neither eigen-decomposition nor matrix inversions. 3. 3.

For the dynamic network case, we propose an extended node replacement strategy to replace offline nodes one-by-one efficiently. 4. 4.

We propose an accompanied biased signal recovery method for sampling that is more robust to large noise. 5. 5.

Experiments demonstrate that our proposed sampling method outperforms several state-of-the-art sampling strategies in artificial and real-world graphs.

The outline of this paper is as follows. In Section II, we overview GSP concepts and terminology used throughout this paper and discuss the A-optimality sampling problem. We propose the augmented A-optimal sampling criterion in Section III. In Section IV, we detail our fast graph sampling algorithm. We discuss the dynamic graph sampling problem in Section V. We develop our biased signal reconstruction method in Section VI. Finally, we present experimental results nad conclusion in Section VII and VIII, respectively.

II Preliminaries

Denote by $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{W})$ a graph with $N$ nodes indexed by $\mathcal{V}=\{1,\ldots,N\}$. $\mathcal{E}$ specifies the set of connected node pairs $(i,j)$, and the $(i,j)$-th entry $w_{i,j}$ of an $N\times N$ adjacency matrix $\mathbf{W}$ is the weight of an edge connecting nodes $i$ and $j$ ($w_{i,j}=0$ if nodes $i$ and $j$ are not connected). We additionally define a diagonal degree matrix $\mathbf{D}$, where $d_{i,i}=\sum_{j}w_{i,j}$. In this paper we focus on connected, undirected graphs with no self-loops, and we adopt the symmetric combinatorial graph Laplacian matrix $\mathbf{L}=\mathbf{D}-\mathbf{W}$ as the variation operator. Suppose that the eigen-decomposition of $\mathbf{L}$ is $\mathbf{L}=\mathbf{V}\boldsymbol{\Sigma}{{\mathbf{V}^{\top}}}$, where $\mathbf{V}={[\mathbf{v}_{1}},...,{\mathbf{v}_{N}}]$ contains a set of $N$ orthonormal eigenvectors as columns corresponding to non-decreasing eigenvalues $0={\lambda_{1}}<{\lambda_{2}}\leq{\lambda_{3}}\leq\ldots\leq{\lambda_{N}}$. Then the GFT of a graph signal $\mathbf{x}\in\mathbb{R}^{N}$ is defined as its expansion on the eigenvector space of $\mathbf{L}$, i.e., $\tilde{\mathbf{x}}={\mathbf{V}^{\top}}\mathbf{x}$, and the inverse GFT is $\mathbf{x}=\mathbf{V}\tilde{\mathbf{x}}$. A graph signal is called $K$-bandlimited if its GFT coefficients $\tilde{\mathbf{x}}$ are non-zero only for the lowest $K$ frequencies 333Authors in [18] defined the $\mathcal{K}$-spectrally sparse graph signal by ${\mathbf{x}={\mathbf{V}_{\mathcal{K}}}{{\tilde{\mathbf{x}}}_{\mathcal{K}}}}$ which is a generalization of our formulation since $\mathcal{K}\in\mathcal{V}$ is not restricted to be the first $K$ elements. We use the more conventional definition here, but our proposed algorithm is applicable to the general case.. A $K$-bandlimited graph signal can be expressed as ${\mathbf{x}={\mathbf{V}_{K}}{{\tilde{\mathbf{x}}}_{K}}}$, where ${\mathbf{V}_{K}}$ means the first ${K}$ columns of ${\mathbf{V}}$, and ${{{\tilde{\mathbf{x}}}_{K}}}$ denotes the first ${K}$ elements of ${\tilde{\mathbf{x}}}$. To describe sampling on bandlimited graph signals, we define a sampling operator at first [29].

Definition 1**.**

To sample $M$ elements from $\mathbf{x}$ to produce ${\mathbf{x}_{\mathcal{S}}=\mathbf{Cx}\in{{\rm{\mathbb{R}}}^{M}}}$ with ${\left|\mathcal{S}\right|=M}$ and ${\mathcal{S}\subseteq\mathcal{V}}$, we define a 0-1 binary sampling matrix $\mathbf{C}$ as

[TABLE]

where $\mathcal{S}$ is the set of sampling indices, $\mathcal{S}(i)$ means the $i$-th element of set $\mathcal{S}$, and ${\left|\mathcal{S}\right|}$ is the number of elements in ${\mathcal{S}}$.

A sampled $K$-bandlimited graph signal can now be written as ${{\mathbf{x}_{\mathcal{S}}}=\mathbf{C}{\mathbf{V}_{K}}{{\tilde{\mathbf{x}}}_{K}}}$. Sampling operators satisfying $\textrm{rank}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$ are called qualified sampling operators in [19], since any $K$-bandlimited graph signal $\mathbf{x}$ can be perfectly recovered from samples $\mathbf{x}_{\mathcal{S}}$ on those operators. Specifically, if $\textrm{rank}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$, then, from a sampled graph signal, the LS solution can provide a perfect and unique reconstruction in noiseless environment, i.e., $\hat{\mathbf{x}}_{\text{LS}}={\mathbf{V}_{K}}{\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{\dagger}}\mathbf{x}_{\mathcal{S}}$, where ${\dagger}$ denotes the pseudo-inverse computation [14]. Further, $\textrm{rank}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$ was demonstrated empirically to be fulfilled with high probability via random node selection when $M\geq K$ [12]. In the sequel, we focus on the sampling region $M\geq K$ and assume $\textrm{rank}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$ is satisfied.

When the samples $\mathbf{x}_{\mathcal{S}}$ are corrupted by additive noise $\mathbf{n}_{\mathcal{S}}$, the LS recovery will produce a minimum variance unbiased estimator $\hat{\mathbf{x}}_{\text{LS}}={\mathbf{V}_{K}}{\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{\dagger}}\mathbf{y}_{\mathcal{S}}={\mathbf{V}_{K}}{\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{{\dagger}}}(\mathbf{x}_{\mathcal{S}}+\mathbf{n}_{\mathcal{S}})=\mathbf{x}+{\mathbf{V}_{K}}{\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{\dagger}}\mathbf{n}_{\mathcal{S}}$ [30]. Assuming that noise $\mathbf{n}_{\mathcal{S}}$ is i.i.d. with zero mean and unit variance, the covariance matrix of the reconstruction error is $\mathbf{R}_{\hat{\mathbf{x}}_{\text{LS}}}=\mathbb{E}\left[(\hat{\mathbf{x}}_{\text{LS}}-\mathbb{E}[\hat{\mathbf{x}}_{\text{LS}}])(\hat{\mathbf{x}}_{\text{LS}}-\mathbb{E}[\hat{\mathbf{x}}_{\text{LS}}])^{\top}\right]={\mathbf{V}_{K}}\left[({\mathbf{C}{\mathbf{V}_{K}})^{\top}}\mathbf{C}{\mathbf{V}_{K}}\right]^{-1}{\mathbf{V}^{\top}_{K}}$. By the theory of optimal experiments design [31], finding a sampling matrix $\mathbf{C}$ to minimize the trace of the covariance matrix leads to the known A-optimality criterion:

[TABLE]

where $\mathbb{F}^{M\times N}$ the set of all $\mathbf{C}$ defined in (1).

Minimizing the A-optimality objective (2) directly using $\mathbf{C}$ is difficult because finding an optimal $\mathbf{C}^{*}$ is NP-hard and evaluating the A-optimality value for a given $\mathbf{C}$ requires expensive computations of $\mathbf{V}_{K}$ and matrix inverse.

For simplicity, in the sequel the complement set of $\mathcal{S}$ is denoted by $\mathcal{S}^{c}$. $|\mathcal{E}|$ means the cardinality of set $\mathcal{E}$, and ${{\mathbf{A}_{{\mathcal{S}_{1}}{\mathcal{S}_{2}}}}}$ is the sub-matrix of a matrix ${\mathbf{A}}$ with rows and columns indexed by ${\mathcal{S}_{1}}$ and ${\mathcal{S}_{2}}$, respectively. ${{\mathbf{A}_{{\mathcal{S}}{\mathcal{S}}}}}$ is simplified to ${{\mathbf{A}_{{\mathcal{S}}}}}$. $\mathbf{I}$ is the identity matrix whose dimension depends on the context.

III Augmented A-optimality Graph Signal Sampling

In this section, we will propose a new sampling objective to approximate the A-optimality criterion. We first review the Neumann series theorem in [25] and one result in [22], both of which will be used for proving a forthcoming Theorem.

Theorem 1**.**

(Neumann series theorem) If the eigenvalues $\lambda_{i}$ of a square matrix $\mathbf{A}$ have the property that $\rho(\mathbf{A})\doteq\mathop{\max}\limits_{i}\left|{{\lambda_{i}}}\right|<1$, then its Neumann series $\mathbf{I}+\mathbf{A}+{\mathbf{A}^{2}}+\cdots$ converges to ${(\mathbf{I}-\mathbf{A})^{-1}}$, i.e., ${(\mathbf{I}-\mathbf{A})^{-1}}=\sum\limits_{l=0}^{\infty}{\mathbf{A}^{l}}$.

Proposition 1**.**

Denote by $\mathbf{\Psi}={\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{\top}}\mathbf{C}{\mathbf{V}_{K}}$, $\mathbf{\Phi}=\mathbf{I}-\mathbf{\Psi}$ and ${\delta_{1}}\leq\ldots\leq{\delta_{K}}$ as the eigenvalues of $\mathbf{\Phi}$. Then,

[TABLE]

if ${{\mathbf{C}{\mathbf{V}_{K}}}}$ is full column rank, i.e., $\emph{{rank}}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$.

Proof.

The proof is in Appendix A. ∎

These two results will support us to prove the Theorem 2 in the next subsection.

III-A Augmented A-optimality Criterion

We propose an augmented objective by adding a weighted identity matrix to (2):

[TABLE]

where $\mu$ is a small weight (shift) parameter with $0<\mu<1$. We discuss its design in details later.

Problem (4) is also difficult to optimize directly like equation (2). Before proposing its alternative simpler form, we first present the next Proposition which will be employed to prove the following Theorem.

Proposition 2**.**

Denote by ${\epsilon_{1}}\leq\ldots\leq{\epsilon_{M}}$ the eigenvalues of ${\mathbf{T}_{\mathcal{S}}}$, then

[TABLE]

where ${\mathbf{T}}={\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}$ is an ideal LP graph filter with cutoff frequency $K$.

Proof.

$\forall\mathbf{x}\in{\mathbb{R}^{M}}$ and $\left\|\mathbf{x}\right\|_{2}=1$,

[TABLE]

where $\mathbf{y}=\mathbf{C}^{\top}{\mathbf{x}}$. Hence $\mathbf{y}_{\mathcal{S}}={\mathbf{x}}$, $\mathbf{y}_{\mathcal{S}^{c}}={\mathbf{0}}$ and $\left\|\mathbf{y}\right\|_{2}=1$.

Since $\mathbf{T}$ is symmetric and its eigen-decomposition can be written as $\mathbf{T}={{\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}}={\mathbf{V}}\text{diag}\{1,...,1,0,...,0\}\mathbf{V}^{\top}$, its eigenvalues are 0 or 1. Then, due to the Rayleigh quotient theorem [25], $\forall\mathbf{y}\in{\mathbb{R}^{N}}$ and $\left\|\mathbf{y}\right\|_{2}=1$, $0\leq{\mathbf{y}^{\top}}({{\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}})\mathbf{y}\leq 1$, which holds also for some specific $\mathbf{y}$. Thus,

[TABLE]

which implies Proposition 2 via the Rayleigh quotient theorem. ∎

Now, we propose and prove the next Theorem by employing the results from Theorem 1, Proposition 1 and Proposition 2.

Theorem 2**.**

The augmented A-optimal objective (4) is equivalent to

[TABLE]

if the selected matrix ${{\mathbf{C}{\mathbf{V}_{K}}}}$ is full column rank.

Proof.

Let $\tilde{\mathbf{\Psi}}={{{(\mathbf{C}{\mathbf{V}_{K}})}^{\top}}\mathbf{C}{\mathbf{V}_{K}}}+\mu\mathbf{I}=\mathbf{\Psi}+\mu\mathbf{I}$, and $\tilde{\mathbf{\Phi}}=\mathbf{I}-\tilde{\mathbf{\Psi}}=\mathbf{\Phi}-\mu\mathbf{I}$. We know from Proposition 1 that $0\leq{\delta_{i}}<1$. Hence the eigenvalues of $\tilde{\mathbf{\Phi}}$ are in [$-\mu$,$1-\mu$). Since the shift parameter $0<\mu<1$, $\rho(\tilde{\mathbf{\Phi}})<1$, and by applying Theorem 1

[TABLE]

Leveraging on a property of the trace operator, we write

[TABLE]

where $(a)$ holds because $\textrm{tr}(\mathbf{\Psi}^{d})=\textrm{tr}({\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}}\mathbf{C}{\mathbf{V}_{K}}...{\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}}\\ \mathbf{C}{\mathbf{V}_{K}})=\textrm{tr}(\mathbf{C}{\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}}...\mathbf{C}{\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}})=\textrm{tr}(\mathbf{T}^{d}_{\mathcal{S}})$.

Hence,

[TABLE]

where equality $(a)$ holds since $0<\mu<1$.

From Proposition 2, the eigenvalues of $\mathbf{I}-({\mathbf{T}_{\mathcal{S}}}+\mu\mathbf{I})$ are in [$-\mu$,$1-\mu$]. Again, since the shift parameter $\mu$ is small and positive, $\rho\left[\mathbf{I}-({\mathbf{T}_{\mathcal{S}}}+\mu\mathbf{I})\right]\leq 1-\mu<1$. Using Theorem 1, we have

[TABLE]

Combining (8), (III-A) and (11),

[TABLE]

which implies Theorem 2 because $M$, $K$ and $\mu$ are all constant in a sampling problem. ∎

Optimizing the new problem (7) requires computation of an ideal LP filter rather than first $K$ eigenvectors of $\mathbf{L}$. We discuss an efficient approximation of the LP filter next.

III-B Ideal Low-pass Graph Filter Approximation

One classical method is to approximate an ideal LP filter $\mathbf{T}$ in (7) by using a Chebyshev matrix polynomial of $\mathbf{L}$, i.e., $\mathbf{T}^{\textrm{Poly}}=\sum_{i=1}^{n}\left(\sum_{j=0}^{q}\beta_{j}\lambda^{j}_{i}\right)\mathbf{v}_{i}\mathbf{v}^{\top}_{i}=\sum_{j=0}^{q}\beta_{j}\mathbf{L}^{j}$ [5]. Another method is to apply the Jacobi eigenvalue algorithm: approximately diagonalize $\mathbf{L}$ with an estimated eigenvector matrix $\tilde{\mathbf{V}}=\mathbf{S}_{1}...\mathbf{S}_{J}$ [26]. Specifically, the goal is to solve the following optimization problem:

[TABLE]

where $\hat{\boldsymbol{\Lambda}}$ is constrained to be a diagonal matrix and $\mathbf{S}_{i}$ are constrained to be Givens rotation matrices.

A truncated Jacobi algorithm [32] was adopted to optimize (12) iteratively. With the optimized $\tilde{\mathbf{V}}$ where $\mathbf{L}=\tilde{\mathbf{V}}\hat{\boldsymbol{\Lambda}}\tilde{\mathbf{V}}^{\top}$, an ideal LP filte $\mathbf{T}$ can be implemented as $\mathbf{T}^{\textrm{FGFT}}=\tilde{\mathbf{V}}_{K}\tilde{\mathbf{V}}^{\top}_{K}$.

We adopt $\mathbf{T}^{\textrm{FGFT}}$ for the following reasons:

1. (1)

[26] claimed the superiority of $\mathbf{T}^{\textrm{FGFT}}$ compared with $\mathbf{T}^{\textrm{Poly}}$ when approximating ideal graph LP filters; specifically, $\frac{\left\|\mathbf{T}-\mathbf{T}^{\textrm{FGFT}}\right\|_{F}}{\left\|\mathbf{T}\right\|_{F}}<\frac{\left\|\mathbf{T}-\mathbf{T}^{\textrm{Poly}}\right\|_{F}}{\left\|\mathbf{T}\right\|_{F}}$ under the same approximation complexity [26]. In our work, $\mathbf{T}$ is an ideal graph LP filter. Moreover, if the bandwidth of graph signals $K$ varies in one graph, we can reuse the same $\tilde{\mathbf{V}}$ when realizing $\mathbf{T}^{\textrm{FGFT}}$. 2. (2)

A sufficient condition of Proposition 2 is that the eigenvalues of $\mathbf{T}$ are in $\{0,1\}$, based on which we can claim the eigenvalues of $\mathbf{I}-({\mathbf{T}_{\mathcal{S}}}+\mu\mathbf{I})$ are in [$-\mu$,$1-\mu$] in Theorem 2. Then, (11) can hold with a small positive parameter $\mu$. Small $\mu$ is also required for (4) to approximate the original A-optimality criterion (2) well. Because $\tilde{\mathbf{V}}$ estimated from Givens matrices is orthogonal, Proposition 2 also holds for $\mathbf{T}^{\textrm{FGFT}}_{\mathcal{S}}$. In contrast, the eigenvalue scope of $\mathbf{T}^{\textrm{Poly}}_{\mathcal{S}}$ is unpredictable.

We now write our final augmented objective as

[TABLE]

which requires just an approximated LP filter operator rather than eigen-decomposition of the graph Laplacian operator. For simplicity, we write $\mathbf{G}=\mathbf{T}^{\textrm{FGFT}}+\mu\mathbf{I}$ in the sequel.

III-C Selection of Shift Parameter $\mu$

We now discuss how to select an appropriate shift $\mu$. To well approximate the original criterion (2), the shift $\mu$ ought to be as small as possible, but a small $\mu$ would cause the matrix inverse computation in (13) to be unstable. To ensure numerical stability, we can bound the condition number $\kappa(\mathbf{G})$ of $\mathbf{G}$, where $\kappa(\mathbf{G})=\lambda_{\text{max}}(\mathbf{G})/\lambda_{\text{min}}(\mathbf{G})$, using $\mu$ as follows [33]. As discussed, Proposition 2 also holds for $\mathbf{T}^{\textrm{FGFT}}_{\mathcal{S}}$, so $\mu\leq\lambda\left(\mathbf{G}_{\mathcal{S}}\right)\leq 1+\mu$. Hence we can bound $\kappa(\mathbf{G}_{\mathcal{S}})$ as follows:

[TABLE]

where $\kappa_{0}$ is the upper bound of an acceptable condition number. From the latter part, $\mu\geq\frac{1}{\kappa_{0}-1}$. Because $\mu$ should be as small as possible, the optimal $\mu$ is $\mu^{*}=\frac{1}{\kappa_{0}-1}$. We will set $\kappa_{0}=100$ in our experiments and compare its performance with other design strategies via simulations.

IV Fast Greedy Sampling

It is difficult to compute an optimal solution of our proposed criterion (13) with reasonable complexity since the problem is combinatorial. Greedy algorithms have been commonly used in the graph sampling literature [14, 19]. However, if (13) is minimized naïvely using a greedy scheme, then the algorithm needs to perform one matrix inversion to evaluate each candidate sample set $\mathcal{S}$. To alleviate this computation burden, we propose an accelerated greedy algorithm to avoid the matrix inverse operation.

IV-A Accelerated Greedy Sampling

Our accelerated algorithm is based on the block matrix inversion formula, reviewed as follows [34].

Lemma 1**.**

The inverse of matrix $\mathbf{M}$ can be computed using the inverse of sub-matrix $\mathbf{A}$ and the inverse of its Schur complement, i.e.,

[TABLE]

where $\mathbf{H}=\mathbf{C-VA^{-1}U}$ is the Schur complement $\mathbf{M}/\mathbf{A}$ of sub-matrix $\mathbf{A}$ of matrix $\mathbf{M}$.

We observe that matrix $\mathbf{G}$ is symmetric, and under some permutation, its sub-matrix $\mathbf{G}_{\mathcal{S}\cup\{i\}}$ can be expressed as

[TABLE]

where $\mathbf{g}_{i}$ denotes the partial vector of $i$-th column of $\mathbf{G}$ indexed by $\mathcal{S}$.

From Lemma 1, we can compute the inverse of $\mathbf{G}_{\mathcal{S}\cup\{i\}}$ as

[TABLE]

where $h=\mathbf{G}_{ii}-\mathbf{g}^{\top}_{i}\mathbf{G}^{-1}_{\mathcal{S}}\mathbf{g}_{i}$ is actually a scalar here.

$\mathbf{G}^{-1}_{\mathcal{S}\cup\{i\}}$ can thus be computed using $\mathbf{G}^{-1}_{\mathcal{S}}$ stored in the last iteration. Instead of computing $\mathbf{G}^{-1}_{\mathcal{S}\cup\{i\}}$ directly, (27) performs only two matrix-vector products and then stacks the results. We outline our algorithm in pseudo-code in Algorithm 1 and call it graph filter submatrix (GFS)-based sampling algorithm. We note that GFS requires neither computation of eigenvectors of $\mathbf{L}$ nor large matrix inversion, and it produces exactly the same sampling set as naïvely optimizing problem (13) with a greedy scheme.

We next study the performance bound of this greedy solution via super-modularity analysis.

IV-B Super-modularity Analysis

Let $f(\mathcal{S})=\textrm{tr}\left({\mathbf{T}^{\textrm{FGFT}}_{\mathcal{S}}}+\mu\mathbf{I}\right)^{-1}$ and $g(\mathcal{S})=\text{tr}[(\mathbf{C}\tilde{\mathbf{V}}_{K})^{\top}\mathbf{C}\tilde{\mathbf{V}}_{K}+\mu\mathbf{I}]^{-1}$. Because Proposition 2 holds for $\mathbf{T}^{\textrm{FGFT}}_{\mathcal{S}}$, following a similar proof to one in Theorem 2, we have $f(\mathcal{S})=g(\mathcal{S})+\frac{M-K}{\mu}$. Hence, the solution of greedily optimizing $f(\mathcal{S})$ is exactly the same as that of greedily optimizing $g(\mathcal{S})$, but with much lower complexity. Since $g(\mathcal{S})$ approximates the real MSE value, it is reasonable to evaluate the sub-optimality of its greedy solution.

Lemma 2**.**

The set function $g(\mathcal{S})=\emph{\text{tr}}[(\mathbf{C}\tilde{\mathbf{V}}_{K})^{\top}\mathbf{C}\tilde{\mathbf{V}}_{K}+\mu\mathbf{I}]^{-1}$ is (i) monotone decreasing and (ii) $\alpha$-supermodular with

[TABLE]

Proof.

The detailed proof is shown in Appendix B ∎

Let $\mathcal{S}^{*}$ be the optimal solution of minimizing $g(\mathcal{S})$ when $|\mathcal{S}|=M$, and $\mathcal{S}$ be the set obtained by applying the GFS algorithm. We have the following result using Lemma 2:

[TABLE]

which is easily derived from Theorem 2 in [18].

From (29), the value of $g(\mathcal{S})$ will be bounded by a function of $g(\mathcal{S^{*}})$, thus the sub-optimality of our GFS solution is up-bounded.

IV-C Complexity Analysis

We here analyze the complexity of the proposed GFS sampling strategy, assuming $M=K$. The overall complexity is divided into two parts: “Preparation” and “Selection”. The preparation step includes the complexity of computing the prior information, e.g., eigen-decomposition, and the selection step has the complexity of collecting samples. GFS needs to compute matrix $\mathbf{G}$ in the preparation step. Authors in [26] has shown that the complexity of computing $\tilde{\mathbf{V}}$ via parallel truncated Jacobi algorithm is $\mathcal{O}(NJ\text{log}N)$. During experiments, the number of Givens rotation matrices is set at $J=\mathcal{O}(N\text{log}N)$. Once we obtain $\mathbf{S}_{1},...,\mathbf{S}_{J}$, $\mathbf{T}^{\textrm{FGFT}}$ can be implemented via $\mathbf{T}^{\textrm{FGFT}}=\mathbf{S}_{1}...\mathbf{S}_{J}\mathbf{B}\mathbf{S}_{J}^{\top}...\mathbf{S}_{1}^{\top}$, where $\mathbf{B}$ is a diagonal matrix with ones on the first $K$ diagonal elements and zeros for the rest. Since every Givens matrix $\mathbf{S}_{i}$ is sparse with 4 non-zero entries, $\mathbf{T}^{\textrm{FGFT}}$ can be computed iteratively with complexity $\mathcal{O}(NJ)$. For the selection step, the complexity for matrix-vector product is $\mathcal{O}(M^{2})$. Considering $|\mathcal{S}|=M$ and $|\mathcal{S}^{c}|=N-M$, the complexity of GFS in the sampling step is at most $\mathcal{O}\left(NM^{3}\right)$ as shown in Table I.

Table I shows the complexity of different sampling algorithms. In the preparation step, the E-optimal and MFN sampling methods require $\mathbf{V}_{K}$, whose complexity is $\mathcal{O}\left({\left({\left|\mathcal{E}\right|M+R{M^{3}}}\right){T_{1}}}\right)$ via the block version of the Rayleigh quotient minimization method [35]. $T_{1}$ denotes the average number of iterations for convergence of this method, and $R$ is a constant for mixing the complexity of two iterative steps in this method [35]. MIA requires the $K$-th eigenvalue of $\mathbf{L}$ and an approximate LP filter $\mathbf{T}^{\text{Poly}}$, whose complexity is $\mathcal{O}(qN|\mathcal{E}|)$, where $q$ is the length of a Chebyshev matrix polynomial. In the selection step, the first eigen-pair of $\left((\mathbf{L}^{\top})^{k}\mathbf{L}^{k}\right)_{\mathcal{S}^{c}}$ is required for SP, and its complexity is $\mathcal{O}\left({k\left|\mathcal{E}\right|M{T_{2}}\left(k\right)}\right)$, where $k$ is the order of the proxies approximation, and $T_{2}(k)$ is the number of iterations for convergence of one single eigen-pair [14]. E-optimal requires SVD, and MFN requires matrix inversion computation for selection, so their complexity in this step are both $\mathcal{O}(NM^{4})$. The matrix multiplication in MIA has an asymptotic complexity of $\mathcal{O}(M^{2.373})$, thus its complexity in the selection step is $\mathcal{O}(NLM^{3.373})$, where $L$ is a truncation parameter [22]. From Table I, we can see that the complexity of the proposed GFS method is lower than the MIA, MFN and E-optimal method in the selection step. Numerical analysis for complexity will be illustrated in the experimental section.

V Dynamic Subset Sampling

Most existing graph sampling works assume a static graph from which to select nodes [17, 19, 14], although, as argued in the Introduction, the availability of sensing nodes varies over time for many practical applications. In this section, we extend our proposed algorithm to address sampling on dynamic graphs: how to select nodes for sampling if the available nodes are a time-varying subset? For concreteness, we first develop a mathematical model to generate a node subset $\mathcal{A}_{t}$ at time $t$ that slowly evolves over time.

At $t=0$, we first randomly generate $P_{0}N$ nodes to form an initial available node subset $\mathcal{A}_{0}$ from $\mathcal{V}$, where $P_{0}$ is a fixed probability. Denote by $z^{t}_{j}$ the state of node $j$ at time $t$, where $z^{t}_{j}=1$ if $j\in\mathcal{A}_{t}$, and $z^{t}_{j}=0$ otherwise, i.e. $j$ is in the complement set $\mathcal{A}^{c}_{t}$. Define the state crossover probability as $P(z^{t+1}_{j}=0|z^{t}_{j}=1)=P(z^{t+1}_{j}=1|z^{t}_{j}=0)=\epsilon$. Given $\mathcal{A}_{t}$, one can generate $\mathcal{A}_{t+1}$ from $\mathcal{A}_{t}$ using the state crossover probability. Note that the crossover probability $\epsilon$ tends to be small in practice. The dynamic graph sampling problem444If $\mathcal{S}_{t}$ at time $t$ is still available at $t+1$, we make no changes. While it is possible to swap newly available nodes at time $t+1$ with selected nodes in $\mathcal{S}_{t}$, we found the performance-complexity tradeoff not worthwhile. can be stated formally as follows: Given a sampling set $\mathcal{S}_{t}$ at time $t$, where $\mathcal{S}_{t}\subset\mathcal{A}_{t}$, how to compute $\mathcal{S}_{t+1}$ if $\mathcal{S}_{t}\not\subset\mathcal{A}_{t+1}$?

We next discuss our sampling strategy for this problem.

V-A Sampling on Dynamic Subsets

For conciseness, we denote $\mathcal{P}=\mathcal{A}_{t+1}\cap\mathcal{S}_{t}$ to be the selected and available (SA) set and $\mathcal{U}=\mathcal{S}_{t}\setminus\mathcal{P}$ as the selected and unavailable (SU) node set. Then, our sampling strategy based on node exchange (NE) can be described as follows:

We replace a node $j\in\mathcal{U}$ with an available and unselected (UA) node $k\in\mathcal{Q}$ where $\mathcal{Q}=\mathcal{A}_{t+1}\setminus\mathcal{P}$ is the UA node set. The optimal node $k^{*}$ is greedily chosen if exchanging $j$ with $k^{*}$ minimizes the objective function in all the $(j,k)$ exchange pairs. Repeating this procedure, all unavailable nodes in SU set are replaced by available nodes, which will form a new sampling set $\mathcal{\hat{S}}$. 2. 2.

To further improve performance, node $p\in\mathcal{\hat{S}}$ will be replaced by node $q\in\mathcal{H}$ if the $(p,q)$-pair exchange can induce a lower objective, where $\mathcal{H}=\mathcal{A}_{t+1}\setminus\mathcal{\hat{S}}$. 3. 3.

Step 2) will return the final sampling set $\mathcal{S}_{t+1}$ until it achieves the iteration number constraint, or it could not find one more node-pair exchange to make the objective decrease in this constraint.

After each node-pair exchange, $\mathbf{G}^{-1}_{\mathcal{S}}$ has to be computed. In the following, we propose one method to relieve this computation burden, in which the node-pair is generally indexed by $(j,k)$.

V-B One-pair Node Exchange

The complexity of computing $\mathbf{G}^{-1}_{\mathcal{S}}$ after each node-pair exchange can be reduced based on Sherman-Morrison formula [28], reviewed below.

Lemma 3**.**

Suppose $\mathbf{A}\in\mathbb{R}^{N\times N}$ is an invertible square matrix and $\mathbf{u},\mathbf{v}\in\mathbb{R}^{N}$ are column vectors. If $\mathbf{A}+\mathbf{u}\mathbf{v}^{\top}$ is invertible, then its inverse is given by

[TABLE]

After exchanging node-pair $(j,k)$, we will have a new sample set $\mathcal{\tilde{S}}=\mathcal{S}\cup\left\{k\right\}\setminus\left\{j\right\}$.

It is observed that under some permutation

[TABLE]

and

[TABLE]

$\mathbf{G}_{\mathcal{\tilde{S}}}$ and $\mathbf{G}_{\mathcal{S}}$ share the same sub-matrix of dimension $(M-1)\times(M-1)$. Hence, we can construct $\mathbf{G}_{\mathcal{\tilde{S}}}$ by modifying only the column and row in $\mathbf{G}_{\mathcal{S}}$ corresponding to node $j$. Assuming that node $j$ is the $i$-th node in $\mathcal{S}$, i.e, $j=\mathcal{S}(i)$, we first replace the $i$-th row in $\mathbf{G}_{\mathcal{S}}$ with $\mathbf{G}^{\top}_{\mathcal{\tilde{S}},\{k\}}$ to form an intermediate matrix $\mathbf{F}$, and then replace the $i$-th column of $\mathbf{F}$ with the column vector $\mathbf{G}_{\mathcal{\tilde{S}},\{k\}}$ to construct $\mathbf{G}_{\mathcal{\tilde{S}}}$. Therefore, we can express $\mathbf{G}_{\tilde{\mathcal{S}}}$ as $\mathbf{G}_{\mathcal{S}}$ with two rank-1 updates:

[TABLE]

where $\mathbf{p}=\mathbf{G}_{{\mathcal{\tilde{S}}},\{k\}}-\mathbf{G}_{\mathcal{S},\{j\}}$, $\mathbf{q}=\mathbf{p}-(\mathbf{G}_{kk}-\mathbf{G}_{jj})\mathbf{e}_{i}$ and $\mathbf{e}_{i}$ is an indicator vector with $\mathbf{e}_{i}(i)=1$ and $\mathbf{e}_{i}(j)=0$ for $j\neq i$. As introduced, $\mathbf{F}=\mathbf{G}_{\mathcal{S}}+\mathbf{e}_{i}\mathbf{p}^{\top}$.

From Lemma 3, the inverse of $\mathbf{G}_{\mathcal{\tilde{S}}}$ can be computed via:

[TABLE]

[TABLE]

composed of simple matrix-vector products only.

We compute the initial set $\mathcal{S}_{0}$ using our proposed GFS algorithm described in Section IV. The details of the dynamic subset sampling strategy are outlined in Algorithm 2, called GFS-NE sampling. $K_{0}$ is the upper bound of exchange iterations to contain complexity. In simulations, we will compare this proposed dynamic subset sampling algorithm to existing greedy sampling algorithms which select nodes from scratch at each time $t$ from time-varying subsets.

VI Accompanied Graph Signal Reconstruction

Extending our derivation in sampling, we propose a biased signal reconstruction method with lower variance and complexity than the conventional LS estimator.

VI-A Graph Filter Sub-matrix-based Reconstruction

Recall that the LS estimator [14] is $\hat{\mathbf{x}}_{\text{LS}}={\mathbf{V}_{K}}[{{{(\mathbf{C}{\mathbf{V}_{K}})}^{\top}}\mathbf{C}{\mathbf{V}_{K}}}]^{-1}{(\mathbf{CV}_{K})^{\top}}\mathbf{y}_{\mathcal{S}}$, where $\mathbf{y}_{\mathcal{S}}$ is an observed noise-corrupted signal. The LS reconstruction requires computation of the first $K$ eigenvectors of the graph operator $\mathbf{L}$ and a matrix inverse operation, and it achieves lowest MSE value among unbiased recovery solution. We next derive a biased signal reconstruction method.

Proposition 3**.**

Adding a small shift ${\beta}>0$ into the LS solution, we obtain a biased estimator $\hat{\mathbf{x}}={\mathbf{V}_{K}}[{{{(\mathbf{C}{\mathbf{V}_{K}})}^{\top}}\mathbf{C}{\mathbf{V}_{K}}}+{\beta}\mathbf{I}]^{-1}{(\mathbf{CV}_{K})^{\top}}\mathbf{y}_{\mathcal{S}}$, which can be approximated by

[TABLE]

where $\mathbf{T}^{\emph{{FGFT}}}$ was obtained during GFS sampling using $\mathbf{T}^{\textrm{FGFT}}=\tilde{\mathbf{V}}_{K}\tilde{\mathbf{V}}^{\top}_{K}$, and ${\mathbf{H}}=\mathbf{T}^{\emph{{FGFT}}}+{\beta}\mathbf{I}$.

Proof.

Similar derivation in (8) and (III-A) can show that the shifted LS solution is equal to

[TABLE]

where $\mathbf{P}={{\mathbf{C}}^{\top}}\mathbf{C}{\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}$ and (a) holds since ${\mathbf{V}_{K}}\mathbf{\Psi}^{d}{\mathbf{V}^{\top}_{K}}={\mathbf{V}_{K}}({\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}}\mathbf{C}{\mathbf{V}_{K}})...({\mathbf{V}^{\top}_{K}}{{\mathbf{C}}^{\top}}\mathbf{C}{\mathbf{V}_{K}}){\mathbf{V}^{\top}_{K}}={\mathbf{V}_{K}}{\mathbf{V}^{\top}_{K}}\mathbf{P}^{d}$.

From the definition of $\mathbf{C}$, we know that {\mathbf{C}^{\top}}\mathbf{C}={\left[\begin{array}[]{l}{\mathbf{I}_{\mathcal{S}}}~{}~{}~{}{\bf{0}}\\ {\bf{0}}~{}~{}~{}~{}{\kern 1.0pt}{\bf{0}}\\ \end{array}\right]} under appropriate permutation, which leads to \mathbf{P}=\left[\begin{array}[]{l}{\mathbf{T}_{\mathcal{SV}}}\\ ~{}\mathbf{0}\\ \end{array}\right]. Hence, \mathbf{I}-(\mathbf{P}+{\beta}\mathbf{I})=\left[\begin{array}[]{l}{\mathbf{I}}-({\mathbf{T}_{\mathcal{S}}}+{\beta}\mathbf{I})~{}~{}-{\mathbf{T}_{\mathcal{S}{\mathcal{S}^{c}}}}\\ ~{}~{}~{}~{}~{}~{}\mathbf{0}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(1-{\beta}){\mathbf{I}}\\ \end{array}\right], resulting in

[TABLE]

where “$\bullet$” denotes a nonzero matrix whose dimension is $M\times(N-M)$.

Then, (VI-A) can be written as

[TABLE]

where equality $(a)$ holds from a similar derivation of (11) with a different shift. Finally, we approximate $\mathbf{T}$ by the $\mathbf{T}^{\textrm{FGFT}}$. ∎

Remark: If we choose $\beta=\mu$, then $\mathbf{H}=\mathbf{G}$. All involved matrices in (40) have been obtained in sampling. However, the value for $\mu$ has been designed in the sampling procedure based on the condition number constraint. If $\beta\neq\mu$, we can customize a $\beta$ based on the bias-variance tradeoff (to be discussed), at the cost of computing the inverse $\mathbf{H}^{-1}_{\mathcal{S}}$ once.

We next prove the robustness of the proposed estimator $\mathbf{\hat{x}}$ and propose one strategy to compute the value of $\beta$.

VI-B Robustness Analysis

Assume that the signal has the same energy for different SNRs, and noise $\mathbf{n}$ is i.i.d. with zero mean and variance $\omega^{2}$ which varies with SNR. We write the following proposition.

Proposition 4**.**

Given $\mathbf{T}^{\emph{{FGFT}}}$ sufficiently approximates the ideal LP filter, MSE of the proposed GFS reconstruction is:

[TABLE]

and the MSE of the LS solution is:

[TABLE]

where $(\sigma_{i},\mathbf{u}_{i})$ is the $i$-th eigen-pair of the matrix $\mathbf{\Psi}={\left({\mathbf{C}{\mathbf{V}_{K}}}\right)^{\top}}\mathbf{C}{\mathbf{V}_{K}}$.

Proof.

The proof is detailed in Appendix C. ∎

When the noise variance $\omega^{2}$ is large enough to dominate MSE, the GFS recovery has lower MSE than the LS solution since $\frac{\sigma_{i}}{(\sigma_{i}+\beta)^{2}}<\frac{1}{\sigma_{i}}$, $\forall i$. Therefore the proposed method is more robust to large noise than the LS reconstruction.

The RHS of formula (44) is consisted of bias and variance of estimator $\mathbf{\hat{x}}$ respectively. For the bias part, $\mathbf{u}_{i}$ and $\mathbf{\tilde{x}}_{K}$ are constant in reconstruction, thus smaller $\beta$ will bring smaller bias. For the variance part, smaller $\beta$ will bring larger variance. Hence, given the noise variance $\omega^{2}$, the optimal $\beta$ can be designed to balance the bias-variance tradeoff to achieve the lowest MSE. Instead of computing the gradient of function (44), we propose an empirical and practical strategy to design $\beta$ based on some prior information. For lower MSE, $\beta$ should be comparable to $\sigma_{i}$ to combat large noise and to maintain reasonable bias. We observe that the average value of $\sigma_{i}$ is

[TABLE]

We adopt $\mathbf{T}^{\text{FGFT}}$ to approximate $\mathbf{T}$ and select $\beta$ to be the lower bound of the above average value

[TABLE]

where $|\mathcal{J}|=M$ and $\mathcal{J}\subset\mathcal{V}$. $\forall i\in\mathcal{J}$ and $\forall j\in\mathcal{J}^{c}$, $\mathbf{T}^{\text{FGFT}}_{ii}\leq\mathbf{T}^{\text{FGFT}}_{jj}$. We will validate the performance of this selection of $\beta$ via simulations.

VII Experimental Results

Our proposed sampling and reconstruction strategies are evaluated via experimental simulations. All experiments were performed in MATLAB R2017b, running on a PC with Intel Core I3 3.7 GHz CPU and 16GB RAM. The simulated artificial graphs 555All of these graphs are generated using GSP open source in [36]. are described as follows:

(G1) Community graphs with 1000 nodes and 31 communities;

(G2) Sensor graphs with 1000 nodes.

(G3) Hype-cube graphs with 1002 nodes in 3-dimension.

We also perform experiments on the Minnesota network, which is a real-world graph with $2642$ nodes. Artificial graph signals are assumed to be bandlimited with bandwidth $K=50$. They are constructed by generating appropriate GFT coefficients: the non-zero GFT coefficients are randomly generated from the distribution $\mathcal{N}\left(1,0.5^{2}\right)$, and the coefficients after $K=50$ are zeros. The generated graph signals are corrupted by additional white Guassian noise (AWGN) with different signal-to-noise ratios (SNRs).

VII-A Static Global Sampling

VII-A1 Reconstruction MSE evaluation

The metric used to evaluate a sampling scheme is MSE of the reconstructed signal with respect to the ground truth signal. We compare our proposed GFS sampling algorithm to competing schemes that employ other criteria: SP in [14], E-optimal in [19], MIA in [22] and MFN in [17]. The truncation degree in MIA is $L=10$. The SGWT toolbox [36] is adopted to approximate the ideal LP filter in the MIA method, where $p=25$ and $\alpha=30$ [37]. For SP, the approximation order is $k=10$. As described in Section III, for GFS sampling, the shift parameter $\mu$ is set to be $1/(\kappa_{0}-1)$, where we set $\kappa_{0}=100$ as the condition number constraint. The number of Given rotations matrices is $J=6N\text{log}N$ for G1–G3 and $J=30N\text{log}N$ for the Minnesota graph.

Fig. 1 shows the reconstruction MSE of different sampling strategies in terms of sample size in G1–G3 at different SNRs. As shown in Fig. 1, our proposed GFS achieves lower MSE than three other competing sampling strategies, and closely approximates MFN’s performance in community, sensor and cube graphs at different SNRs. As previously discussed, MFN greedily minimizes MSE one node at a time directly, thus it needs to compute the first $K$ eigenvectors once and perform matrix inverse for each metric evaluation, as detailed in Section IV. In contrast, our proposed strategy can obtain comparable performance with significantly lower complexity.

Fig. 2 shows the reconstruction MSE of different sampling methods for the real-world Minnesota network, showing that our proposed sampling also outperforms SP and E-optimal sampling, where the performance of MIA is too poor to be shown. For visualization, we further present one specific demonstration in Fig. 3. Fig. 3 (a) and (b)–(f) show the original graph signal and the interpolated signal of different sampling methods, respectively. Fig. 3 (b)–(f) show that the proposed GFS achieves the lowest MSE value among competing strategies, and the reconstructed signal is visually smoother with respect to the graph topology compared with signal recovered by SP and MIA.

VII-A2 Numerical complexity comparison

We evaluate the complexity of GFS via numerical simulations. Specifically, we compute the numerical complexity order value in Table I assuming $M=0.05N$. In these experiments, we set $R=10$, $T_{1}=T_{2}(k)=100$ and $|\mathcal{E}|=\mathcal{O}(N)$. Fig. 4 shows the numerical complexity value and shows that GFS has lower complexity than MFN, E-optimal and MIA, especially for large graphs.

VII-A3 Appropriateness of the shift parameter $\mu$

The augmented A-optimality function will well approximate the original criterion (2) if $\mu$ is sufficienty small. However, we customize a $\mu$ based on condition number in the GFS sampling which would bring some approximation error. Table II shows the performance of the designed $\mu=1/(\kappa_{0}-1)$ approaches that of a extremely small $\mu$. This indicates the proposed $\mu$ for stable computation won’t bring a large performance gap between the original A-optimality value and the augmented one. In Section VI, we propose a new shift parameter $\beta$ based on equation (47). Table II tells that if we apply $\mu=\beta$, the sampling performance will also be comparable.

VII-B Static Graph Signal Reconstruction

Simulations on evaluating the proposed GFS reconstruction algorithm are performed subsequently, where the samples were all collected by the GFS sampling algorithm. Fig. 5 shows the reconstruction MSE of different recovery strategies in G1–G3 at 0dB. As illustrated in Fig. 5, the GFS reconstruction outperforms the LS reconstruction for different graphs when the shift parameter $\beta$ is designed by (47), which validates the robustness of the proposed GFS reconstruction compared with the LS solution.

Suppose that $\beta=\mu=1/(\kappa_{0}-1)$, the performance of the proposed GFS sometimes achieves better performance than the LS solution. In this case, the GFS reconstruction does not compute the matrix inverse, thus enjoys lower complexity in recovery, as remarked in Section VI. Based on the experimental results in Fig. 5 and Table II, we can design the shift $\mu$ in sampling and $\beta$ in reconstruction both from equation (47).

VII-C Dynamic Subset Sampling

We perform simulations to evaluate the proposed GFS-NE dynamic sampling strategy. In our experiments, we set $J=18N\text{log}(N)$, $P_{0}=0.8$, $\epsilon=2\%$ and $K_{0}=50$. We first randomly pick up $0.8N$ nodes to form the initial available subset, then this available set evolves over time based on the model defined in Section V 666 There sometimes appears a bad initial subset, then all sampling algorithms perform extremely bad on it, which is meaningless to do sampling work on such subset. We propose one strategy to recognize bad subsets, see appendix D for details, on which we do not collect samples.. For comparison, we simulate the GFS, the E-optimal [19] and SP [14] sampling algorithms, all of which select samples from scratch from new available subset. Fig. 6 shows the reconstruction MSE of different sampling strategies in terms of time evolution at 0dB with fixed sample size. It shows that our proposed GFS-NE ourperforms SP and E-optimal and is comparable to GFS selecting samples from scratch at each time $t$ 777Note that a competing greedy strategy can also replace each unavailable node in $\mathcal{S}_{t}$ one-by-one instead of selecting the entire set at time $t+1$ from scratch. However, it has to compute the initial objective value of the SA set to proceed sampling, and its performance will be necessarily be no better than the “start-over” approach. Hence our method will also outperform competing schemes with unavailable node replacement. .

VIII Conclusion

Graph sampling with noise remains a challenging problem: MMSE leads to the known A-optimality criterion for independent noise, which is expensive to evaluate and difficult to optimize. In this paper, we propose an augmented objective based on Neumann series expansion to approximate the A-optimality criterion, which can be expressed as a function of an ideal LP graph filter, efficiently approximated via fast graph Fourier transform. Using the augmented objective, we select nodes greedily without any matrix inverse computation based on a matrix inversion lemma. Further, we extend our sampling scheme to the dynamic network case, where the availability of nodes is time-varying. For signal recovery, we design an accompanied signal reconstruction strategy to obtain a biased but robust estimator. Experimental results validate the superiority of the proposed sampling strategy compared with existing schemes and demonstrate the effectiveness of our biased recovery algorithm over unbiased LS reconstruction.

Appendix A Proof of Proposition 1

From the definition of $\mathbf{C}$, {\mathbf{C}^{\top}}\mathbf{C}={\left[\begin{array}[]{l}{\mathbf{I}_{\mathcal{S}}}~{}~{}~{}{\bf{0}}\\ {\bf{0}}~{}~{}~{}~{}{\kern 1.0pt}{\bf{0}}\\ \end{array}\right]} under appropriate permutation. Then, $\forall\mathbf{x}\in{\mathbb{R}^{K}}$ and $\left\|\mathbf{x}\right\|_{2}=1$,

[TABLE]

where $\mathbf{b}={{\mathbf{V}_{K}}\mathbf{x}}$.

Since ${\mathbf{b}^{\top}}\mathbf{b}={\left({{\mathbf{V}_{K}}\mathbf{x}}\right)^{\top}}\left({{\mathbf{V}_{K}}\mathbf{x}}\right)=1$, $0\leq{\mathbf{x}^{\top}}\mathbf{\Psi}\mathbf{x}\leq 1$. Because $\textrm{rank}{{{\left(\mathbf{C}{\mathbf{V}_{K}}\right)}}}=K$, $\mathbf{\Psi}$ is positive definite, which results in $0<{\mathbf{x}^{\top}}\mathbf{\Psi}\mathbf{x}\leq 1$ and $0\leq{\mathbf{x}^{\top}}\mathbf{\Phi}\mathbf{x}<1$. Due to the Rayleigh quotient theorem, $0\leq{\delta_{i}}<1$.

Appendix B Proof of Lemma 2

If we view matrix $\tilde{\mathbf{V}}_{K}$ as $\tilde{\mathbf{V}}_{K}=[\mathbf{t}_{1}~{}\mathbf{t}_{2}~{}...\mathbf{t}_{N}]^{\top}$ with $\mathbf{t}_{i}\in\mathbb{R}^{K}$, then $(\mathbf{C}\tilde{\mathbf{V}}_{K})^{\top}\mathbf{C}\tilde{\mathbf{V}}_{K}=\sum_{i\in\mathcal{S}}\mathbf{t}_{i}\mathbf{t}^{\top}_{i}$. Define $\mathbf{Z}(\mathcal{S})=\sum_{i\in\mathcal{S}}\mathbf{t}_{i}\mathbf{t}^{\top}_{i}+\mu\mathbf{I}$, we know that $\mathbf{Z}\left(\mathcal{S}\cup\{j\}\right)=\mathbf{Z}(\mathcal{S})+\mathbf{t}_{j}\mathbf{t}^{\top}_{j}$ for any $j\not\in\mathcal{S}$ and $g(\mathcal{S})=\text{tr}[\mathbf{Z}(\mathcal{S})^{-1}]$.

(i) Monotonic decreasing

[TABLE]

for any $j\not\in\mathcal{S}$.

From Lemma 3, we know that

[TABLE]

It is easy to prove the eigenvalues of $\mathbf{Z}(\mathcal{S})$ are in $(\mu,1+\mu]$ for any $\mathcal{S}$ from the similar derivation in Proposition 1 [22]. Therefore, $\mathbf{t}^{\top}_{j}\mathbf{Z}(\mathcal{S})^{-1}\mathbf{t}_{j}\geq\lambda^{-1}_{\text{max}}[\mathbf{Z}(\mathcal{S})]\cdot\|\mathbf{t}_{j}\|^{2}_{2}\geq\frac{\|\mathbf{t}_{j}\|^{2}_{2}}{1+\mu}>0$ from Rayleigh quotient theorem. Using this result and combining it with the above two equations, we know that $g(\mathcal{S}\cup\{j\})-g(\mathcal{S})\leq 0$ for any $j\not\in\mathcal{S}$ which implies the set function $g$ is monotonic decreasing.

(ii) $\alpha$-supermodular

We first present the definition of $\alpha$-supermodularity introduced in paper [18].

Definition 2**.**

A set function $g:~{}2^{\mathcal{V}}\rightarrow\mathbb{R}$ is $\alpha$-supermodular if for all sets $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{V}$ and all $j\not\in\mathcal{B}$, the following equation holds for some $\alpha\geq 0$

[TABLE]

$\alpha$-supermodularity is only of interest when $\alpha$ takes the largest value [18]:

[TABLE]

Since $\lambda^{-1}_{\text{max}}[\mathbf{Z}(\mathcal{S})]\|\mathbf{t}_{j}\|^{2}_{2}\leq\mathbf{t}^{\top}_{j}\mathbf{Z}(\mathcal{S})^{-1}\mathbf{t}_{j}\leq\lambda^{-1}_{\text{min}}[\mathbf{Z}(\mathcal{S})]\|\mathbf{t}_{j}\|^{2}_{2}$, the lower bound of $\alpha$ is

[TABLE]

As we claimed, the eigenvalues of $\mathbf{Z}(\mathcal{S})$ are in $(\mu,1+\mu]$, so

[TABLE]

It is easy to drive that

[TABLE]

which demonstrates the function $\alpha^{\prime}$ is decreasing in terms of $\|\mathbf{t}_{j}\|^{2}_{2}$. Since $\mathbf{t}_{j}$ is a row of $\mathbf{V}_{K}$, i.e., $\|\mathbf{t}_{j}\|^{2}_{2}\leq 1$, we have

[TABLE]

Hence, there exists $\alpha>0$ to fulfill equation (50) and its value can be bounded by $\frac{\mu^{3}(\mu+2)}{(\mu+1)^{4}}$, which means this function is approximately supermodular with parameter $\alpha$.

Appendix C Proof of Proposition 4

Assume the eigen-decomposition of $\mathbf{\Psi}$ is $\mathbf{\Psi}=\mathbf{U}\mathbf{\Lambda}{{\mathbf{U}^{\top}}}$ where $\mathbf{\Lambda}=\textrm{diag}\{\sigma_{i}\}$ and $\mathbf{U}\mathbf{U}^{\top}=\mathbf{I}$. We have derived in Appendix.B in the paper [22] that

[TABLE]

Rewrite the GFS reconstruction $\hat{\mathbf{x}}={\mathbf{V}_{K}}[{{{(\mathbf{C}{\mathbf{V}_{K}})}^{\top}}\mathbf{C}{\mathbf{V}_{K}}}+\beta\mathbf{I}]^{-1}{(\mathbf{CV}_{K})^{\top}}\mathbf{y}_{\mathcal{S}}$ by $\hat{\mathbf{x}}=\mathbf{H}\mathbf{y}_{\mathcal{S}}=\mathbf{H}\mathbf{C}(\mathbf{x}+\mathbf{n})$, then the bias of this estimator is

[TABLE]

Recall that the original graph signal was assumed $K$-bandlimited, so

[TABLE]

in which

[TABLE]

Therefore,

[TABLE]

The variance part of this solution can be derived from its covariance matrix: $\text{cov}(\mathbf{\hat{x}})=\mathbb{E}\left[(\mathbf{\hat{x}}-\mathbb{E}[\mathbf{\hat{x}}])(\mathbf{\hat{x}}-\mathbb{E}[\mathbf{\hat{x}}])^{\top}\right]=\mathbb{E}\left[\mathbf{HCnn}^{\top}\mathbf{C}^{\top}\mathbf{H}^{\top}\right]$. Because $\mathbb{E}[\mathbf{nn}^{\top}]=\omega^{2}\bf{I}$ and $\mathbf{C}\mathbf{C}^{\top}=\mathbf{I}$, the variance of $\mathbf{\hat{x}}$ is

[TABLE]

From the definition of $\bf{H}$ and the property of trace operator,

[TABLE]

Therefore, the variance of $\mathbf{\hat{x}}$ is

[TABLE]

Appendix D Recognize the Bad Initial Available Subset

Space $PW_{\eta}(G)$ contains all bandlimited graph signals with bandwidth less than $\eta$ and space $L_{2}(\mathcal{S}^{c})$ covers graph signals fulfilling $\mathbf{f}(\mathcal{S})=0$ [15]. Authors in [14] claimed that for set $\mathcal{S}$ to be a uniqueness set of space $PW_{\eta}(G)$, $\eta$ needs to be less than the minimum bandwidth of all signals in space $L_{2}(\mathcal{S}^{c})$, which is defined as the cutoff frequency of $L_{2}(\mathcal{S}^{c})$. Thus, a set $\mathcal{S}$ to be good means its cutoff frequency should be as large as possible. They defined the cutoff frequency of space $L_{2}(\mathcal{S}^{c})$ as

[TABLE]

where $\eta(\mathbf{f})$ is the frequency of a specific graph signal $\mathbf{f}$.

A strategy was proposed in [14] to approximate the cutoff frequency by

[TABLE]

where $\lambda_{1,k}$ is the smallest eigenvalue of the reduced matrix $\left((\mathbf{L}^{\top})^{k}\mathbf{L}^{k}\right)_{\mathcal{S}^{c}}$.

This cutoff frequency can be used for judging if a set is good for initialization because a reasonable set won’t have too small $\Omega_{k}(\mathcal{S})$. We take this as a criterion to detour a bad initial subset. During experiments, the size of the available set is $M=800$, so we randomly pick 50 sets off-line with $M=800$ to compute their approximate cutoff frequency $\Omega_{k}(\mathcal{S})$ and sort these values in ascending order. We take the fifth value in this rank, denoted by $\tilde{\Omega}$, as the bound for define a bad set. For an initial set $\mathcal{S}$, once $\Omega_{k}(\mathcal{S})>\tilde{\Omega}$, we define it to be a good set and then perform sampling on it. The following available sets are generated from this initial set with small crossover probability, so their reliability can be guaranteed.