Fast Graph Sampling Set Selection Using Gershgorin Disc Alignment

TL;DR

This paper introduces a fast, eigen-decomposition-free graph sampling method based on Gershgorin disc alignment, improving efficiency and accuracy in reconstructing smooth graph signals.

Contribution

It proposes a novel Gershgorin disc alignment approach for sample selection, avoiding expensive eigenvector computations and providing a fast approximation scheme.

Findings

Runs substantially faster than existing schemes

Outperforms other eigen-decomposition-free methods in reconstruction error

Effective in selecting samples for smooth graph signal reconstruction

Abstract

Graph sampling set selection, where a subset of nodes are chosen to collect samples to reconstruct a smooth graph signal, is a fundamental problem in graph signal processing (GSP). Previous works employ an unbiased least-squares (LS) signal reconstruction scheme and select samples via expensive extreme eigenvector computation. Instead, we assume a biased graph Laplacian regularization (GLR) based scheme that solves a system of linear equations for reconstruction. We then choose samples to minimize the condition number of the coefficient matrix---specifically, maximize the smallest eigenvalue $\lambda_{\min}$. Circumventing explicit eigenvalue computation, we maximize instead the lower bound of $\lambda_{\min}$, designated by the smallest left-end of all Gershgorin discs of the matrix. To achieve this efficiently, we first convert the optimization to a dual problem, where we minimize the number of samples needed to align all Gershgorin disc left-ends at a chosen lower-bound target $T$. Algebraically, the dual problem amounts to optimizing two disc operations: i) shifting of disc centers due to sampling, and ii) scaling of disc radii due to a similarity transformation of the matrix. We further reinterpret the dual as an intuitive disc coverage problem bearing strong resemblance to the famous NP-hard set cover (SC) problem. The reinterpretation enables us to derive a fast approximation scheme from a known SC error-bounded approximation algorithm. We find an appropriate target $T$ efficiently via binary search. Extensive simulation experiments show that our disc-based sampling algorithm runs substantially faster than existing sampling schemes and outperforms other eigen-decomposition-free sampling schemes in reconstruction error.

Tables

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

	E-optimal	SP	MFN	MIA	Ed-free	proposed BS-GDA
Preparation	$𝒪\left(\left(MK+R{K}^3\right)T_1\right)$	$𝒪\left(kMK{T}_2\right)$	$𝒪\left(\left(MK+R{K}^3\right)T_1\right)$	$𝒪\left(qNM\right)$	$𝒪\left(M\overline{q}+N\overline{Q}\right)$	$𝒪\left(N\left(P+Q\right){\log }_2\frac{1}{ϵ}\right)$
Sampling	$𝒪\left(N{K}^4\right)$	$𝒪\left(NK\right)$	$𝒪\left(N{K}^4\right)$	$𝒪\left(NF{K}^{3.373}\right)$	$𝒪\left(NK\overline{Q}\right)$	$𝒪\left(NKP{\log }_2\frac{1}{ϵ}\right)$

Table 2. TABLE II: Speedup factors of our algorithm with respect to other sampling algorithms for N = 3000 𝑁 3000 N=3000

Sampling Algorithms	Sensor	Community
Random [27]	0.22	0.21
E-optimal [24]	2812.77	1360.76
SP [16]	174.09	466.18
MFN [22]	2532.91	1184.23
MIA [20]	1896.19	964.65
Ed-free [9]	1.82	8.11

Equations

$$\mathbf{L}\triangleq \mathbf{D}-\mathbf{W}$$

$$\mathbf{L}=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^{\top }$$

$$\mathbf{\alpha }={\mathbf{U}}^{\top }\mathbf{x}=[{\mathbf{u}}_1,...,{\mathbf{u}}_N{]}^{\top }\mathbf{x}$$

$$\mathbf{y}={\mathbf{x}}_{\mathcal{S}}+\mathbf{n}=\mathbf{H}\mathbf{x}+\mathbf{n}$$

$${\mathbf{H}}_{ij}=\{\begin{array}{l}1,j={\mathcal{S}}_i;\\ 0,\text{mbox}otherwise.\end{array}$$

$$L_2(\mathcal{S})=\{\mathbf{x}\in {\mathbb{R}}^N|{\mathbf{x}}_{{\mathcal{S}}^c}=0\}$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\Vert \mathbf{H}\mathbf{x}-\mathbf{y}{\Vert }_2^2+\mu {\mathbf{x}}^{\top }\mathbf{L}\mathbf{x}$$

$$({\mathbf{H}}^{\top }\mathbf{H}+\mu \mathbf{L})\hat{\mathbf{x}}={\mathbf{H}}^{\top }\mathbf{y}.$$

$${\lambda }_{\min }$$

$$=\mu \frac{{\alpha }_1^2{\theta }_1+\dots +{\alpha }_k^2{\theta }_k}{\Vert \mathbf{\alpha }{\Vert }_2^2}\le \mu \omega .$$

$$a_i=\{\begin{array}{l}1,i\in \mathcal{S},\\ 0,\text{mbox}otherwise.\end{array}$$

$$\underset{\mathbf{a}}{\max }$$

$$\mathbf{B}=\mathbf{A}+\mu \mathbf{L},$$

$$\mathbf{A}=\mathrm{diag}(\mathbf{a}),a_i\in \{0,1\},\sum _{i=1}^N{a}_i\le K.$$

$$\hat{\mathbf{x}}$$

$$=({\mathbf{H}}^{\top }\mathbf{H}+\mu \mathbf{L}{)}^{-1}{\mathbf{H}}^{\top }(\mathbf{H}\mathbf{x}+\mathbf{n})$$

$$=({\mathbf{H}}^{\top }\mathbf{H}+\mu \mathbf{L}{)}^{-1}{\mathbf{H}}^{\top }\mathbf{H}(\mathbf{x}+\tilde{\mathbf{n}})$$

$$\hat{\mathbf{x}}={\mathbf{B}}^{-1}(\mathbf{B}-\mu \mathbf{L})(\mathbf{x}+\tilde{\mathbf{n}})=(\mathbf{I}-\mu {\mathbf{B}}^{-1}\mathbf{L})(\mathbf{x}+\tilde{\mathbf{n}})$$

$$\Vert \hat{\mathbf{x}}-\mathbf{x}{\Vert }_2$$

$$=\Vert -\mu {\mathbf{B}}^{-1}\mathbf{L}(\mathbf{x}+\tilde{\mathbf{n}})+\tilde{\mathbf{n}}{\Vert }_2$$

$$\le \Vert -\mu {\mathbf{B}}^{-1}\mathbf{L}(\mathbf{x}+\tilde{\mathbf{n}}){\Vert }_2+\Vert \tilde{\mathbf{n}}{\Vert }_2$$

$$\le \mu \Vert {\mathbf{B}}^{-1}{\Vert }_2\Vert \mathbf{L}(\mathbf{x}+\tilde{\mathbf{n}}){\Vert }_2+\Vert \tilde{\mathbf{n}}{\Vert }_2$$

$${\Psi }_i(x_{ii},R_i)=\{z\in \mathbb{C}||z-x_{ii}|\le R_i\},i\in \{1,\dots ,n\}.$$

$$b_{ii}-R_i\le \lambda \le b_{ii}+R_i,$$

$$\underset{i}{\min }\{b_{ii}-R_i\}=\underset{i}{\min }{a}_i\le {\lambda }_{\min }(\mathbf{B}).$$

$$a_i:0\to 1,i\in \mathcal{S}$$

$$\mathbf{C}=\mathbf{S}\mathbf{B}{\mathbf{S}}^{-1},$$

$$\underset{\mathbf{a},\mathbf{s}}{\max }$$

$$\mathbf{C}=\mathbf{S}\left(\mathbf{A}+\mu \mathbf{L}\right){\mathbf{S}}^{-1}$$

$$\mathbf{A}=\mathrm{diag}(\mathbf{a}),a_i\in \{0,1\},\sum _{i=1}^N{a}_i\le K,$$

$$\mathbf{S}=\mathrm{diag}(\mathbf{s}),s_i>0.$$

$$\underset{\mathbf{a},\mathbf{s}}{\min }$$

$$\mathbf{C}=\mathbf{S}\left(\mathbf{A}+\mu \mathbf{L}\right){\mathbf{S}}^{-1},c_{ii}-\sum _{j\ne i}|c_{ij}|\ge T,\forall i$$

$$\mathbf{A}=\mathrm{diag}(\mathbf{a}),a_i\in \{0,1\},$$

$$\mathbf{S}=\mathrm{diag}(\mathbf{s}),s_i>0.$$

$$\underset{i}{\min }{\ell }_i\le {\lambda }_{\min }(\mathbf{B})\le \underset{i}{\max }{\ell }_i$$

$$\underset{\mathbf{s}}{\max }$$

Full text

I Introduction

Sampling of smooth signals—and subsequent signal reconstruction from collected samples—are fundamental problems in signal processing. Graph sampling extends the well-known studies of Nyquist sampling on regular data kernels to irregular data kernels described by graphs. Graph sampling can be studied under different contexts: aggregation sampling [1, 2], local measurement [3] and subset sampling [4]. In this paper, we focus on graph sampling set selection: a subset of nodes are chosen to collect samples to reconstruct a graph signal with a well-defined notion of smoothness.

For the noiseless case, graph sampling set selection is well investigated [5, 6, 7]; in particular, [7] shows that random selection of $K$ samples on a graph for a target $K$-bandlimited signal can enable perfect reconstruction with high probability. The noisy case where acquired samples are corrupted by additive noise is more challenging, and there exists many proposals. Assuming that the target signal is bandlimited, a class of schemes compute one or more extreme eigenvectors of a graph Laplacian matrix (or sub-matrix) and examine the per-node magnitudes to greedily determine a sample set. While these methods show good sampling performance using a least-squares (LS) signal reconstruction scheme, their complexity is high. On the other hand, there are faster sampling methods based on the notion of local operators [8, 9], but these methods are ad-hoc in nature, and there are no global notions of optimality associated with the selected sampling sets.

In practice, real-world graph signals, such as temperature in sensor network [10], category labels in graph-based classifier [11] and ratings in recommendation system [12], are typically not strictly bandlimited with known bandwidths. Instead, generally “smooth” signals with respect to known graph structures are relatively more common. Different from standard least-squares reconstruction only considering spectrum energy under a given cutoff frequency, we assume a different signal smoothness notion defined using graph Laplacian regularization (GLR), and then propose a corresponding sampling scheme with roughly linear complexity. Further, our scheme can be proven to minimize the global worst case mean squared error (MSE); thus our sampling scheme can gracefully scale to very large graphs and has excellent signal reconstruction performance. Specifically, we first assume a biased GLR based scheme that solves a system of linear equations for reconstruction given collected samples [13]. We then choose samples to minimize the condition number of the coefficient matrix; i.e., maximize the smallest eigenvalue $\lambda_{\min}$. Circumventing computation of extreme eigenvalues, we maximize instead the lower bound of $\lambda_{\min}$, designated by the smallest left-end of all Gershgorin discs of the matrix [14].

To achieve this efficiently, we first convert the optimization to a dual problem, where we minimize the number of samples needed to align all Gershgorin disc left-ends at a pre-selected lower-bound target $T$. Algebraically, the dual problem amounts to two disc operations: i) shifting of disc centers due to sampling, and ii) scaling of disc radii due to a similarity transformation of the matrix. We further reinterpret the dual as an intuitive disc coverage (DC) problem, which bears strong resemblance to the famous set cover (SC) problem [15]. Though DC remains NP-hard (we provide a proof by reduction), the reinterpretation enables us to derive a fast approximation scheme from a known SC error-bounded approximation algorithm [15]. An appropriate target $T$ can be located quickly via simple binary search. Extensive simulation experiments show that our disc-based sampling algorithm runs substantially faster than existing eigen-decomposition-based schemes (hundreds to thousands times faster) for medium-sized graphs and outperforms other eigen-decomposition-free schemes in reconstruction error for large-sized graphs.

The outline of the paper is as follows. We first overview related works in Section II. We introduce notations and a biased reconstruction method in Section III. In Section IV, we first propose our sampling formulation to minimize the upper bound of MSE and then derive its dual problem based on Gershgorin circle theorem. We develop a fast sampling set selection algorithm to tackle our sampling formulation with roughly linear complexity in Section V. Experimental results are presented in Section VI. We conclude in Section VII.

II Related Works

The problem of graph sampling can be defined from different contexts. From an aggregation sampling context, the authors in [1, 2] proposed to aggregate the sampled signal at one node (or a subset) but with the graph signal shifted successively by a known shift operator. From a local weighted sampling context, [3] proposed to collect samples via a weighted sum of signals from different local sets. In contrast, we focus on graph sampling in the context of set selection in this paper: how to select a subset of nodes to collect samples to reconstruct a smooth signal in high fidelity. We next focus on detailed reviews of related works under the umbrella of set selection-based graph sampling.

II-A Deterministic Graph Sampling Set Selection

Sampling theory for noise-free graph signals was first studied in [5], where a uniqueness set is defined for perfect recovery of bandlimited graph signals. Starting from this definition, [16, 17, 18] proved a sufficient and necessary condition for a uniqueness set, and proposed a lightweight sampling set selection strategy based on the notion of graph spectral proxies.

In the presence of noise, assuming an unbiased least-squares (LS) signal reconstruction with independent and identically distributed (i.i.d) noise model, the minimum mean square error (MMSE) function for sampling is actually the known A-optimal criterion [19, 20, 21]. In response, the authors in [22, 23] proposed to directly optimize the A-optimal function with greedy scheme. Alternatively, [24] used the E-optimality criterion for graph sampling, which can be interpreted as minimizing the worst case of MSE. However, the criteria in [22, 23, 24] all required extreme eigen-pair computation of the graph variation operator to choose samples, which is not scalable for large graphs. Though [20] avoided eigen-pair computation by developing a sampling method based on truncated Neumann series, it must compute many matrix-series multiplications for good approximation, which is expensive.

To tackle the complexity issue in large scale graphs, very recently, [8, 9] proposed an eigen-decomposition-free graph sampling method by successively maximizing the coverage of the localization operators, which was implemented via Chebyshev polynomial approximation [25]. Nonetheless, this sampling idea utilizes local message heuristically and has no global performance guarantee. Orthogonally, our early work [26] proposed a fast graph sampling algorithm via Gershgorin disc alignment without any eigen-pair computation, but its performance was sub-optimal due to the node sampling strategy simply based on breadth first search (BFS). It is known that random sampling [27, 28] can lead to very low computational complexity, but it typically requires more samples for the same signal reconstruction quality compared to its deterministic counterparts.

In this paper, we propose a roughly linear-time deterministic graph sampling method without any eigen-pair computation, while its reconstruction performance remains superior. Unlike previous fast sampling works [16, 9, 26], our method can be proven to optimize the robustness of the linear reconstruction system resulting from GLR and reduce the upper-bound of the global reconstruction MSE, while each sampling step is performed locally for speed.

II-B Sampling of GMRF Graph Signals

Besides the bandlimitedness assumption, graph signal (label information) in classical semi-supervised learning is interpreted as a Gaussian Markov random field (GMRF), where the signal is smooth with respect to the graph Laplacian operator [29]. With this property, active semi-supervised learning can be interpreted as graph sampling on a nearest neighbor graph constructed by node features [11]. An adaptive graph sampling method was proposed [30] to collect samples based on both successively modified graph and selected samples. To improve efficiency, authors in [31] proposed a hybrid adaptive and non-adaptive sampling method to combine the graph sampling and active learning methods. Under the GMRF model, [32] showed that the MSE function for graph sampling is the trace of the posterior covariance matrix, based on which they developed a greedy sampling method. However, this method is very inefficient in large graphs because of the matrix inverse computation in each greedy step. In this paper, we derive a fast sampling strategy for GMRF graph signals with smoothness prior, but it is also applicable for bandlimited graph signals.

III Preliminary

III-A Background and Notations

A graph ${\cal G}$ is defined by a triplet ${\cal G}({\cal V},{\cal E},{\mathbf{W}})$, where ${\cal V}$ and ${\cal E}$ represent sets of $N$ nodes and $M$ edges in the graph, respectively. Associated with each edge $(i,j)\in{\cal E}$ is a weight $w_{ij}$, which reflects the correlation or similarity between two connected nodes $i$ and $j$. In this paper, we assume a connected, undirected graph; i.e., $w_{ij}=w_{ji},\forall i,j\in{\cal V}$. ${\mathbf{W}}$ is an adjacency matrix with $w_{ij}$ as the $(i,j)$-th entry of the matrix. We assume a non-negative weight matrix with elements $w_{ij}>0$ if $\forall(i,j)\in{\cal E}$, and $w_{ij}=0$ otherwise.

Given ${\mathbf{W}}$, the combinatorial graph Laplacian matrix ${\mathbf{L}}$ is computed as [33, 34]:

[TABLE]

where $\mathbf{D}=\operatorname*{diag}({\mathbf{W}}\mathbf{1})$ is a diagonal degree matrix. $\mathbf{1}$ is a vector of all 1’s and $\operatorname*{diag}(\cdot)$ is an operator that returns a diagonal matrix with the elements of an input vector on the main diagonal. As ${\mathbf{L}}$ is real and symmetric, it can be eigen-decomposed via the Spectral Theorem [35],

[TABLE]

where ${\boldsymbol{\Lambda}}$ is a diagonal matrix containing real eigenvalues $\theta_{k}$ along its diagonal, and ${\mathbf{U}}$ is an orthonormal matrix composed of real eigenvectors ${\mathbf{u}}_{k}$’s as columns. Since $w_{ij}\geq 0$ for $\forall(i,j)\in{\cal E}$, ${\mathbf{x}}^{\top}{\mathbf{L}}{\mathbf{x}}=\sum_{i}\sum_{j}w_{ij}(x_{i}-x_{j})^{2}\geq 0$ for $\forall{\mathbf{x}}\in{\mathbb{R}}^{N}$. Thus, ${\mathbf{L}}$ is a positive semi-definite (PSD) matrix. The non-negative eigenvalues $\theta_{k}$’s are interpreted as graph frequencies. By sorting $\theta_{k}$’s in a non-decreasing order, i.e., $0=\theta_{1}\leq\theta_{2}\leq\ldots\leq\theta_{N}$, $\theta_{k}$’s are ordered from lowest to highest frequencies. $\theta_{1}$ represents the lowest DC graph frequency and $\theta_{N}$ represents the highest AC graph frequency. The corresponding eigenvectors in ${\mathbf{U}}$ are interpreted as Fourier-like graph frequency basis; together they form the graph spectrum for a graph ${\cal G}$.

Given a graph frequency basis ${\mathbf{U}}$, a graph signal ${\mathbf{x}}$ can be transformed to its graph frequency domain representation ${\boldsymbol{\alpha}}$ via graph Fourier transform (GFT) [33, 34]:

[TABLE]

A graph signal ${\mathbf{x}}$ is called $\omega$-bandlimited if $\alpha_{i}=0$, for $\forall i\in\{i~{}|~{}\theta_{i}>\omega\}$.

A graph signal ${\mathbf{x}}\in{\mathbb{R}}^{N}$ is defined on the nodes of a graph, i.e., one scalar value $x_{i}$ is assigned to the $i$-th node. Given a graph signal ${\mathbf{x}}$ and subset ${\cal S}\subset\{1,...,N\}$ with cardinality $|{\cal S}|=K<N$, we define the sampling process as [5, 36, 24, 16, 27, 20, 22, 23, 9, 7]:

[TABLE]

where ${\mathbf{y}}\in{\mathbb{R}}^{K}$ is a sampled observation of length $K$. ${\mathbf{x}}_{\cal S}\in{\mathbb{R}}^{K}$ denotes a sub-vector of ${\mathbf{x}}$ consisting of elements indexed by ${\cal S}$ and ${\mathbf{n}}$ is an additive noise term. The sampling matrix ${\mathbf{H}}\in\{0,1\}^{K\times N}$ is defined as [20]:

[TABLE]

where ${\cal S}_{i}$ denotes the $i$-th element of set ${\cal S}$. The complement set of ${\cal S}$ on ${\cal V}$ is denoted by ${\cal S}^{c}={\cal V}\setminus{\cal S}$. We define $L_{2}({\cal S})$ as the space of all graph signals that are zero except on the subset of nodes indexed by ${\cal S}$ [16]:

[TABLE]

III-B Signal Reconstruction with Graph Laplacian Regularizer

Instead of assuming an unbiased least-squares (LS) signal reconstruction scheme from sparse samples [24, 22, 23, 20], we adopt a popular biased graph signal reconstruction scheme based on graph Laplacian regularization (GLR) [13], which has demonstrated its effectiveness in numerous practical applications, such as image processing [13, 37, 38, 39, 40, 41] , computer graphics [42] and semi-supervised learning [43, 44]. Biased signal reconstruction using GLR leads to an unconstrained $l_{2}$-$l_{2}$-norm minimization, whose solution can be obtained by solving a system of linear equations, efficiently computed using state-of-the-art numerical linear algebra algorithms such as conjugate gradient (CG) [45]. For large noise variance, it is shown that the biased estimation has smaller reconstructed errors on average than the unbiased counterpart [20, 46].

Specifically, given a noise-corrupted sampled observation ${\mathbf{y}}\in{\mathbb{R}}^{K}$ on a graph ${\cal G}$, one can formulate an optimization for the target signal $\hat{{\mathbf{x}}}\in\mathbb{R}^{N}$ using GLR as follows:

[TABLE]

where ${\mathbf{H}}$ is a sampling matrix defined in (5). $\mu>0$ is a tradeoff parameter to balance GLR against the $l_{2}$-norm data fidelity term. The law of selecting an optimal $\mu$ with respect to both noise variance and graph spectrum is studied in [46]. The focus of this paper is on graph sampling rather than graph signal reconstruction, and interested readers can refer to [46] for a discussion on the appropriate selection of $\mu$.

Since objective (7) is quadratic, the optimal solution ${\hat{\mathbf{x}}}$ can be obtained by solving a system of linear equations:

[TABLE]

For ${\mathbf{L}}$ constructed on a connected graph ${\cal G}$, DC graph frequency ${\mathbf{u}}_{1}=\mathbf{1}/{\sqrt{N}}$ leads to ${\mathbf{u}}^{\top}_{1}{\mathbf{L}}{\mathbf{u}}_{1}=0$ and $\|{\mathbf{H}}{\mathbf{u}}_{1}\|_{2}>0$ for $\forall{\cal S}$ with $K\geq 1$; for any nonzero graph signal ${\mathbf{x}}\in{\mathbb{R}}^{N}\neq{\mathbf{u}}_{1}$, ${\mathbf{x}}^{\top}{\mathbf{L}}{\mathbf{x}}>0$ and $\|{\mathbf{H}}{\mathbf{x}}\|_{2}\geq 0$. Therefore, $\forall{\mathbf{x}}\neq\mathbf{0}$, ${\mathbf{x}}^{\top}({\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}){\mathbf{x}}=\|{\mathbf{H}}{\mathbf{x}}\|^{2}_{2}+\mu{\mathbf{x}}^{\top}{\mathbf{L}}{\mathbf{x}}>0$. Hence, coefficient matrix ${\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}$ must be positive definite (PD), and (8) has a unique solution ${\hat{\mathbf{x}}}=({\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}})^{-1}{\mathbf{H}}^{\top}{\mathbf{y}}$.

IV Problem Formulation

IV-A Reconstruction Stability and Sampling Set

Although ${\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}$ is a PD matrix, (8) can potentially be poorly conditioned, because the minimum eigenvalue $\lambda_{\min}$ of ${\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}$ can be very close to zero. For instance, if ${\mathbf{H}}$ is chosen such that there exists a low frequency graph signal ${\mathbf{x}}\in L_{2}({\cal S}^{c})$ and its bandlimit $\omega$ is very small, then

[TABLE]

where $\mathop{=}\limits^{(1)}$ holds because of the definition in (6). For small $\mu$, $\lambda_{\min}$ is then upper-bounded by a small number.

This means that condition number $\rho=\lambda_{\max}/\lambda_{\min}$ for coefficient matrix ${\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}$ can be very large, resulting in instability when $\hat{{\mathbf{x}}}$ is numerically computed in (8), and more importantly, $\hat{{\mathbf{x}}}$ is sensitive to the noise in the observation ${\mathbf{y}}$. Our goal then is to maximize $\lambda_{\min}$ by appropriately choosing sampling matrix ${\mathbf{H}}$, so that $\rho$ is minimized111Largest eigenvalue $\lambda_{\max}$ is upper-bounded by $1+2\mu d_{\max}\ll\infty$, where $d_{\max}$ is the maximum node degree of the graph $\mathcal{G}$. Thus to minimize condition number $\rho=\lambda_{\max}/\lambda_{\min}$, it is sufficient to maximize $\lambda_{\min}$..

Denote a diagonal matrix by ${\mathbf{A}}={\mathbf{H}}^{\top}{\mathbf{H}}=\operatorname*{diag}({\mathbf{a}})=\operatorname*{diag}(a_{1},...,a_{N})\in{\mathbb{R}}^{N\times N}$, where sampling vector ${\mathbf{a}}$ satisfies

[TABLE]

In other words, ${\mathbf{a}}$ is [math]-$1$ bit array encoded sampling set ${\cal S}$. Denote by ${\mathbf{B}}={\mathbf{A}}+\mu{\mathbf{L}}$.

Given matrix ${\mathbf{L}}$, parameter $\mu$ and a sampling budget $K$ where $0<K<N$, our objective is to find a sampling set ${\cal S}$ with $|{\cal S}|\leq K$ (or, equivalently, a sampling vector ${\mathbf{a}}$ with $\sum_{i=1}^{N}a_{i}\leq K$), to maximize $\lambda_{\min}$ of ${\mathbf{B}}$, i.e.,

[TABLE]

We next show that maximizing $\lambda_{\min}({\mathbf{B}})$ in (11) is equivalent to minimizing the worst case of MMSE.

Proposition 1**.**

Maximizing $\lambda_{\min}({\mathbf{B}})$ in (11) minimizes the upper bound of MSE between original signal ${\mathbf{x}}$ and reconstructed signal ${\hat{\mathbf{x}}}$.

Proof.

Reconstructed signal ${\hat{\mathbf{x}}}$ from a noisy sampled signal ${\mathbf{y}}$ using (8) is

[TABLE]

where ${\tilde{\mathbf{n}}}\in{\mathbb{R}}^{N}$ is an $N$-dimensional additive noise. Denote by ${\mathbf{B}}={\mathbf{A}}+\mu{\mathbf{L}}={\mathbf{H}}^{\top}{\mathbf{H}}+\mu{\mathbf{L}}$. We have

[TABLE]

To minimize MSE between original signal ${\mathbf{x}}$ and reconstructed signal ${\hat{\mathbf{x}}}$ is equivalent to minimizing $\|{\hat{\mathbf{x}}}-{\mathbf{x}}\|_{2}$.

[TABLE]

Since ${\mathbf{B}}$ is a symmetric PD matrix, $\|{\mathbf{B}}^{-1}\|_{2}=\frac{1}{\lambda_{\min}({\mathbf{B}})}$. Thus, maximizing $\lambda_{\min}({\mathbf{B}})$ minimizes the upper bound of MSE, given fixed $\mu$, ${\mathbf{L}}$ and ${\mathbf{x}}$. ∎

IV-B Maximizing the lower bound of $\lambda_{\min}({\mathbf{B}})$

Solving (11) directly is a challenging combinatorial optimization problem. To begin with, for large graphs, computing $\lambda_{\min}({\mathbf{B}})$ given ${\mathbf{a}}$ using a state-of-the-art numerical linear algebra algorithm like the Lanczos algorithm [47] or the LOBPCG algorithm [48] is still computation-expensive. Moreover, we need to compute $\lambda_{\min}({\mathbf{B}})$ for exponential number of candidate sampling vectors ${\mathbf{a}}$ where $\sum_{i=1}^{N}a_{i}=K$.

In order to study $\lambda_{\min}({\mathbf{B}})$ without performing expensive eigen-pair computation, we design an efficient algorithm based on Gershgorin Circle Theorem (GCT) [35]. We first describe GCT below:

Theorem 1** (Gershgorin Circle Theorem).**

Let ${\mathbf{X}}$ be a complex $n\times n$ matrix with entries $x_{ij}$. For $i\in\{1,...,n\}$, let $R_{i}=\sum_{j\neq i}|x_{ij}|$ be the sum of the absolute values of the off-diagonal entries in the $i$-th row. Consider the $n$ Gershgorin discs

[TABLE]

where $x_{ii}$ and $R_{i}$ are the center and radius of disc $i$, respectively. Each eigenvalue of ${\mathbf{X}}$ lies within at least one Gershgorin disc $\Psi_{i}(x_{ii},R_{i})$.

For our specific real-valued coefficient matrix ${\mathbf{B}}$, each real eigenvalue $\lambda$ lies within at least one Gershgorin disc $\Psi_{i}(b_{ii},R_{i})$ with disc center $b_{ii}$ and radius $R_{i}$, i.e.,

[TABLE]

where $R_{i}=\sum_{j\neq i}|b_{ij}|=\mu\sum_{j}w_{ij}=\mu d_{i}$, and $d_{i}$ is the degree of node $i$. The third equality is true since there are no self-loops in ${\cal G}$. Center of disc $i$ is $b_{ii}=\mu d_{i}+a_{i}$.

Using GCT, we can compute the lower bound of $\lambda_{\min}({\mathbf{B}})$ by examining left-ends of Gershgorin discs. Specifically, the lower bound of $\lambda_{\min}({\mathbf{B}})$ is computed as the smallest left-end of all Gershgorin discs, i.e.,

[TABLE]

For intuition, we examine the two extreme sampling cases. One extreme case is when no nodes are sampled, i.e., ${\cal S}=\varnothing$, and both the lower bound and $\lambda_{\min}$ of ${\mathbf{B}}=\mu{\mathbf{L}}$ are [math]. The other extreme case is when all nodes are sampled, i.e., ${\cal S}={\cal V}$, and both the lower bound and $\lambda_{\min}$ of ${\mathbf{B}}={\mathbf{I}}+\mu{\mathbf{L}}$ are $1$. Note that in both extreme cases, the smallest Gershgorin disc left-end (lower bound) is also exactly $\lambda_{\min}$. The observation inspires us to maximize the lower bound of $\lambda_{\min}({\mathbf{B}})$ instead of maximizing $\lambda_{\min}({\mathbf{B}})$ directly.

In general, the tightness of the Gershgorin lower bound of $\lambda_{\min}$ does not hold for any sampling sets. In order to tighten and maximize the lower bound of $\lambda_{\min}({\mathbf{B}})$, we introduce two basic operations to manipulate Gershgorin discs:

Disc Shifting: The first operation is disc shifting via sampling. In (16), the left-end of the $i$-th Gershgorin disc $\Psi_{i}$ in matrix ${\mathbf{B}}$ is $b_{ii}-R_{i}=a_{i}$. When node $i$ is sampled, we have

[TABLE]

The corresponding left-end of disc $\Psi_{i}$ shifts from [math] to $1$.

Disc Scaling: The second operation is disc scaling via similarity transformation. For $0<K<N$, there is at least one $a_{i}=0$, $i\in\mathcal{S}^{c}$, and the lower bound of $\lambda_{\min}({\mathbf{B}})$ in (17) is always [math]. We introduce disc scaling to pull left-ends of unsampled discs away from [math], resulting in a tighter lower bound for $\lambda_{\min}({\mathbf{B}})$. Specifically, disc scaling is defined as

[TABLE]

where ${\mathbf{S}}$ is a chosen invertible square matrix, named scaling matrix. Matrix ${\mathbf{C}}$ is the similarity transformation of the square matrix ${\mathbf{B}}$ and shares the same eigenvalues with ${\mathbf{B}}$ (but with different eigenvectors) [14]. Thus, we can examine Gershgorin disc left-ends of transformed matrix ${\mathbf{C}}$ instead of original ${\mathbf{B}}$.

To simplify the design of ${\mathbf{S}}$, we employ a PD diagonal matrix ${\mathbf{S}}=\operatorname*{diag}({\mathbf{s}})$, $s_{i}>0$, $\forall i\in\{1,\ldots,N\}$, to scale Gershgorin discs of ${\mathbf{B}}$, so that left-ends of Gershgorin discs of resulting transformed ${\mathbf{C}}$ can be moved right. Specifically, the $i$-th diagonal term $s_{i}$ of ${\mathbf{S}}$ in (19) and its corresponding element $s_{i}^{-1}$ in ${\mathbf{S}}^{-1}$ are used to scale the radius $R_{i}$ of $\Psi_{i}$ and the radii of its neighbors’ discs $\Psi_{j}$ respectively, where $j\in\mathcal{N}_{i}=\{j~{}|~{}w_{ij}>0\}$. For example, if a scalar $s_{i}>1$ is used to expand $R_{i}$, then its neighbors’ discs are shrunk with $s_{i}^{-1}<1$. Since $s_{i}$ is always offset by $s_{i}^{-1}$ on the main diagonal, the center $b_{ii}$ of disc $\Psi_{i}$ is unchanged.

5-node Path Graph Example: To illustrate the two basic operations (shifting and scaling) on Gershgorin discs, as an example we consider a $5$-node path graph and start by sampling node $3$, as shown in Fig. 1. Assuming $\mu=1$, the graph’s coefficient matrix ${\mathbf{B}}$ with $(3,3)$-th entry updated is shown in Fig. 2a. The left-end of disc $\Psi_{3}$ is shifted from [math] to $1$, as shown in Fig. 2d. After disc $\Psi_{3}$ is shifted, we apply a scalar $s_{3}>1$ on the third row of ${\mathbf{B}}$ to expand radius and center of disc $\Psi_{3}$ in Fig. 2b and 2e. Afterwords, its neighbors’ discs $\Psi_{2}$ and $\Psi_{4}$ are shrunk when applying $s_{3}^{-1}<1$ on the third column of ${\mathbf{B}}$, as shown in Fig. 2c and 2f. Note that $s_{3}$ is offset by $s_{3}^{-1}$ on the main diagonal of disc $\Psi_{3}$, so its center $b_{33}=1+d_{3}$ is unchanged.

With the above two operations, we can select a sampling set and shift left-ends of the corresponding Gershgorin discs from [math] to $1$. At the same time, we scale radii of all Gershgorin discs to align their left-ends, resulting in a tighter lower bound of $\lambda_{\min}({\mathbf{B}})$. Given a graph ${\cal G}$ and a sampling budget $K$, maximizing the lower bound using the two disc operations (encoded in vectors ${\mathbf{a}}$ and ${\mathbf{s}}$) can be formulated as:

[TABLE]

where the objective is to maximize the smallest left-end $\min_{i\in\{1,\ldots,N\}}~{}c_{ii}-\sum_{j\neq i}|c_{ij}|$ of all Gershgorin discs of similarity-transformed matrix ${\mathbf{C}}$ of original ${\mathbf{A}}+\mu{\mathbf{L}}$. We call this formulated optimization (20) the primal problem.

IV-C Formulating the Dual Problem

The primal problem (20) is difficult to solve because of the combinatorial nature of $a_{i}$’s and the max-min objective. Instead, we formulate a dual problem222Though the dual problem is not exactly the same definition as classical optimization setups in linear programming [49] and convex programming [50], the essence of exchanging the roles of objective function and constraints to obtain a related optimization problem is identical. (21) that is essentially a rewriting of (20) with the objective and the constraint on $\sum_{i=1}^{N}a_{i}$ swapping roles. Specifically, given a graph ${\cal G}$ and a lower bound target $T$, (21) minimizes the number of sampled nodes subject to a constraint to lower-bound all left-ends:

[TABLE]

Instead of a max-min criterion, the objective in (21) is now linear in $a_{i}$’s.

Optimal solutions to the primal problem (20) and the dual problem (21) are related via the following proposition:

Proposition 2**.**

If $T$ is chosen such that there exists at least one optimal solution $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ to the dual problem (21) satisfying $\sum_{i}\hat{a}_{i}=K$, then there exists one optimal solution to the dual problem (21) at lower bound target $T$ that is also an optimal solution to the primal problem (20).

Proof.

For feasible solutions $({\mathbf{a}},{\mathbf{s}})$ to the dual problem (21) with $\sum_{i}a_{i}>K$, they are not feasible to the primal problem (20), and thus can be ignored. For each solution $({\mathbf{a}},{\mathbf{s}})$ to the dual problem (21) with $\sum_{i}a_{i}<K$ satisfying all constraints except $\min_{i}c_{ii}-\sum_{j\neq i}|c_{ij}|\geq T$, they must have $\min_{i}c_{ii}-\sum_{j\neq i}|c_{i,j}|<T$, otherwise $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ would not be the optimal solution with $\sum_{i}\hat{a}_{i}=K$. Thus this solution $({\mathbf{a}},{\mathbf{s}})$, while feasible in the primal, has a worse objective than $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ in the primal problem. So the optimal solution to the primal problem (20) must be a solution to the dual problem (21) with $\sum_{i}a_{i}=K$ (though not necessarily the one $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ if multiple optimal solutions exist). ∎

Denote by $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ an optimal solution to the dual problem (21) given lower bound target $T$, and ${\hat{K}}_{T}$ the corresponding objective function value. We can further prove that ${\hat{K}}_{T}$ is non-decreasing with respect to $T$.

Proposition 3**.**

Let ${\hat{K}}_{T}$ be the optimal objective function value to (21) given lower bound $T$. Then, we have ${\hat{K}}_{T_{1}}\geq{\hat{K}}_{T_{2}}$, $\forall T_{1}\geq T_{2}$.

Proof.

Since $T_{1}\geq T_{2}$, ${\hat{K}}_{T_{1}}$ is a feasible objective function value to (21) given lower bound $T_{2}$. Hence, ${\hat{K}}_{T_{1}}\geq{\hat{K}}_{T_{2}}$ must be satisfied, since ${\hat{K}}_{T_{2}}$ is the optimal (minimum) objective function value to (21) given $T_{2}$. ∎

Proposition 2 and Proposition 3 provide us with an important insight:

Instead of solving the primal problem (20) directly, we can solve the dual problem (21) iteratively for different $T\in(0,1)$, where we employ a binary search to seek the largest ${\hat{T}}$ such that the objective value of (21) is $K$ exactly.

When this largest ${\hat{T}}$ is found, we achieve an optimal solution $({\hat{\mathbf{a}}},{\hat{\mathbf{s}}})$ to (20), and ${\hat{T}}$ is the corresponding optimal objective value. When solving the dual problem (21), we select the smallest sampling set and scale radii of all Gershgorin discs, so that the discs’ left-ends (the lower bound of $\lambda_{\min}({\mathbf{B}})$) is no smaller than target $T$. Hence, we refer to the dual problem (21) also as a disc alignment problem.

When left-ends of all discs are aligned to $T$ as closely as possible, not only is the smallest left-end of Gershgorin discs a lower bound for $\lambda_{\min}({\mathbf{B}})$, as stated in GCT, but we can also prove that the largest left-end gives an upper bound:

Proposition 4**.**

Let $\Gamma_{i}$ be the Gershgorin disc from the $i$-th row of matrix ${\mathbf{C}}={\mathbf{S}}{\mathbf{B}}{\mathbf{S}}^{-1}$, where ${\mathbf{S}}$ is a non-singular diagonal matrix satisfying ${\mathbf{S}}=\operatorname*{diag}({\mathbf{s}})$ where ${\mathbf{s}}>\mathbf{0}$. Denote by $\ell_{i}$ the left-end of $\Gamma_{i}$, we have

[TABLE]

In words, while the smallest left-end $\ell_{i}$ gives a lower bound for $\lambda_{\min}({\mathbf{B}})$, the largest left-end $\ell_{i}$ gives an upper bound.

The proof of Proposition 4 is in the Appendix A. Note that $\lambda_{\min}({\mathbf{B}})$ can be exactly the lower bound $T$, if the left-ends of all Gershgorin discs are aligned at $T$ satisfying $\min_{i}\ell_{i}=\max_{i}\ell_{i}$.

V Sampling Set Selection Algorithm Development

Although, unlike the primal (20), the dual (21) has a linear objective, it is nonetheless difficult to solve optimally. Instead of manipulating all discs simultaneously, we first define a coverage subset $\Omega_{i}$ for each node $i$—the maximum subset of discs whose left-ends can be aligned at or beyond $T$ by sampling only node $i$. We introduce an approximation algorithm to pre-compute each coverage subset $\Omega_{i}$ (Algorithm 1) in Section V-A.

Pre-computed coverage subsets lead to a coarse but intuitive reinterpretation of the dual (21) as the disc coverage problem, which we prove to be NP-hard via reduction from the famous set cover (SC) problem in Section V-B. The reinterpretation enables us to derive a fast algorithm for the dual (21) based on a known SC approximation algorithm (Algorithm 2) in the combinatorial optimization literature in Section V-C.

We further employ a binary search to seek the largest ${\hat{T}}$ such that the objective of the dual (21) is $K$ exactly. The complete sampling algorithm (Algorithm 3) is proposed in Section V-D, which is summarized in Fig. 3. Finally, we discuss the computational complexities of different graph sampling approaches in Section V-E.

V-A Coverage Subset of Each Sampled Node

Before we approximately solve (21)—align all Gershgorin discs’ left-ends at or beyond target $T\in(0,1)$ using as few sampled nodes as possible—we first define a coverage subset $\Omega_{i}$ for each node $i$.

Definition 1** (Coverage subset).**

A coverage subset $\Omega_{i}$ for each node $i$ is the optimal subset ${\cal X}$ satisfying

[TABLE]

In words, a coverage subset $\Omega_{i}$ is the maximum subset of discs whose left-ends can be aligned at or beyond $T$ by sampling only node $i$.

We propose an approximation algorithm to estimate each coverage subset $\Omega_{i}$. As described in Section IV-B, disc scaling operation can balance left-ends of discs in the local neighbourhood of a sampled node. Starting from the sampled node $i$, we progressively employ disc scaling operations to align left-ends of discs at $T$ with increasing hops from the node $i$, as illustrated in Fig. 4. The algorithm stops when farther nodes cannot be aligned.

Specifically, we first sample node $i$ (thus moving the corresponding disc $\Psi_{i}$’s center $b_{ii}$ from $\mu d_{i}$ to $1+\mu d_{i}$, and $\Psi_{i}$’s left-end $a_{i}$ from [math] to $1$). We next apply scalar $s_{i}$ to expand $\Psi_{i}$’s radius $R_{i}$ and align its left-end at exactly $T$. Scalar $s_{i}$ must hence satisfy

[TABLE]

where initially $s_{j}=1,\forall j$. ${\cal N}_{i}$ denotes the neighbors of node $i$ satisfying $w_{ij}>0$. Solving for $s_{i}$ in (24), we obtain

[TABLE]

Computing $s_{i}$ for sample node $i$ using (25) implies that $s_{i}>1$ (i.e., expansion of $\Psi_{i}$’s radius), which means that node $i$’s neighbors’ discs $\Psi_{j}$’s radii will shrink due to $s_{i}^{-1}$. Specifically, disc left-end $b_{jj}-R_{j}$ of an unsampled neighbor $j$’s disc (thus $a_{j}=0$) is now:

[TABLE]

If a neighboring disc $\Psi_{j}$’s left-end remains smaller than $T$, then nothing further needs to be done. However, if $\Psi_{j}$’s left-end is larger than $T$, then $j$ belongs to $\Omega_{i}$. We then subsequently expand its radius to align its left-end at $T$ using (25). This shrinks the disc radii of node $j$’s neighbors and so on; The sequence of nodes explored starting from sample $i$ is determined using Breadth First Search (BFS) [15]. Expansion factor $s_{j}$ decreases with hops away from node $i$. Since we always expand a current disc ($s_{j}\geq 1$) leading to shrinking of neighboring discs (${s_{j}}^{-1}\leq 1$) in each step, the left-end of each already scaled node remains larger than or equal to $T$.

Note that $1\leq|\Omega_{i}|\leq N$ for any value of $T\in(0,1)$. When $T\rightarrow 0$, $|\Omega_{i}|\rightarrow N$, and we need to visit all $N$ nodes using BFS. The worst-case time complexity for computing a single $\Omega_{i}$ is ${\mathcal{O}}(N+M)$, where $M=|{\cal E}|$. Typically, sampling budget $K$ is much larger than $1$, and it is not necessary to compute $\Omega_{i}$ corresponding to a very small $T$. To enable faster computation, we set a parameter $p$ to estimate $\Omega_{i}$ only within a $p$-hop neighborhood from node $i$. Thus, the time complexity of estimating $\Omega_{i}$ within $p$ hops is ${\mathcal{O}}(P+Q)$, where $P$ is the number of nodes within $p$ hops from node $i$, and $Q$ is the number of edges within $p$ hops from node $i$. Assuming a maximum degree per node $d_{\max}$, we have $P\leq Q\leq d^{p}_{\max}$. $d^{p}_{\max}$ is achieved only if graph ${\cal G}$ is a tree. Since $p$ is typically set much smaller than the number of hops required to traverse the entire graph with BFS, $P\ll N$ and $Q\ll M$. The procedure of estimating coverage subset $\Omega_{i}$ is sketched in Algorithm 1.

Illustrative Example on 5-node Path Graph: We use an example to illustrate the operations of the above algorithm, shown in Fig. 5. We assume the same $5$-node path graph in Fig. 1 and choose to sample node $3$. Assuming $\mu=1$, the graph’s coefficient matrix ${\mathbf{B}}$ with entry $(3,3)$ updated is shown in Fig. 2a. The left-end of node 3’s Gershgorin disc shifts from [math] to $1$, as shown in Fig. 2d.

We next perform disc scaling to align the left-end of disc $\Psi_{3}$ to $T$. As shown in Fig. 5a, scalar $s_{3}$ is applied to the third row of ${\mathbf{B}}$, and thus the radius of disc $\Psi_{3}$ is expanded by $s_{3}$, where $s_{3}>1$. Simultaneously, scalar $s_{3}^{-1}$ is applied to the third column, and thus the radii of discs $\Psi_{2}$ and $\Psi_{4}$ are shrunk due to the scaling of $w_{23}$ and $w_{43}$ respectively by $s_{3}^{-1}$. Note that entry $(3,3)$ of ${\mathbf{B}}$ (and $\Psi_{3}$’s disc center) is unchanged, since scalar $s_{3}$ is offset by $s_{3}^{-1}$. We see that due to the expansion of disc $\Psi_{3}$ radius, the disc left-ends of its neighboring nodes (nodes $2$ and $4$) move beyond threshold $T$, as shown in Fig. 5d.

We next apply scalar $s_{2}$ to disc $\Psi_{2}$ to expand its radius by $s_{2}$, where $s_{3}>s_{2}>1$, and the radii of discs $\Psi_{1}$ and $\Psi_{3}$ are shrunk due to the scaling of $w_{12}$ and $w_{32}$ by $s_{2}^{-1}$, as shown in Fig. 5b. $s_{2}$ must be smaller than $s_{3}$ for the left-end of $\Psi_{2}$ not to move past [math]. The discs are shown in Fig. 5e. Scaling $\Psi_{2}$ does not move the left-end of disc $\Psi_{1}$ beyond threshold $T$. Thus, we stop aligning the further-away disc $\Psi_{1}$.

Subsequently, similar disc operations can be performed on $\Psi_{4}$ as shown in Fig. 5c and 5f. Scaling $\Psi_{4}$ does not move the left-end of disc $\Psi_{5}$ beyond $T$. Thus, we stop aligning the further-away disc $\Psi_{5}$. Finally, the left-ends of discs $\Psi_{2}$, $\Psi_{3}$ and $\Psi_{4}$ are moved beyond threshold $T$, i.e., $\Omega_{3}=\{2,3,4\}$.

Analogy of Throwing Pebbles in a Pond: To impart intuition, we consider the following real-world analogy of computing subset $\Omega_{i}$ for sample $i$. Suppose a planar graph $\mathcal{G}$ models the water surface of a pond, where nodes denote water surface locations. Sampling a node $i$ is analogous to throwing a pebble at a water surface location $i$ in the pond. The thrown pebble causes a ripple to adjacent locations $j\in\mathcal{N}_{i}$, where the magnitude of the ripple at location $j$ is smaller than at $i$. Analogously, when aligning left-ends of Gershgorin discs at $T$, scalar $s_{j}$ of neighbor $j\in\mathcal{N}_{i}$ is necessarily smaller than $s_{i}$. The ripple in the pond dissipates away from the thrown location $i$. Analogously, scalars $s_{j}$ becomes smaller further away in hop count from sample $i$. The ripple stops increasing in size when the dissipating fringes become unnoticeable. Analogously, subset $\Omega_{i}$ stops growing when left-ends of neighboring discs fail to move beyond target $T$.

V-B NP-Hardness of Disc Coverage Problem

After pre-computing coverage subsets, we can reinterpret the dual (21) more coarsely but in a more intuitive manner: given pre-computed coverage subsets $\Omega_{i},\forall i\in\{1,\ldots,N\}$, how to select minimum-cost sampling vector ${\mathbf{a}}$ so that every node $j$ belongs to at least one sample’s coverage subset $\Omega_{i}$, i.e.,

[TABLE]

We term this the disc coverage (DC) problem. DC (27) is related to the dual (21) in the sense that: i) they have the same optimization objective, and ii) every feasible solution ${\mathbf{a}}$ to (27) is also a feasible solution to (21), where the scalars ${\mathbf{s}}$ corresponding to ${\mathbf{a}}$ are the ones calculated when coverage subsets $\Omega_{i}$’s are pre-computed. Thus solving DC provides an approximate feasible solution to the dual.

Unfortunately, even this DC problem (27) is NP-hard. A simple NP-hardness proof via a reduction from the famous NP-hard set cover (SC) problem [15] is provided in Appendix B. However, the resemblance of the DC problem to SC enables us to derive an efficient approximation algorithm to the dual (21), which we describe next.

V-C Disc Alignment via Greedy Set Cover

Given nodes ${\cal U}=\{1,2,\ldots,N\}$ and a collection of pre-computed $N$ coverage subsets $\Omega=\{\Omega_{1},\ldots,\Omega_{N}\}$, we derive a fast algorithm to DC (27) based on an error-bounded approximation algorithm to SC. We first describe SC formally.

Definition 2** (Set cover problem).**

An instance $({\cal X},{\cal F})$ of the set cover problem consists of a finite element set ${\cal X}$ and a family ${\cal F}$ of subsets of ${\cal X}$, such that every element $x\in{\cal X}$ belongs to at least one subset ${\cal S}\in{\cal F}$:

[TABLE]

We say that a subset ${\cal S}\in{\cal F}$ covers its elements. The problem is to find a minimum-size subset collection ${\cal C}\subseteq{\cal F}$ whose members cover all of ${\cal X}$:

[TABLE]

Given the resemblance between DC and SC, we borrow a known greedy approximation algorithm for SC with error-bounded performance [15]: at each iteration, we select a yet-to-be-chosen subset ${\cal S}$ that covers the largest number of uncovered elements in ${\cal X}$. Similarly, for our algorithm for DC, we select sample $i$ and coverage subset $\Omega_{i}$ with the largest uncovered nodes at each iteration. When a sampling vector ${\mathbf{a}}$ satisfying the condition that each node belongs to at least one sample coverage subset $\Omega_{i}$ is computed, the greedy algorithm can achieve a $H(\max_{i}|\Omega_{i}|)$-approximate solution to DC [15], i.e.,

[TABLE]

where $H(n)=\sum_{i=1}^{n}1/i$ is the harmonic series and $\hat{{\mathbf{a}}}$ is the optimal solution to DC. $|{\mathbf{a}}|$ denotes $\sum_{i}a_{i}$.

Since we are given a sampling budget $K$, we terminate the algorithm when the number of selected samples reaches $K$ for practical implementation. The resulting sampling vector ${\mathbf{a}}$ is valid only if $|{\mathbf{a}}|\leq K$ and each node belongs to at least one sample coverage subset $\Omega_{i}$. The procedure of disc alignment via greedy set-covering is outlined in Algorithm 2.

The time complexity of $N$ subsets preparation within $p$ hops is ${\mathcal{O}}(PN+QN)$. We represent uncovered elements in ${\cal U}$ using [math]-$1$ bit array. The uncovered elements in $\Omega_{i}$ to ${\cal U}$ can be computed through at most ${\mathcal{O}}(P)$ times fast bit operations, since the number of “1” in the bit array of $\Omega_{i}$ is at most $P$. We compute uncovered elements of $N$ subsets at each subset selection stage. We employ no more than $K$ times subset selections, thus the time complexity of the greedy subset selection algorithm is ${\mathcal{O}}(NKP)$. The total time complexity of Algorithm 2 comprises ${\mathcal{O}}\left(PN+QN\right)$ node search plus ${\mathcal{O}}\left(NKP\right)$ bit operations.

V-D Binary Searching the Maximum Lower Bound

According to Proposition 3, the number of sampled nodes $K$ is non-decreasing with respect to lower bound $T$ in the dual problem (21). Hence, we can perform binary search to find the maximum lower bound ${\hat{T}}\in(0,1)$ corresponding to $K$ sampled nodes. We name the proposed algorithm Binary Search with Gershgorin Disc Alignment (BS-GDA), as outlined in Algorithm 3.

At each iteration, we estimate validity flag $VF$ and sampling set ${\cal S}$ given a certain lower bound $T$ using Algorithm 2. If $VF=false$, then the lower bound $T$ is set too large, and we update $right$ to reduce $T$. Otherwise, $VF=true$, then $T$ may be still too small, leading to $|{\cal S}|<K$ or not the maximum $T$ corresponding to $|{\cal S}|=K$. We update $left$ to increase $T$. When $right-left\leq\epsilon$, BS-GDA converges and we find the maximum lower bound ${\hat{T}}$ with numerical error lower than $\epsilon$.

In order to achieve numerical precision $\epsilon$ in BS-GDA, we need to employ ${\mathcal{O}}(\log_{2}\frac{1}{\epsilon})$ times binary search. For instance, to achieve numerical precision $\epsilon=10^{-3}$, we employ $10$ times binary search. The total time complexity of Algorithm 3 comprises ${\mathcal{O}}\left(N(P+Q)\log_{2}\frac{1}{\epsilon}\right)$ node search plus ${\mathcal{O}}\left(NKP\log_{2}\frac{1}{\epsilon}\right)$ bit operations.

V-E Complexity Analysis

We here summarize the theoretical computational complexities of different graph sampling approaches, and list them in Table I, assuming that the signal bandwidth $\omega$ is the $K$-th eigenvalue $\theta{{}_{K}}$ of ${\mathbf{L}}$ and the graph is sparse. We separate the complexity into two parts: preparation (for eigen-decomposition or computing initial preparation information) and sampling (during greedy sampling step). The competing algorithms are referred to as E-optimal [24], Spectral proxies (SP) [16], MFN [22], MIA [20] and Eigendecomposition-free (Ed-free) [9], respectively. Complexities of methods E-optimal, SP, MFN and MIA are directly borrowed from paper [20], with modified notations. For E-optimal and MFN, $T_{1}$ is the convergence steps for computing the first $K$ eigen-pairs of Laplacian ${\mathbf{L}}$ and $R$ is a constant for combining complexities. In strategy SP, $k$ is the degree of spectral proxies and $T_{2}$ is the convergence steps for obtaining the first eigen-pair of ${\mathbf{L}}$ via LOBPCG [16]. The parameters $q$ and $F$ in MIA are the degree of Chebyshev polynomials for a low-pass filter approximation and the truncated parameter of a Neuamnn series, respectively. As for Ed-free method, paper [9] claimed that its complexity was dominated by the non-zero entries in matrix ${\mathbf{L}}^{\bar{q}}$, i.e, $J$ in paper [9], where $\bar{q}$ is degree of Chebyshev polynomial approximation. Similarly, the non-zero entries in ${\mathbf{L}}^{\bar{q}}$ is at most $Nd_{\max}^{\bar{q}}$ in a degree-bounded graph, so we use $\bar{Q}$ to denote the largest number of edges in $\bar{q}$-hop neighborhood, which implies $J={\mathcal{O}}(N\bar{Q})$. From Table I, we can see that the complexities of Ed-free and our proposed BS-GDA method are both roughly linear in terms of graph size $N$, which is substantially lower than other eigen-decomposition-based sampling methods and implies those two methods are applicable in fairly large graphs. The explicit execution time comparison will be demonstrated in the coming experimental section.

VI Experiments

VI-A Experimental Setup

We demonstrate the efficacy of our proposed fast graph subset sampling algorithm (BS-GDA) via extensive simulations. Our experimental platform is Windows 10 desktop computer with Intel i5-4670K CPU and 24GB memory. All algorithms are performed on MATLAB R2016b.

We use four types of graphs in GSPBOX [51] for testing:

Random sensor graph with $N$ nodes. 2. 2.

Community graph with $N$ nodes and $\lfloor\sqrt{N}/2\rfloor$ random communities. 3. 3.

Barabási-Albert graph with $N$ nodes. Each graph is constructed by adding new nodes, each with $m=1$ edge, that are preferentially attached to existing nodes with high degrees. 4. 4.

Minnesota road graph with fixed $N=2642$ nodes.

Random sensor graphs are weighted sparse graphs generated by GSPBOX, where each node connects to its six nearest neighbours. Edge weights are computed using:

[TABLE]

where ${\mathbf{x}}_{i}$ is the 2D coordinate of node $i$, and $\sigma_{x}$ is set automatically by the toolbox. Community graphs, Barabási-Albert graphs and Minnesota road graph are unweighted sparse graphs generated by GSPBOX, where we manually compute weight for each edge. For community graphs and Minnesota road graph, we compute edge weights using (31). We set parameter $\sigma_{x}=1$ for community graphs and set parameter $\sigma_{x}=0.1$ for Minnesota road graphs. Since Barabási-Albert graphs are scale-free graphs without 2D coordinates, we randomly generate edge weights using a uniform probability distribution with interval $(0,1)$. Illustrations of four types of graphs are shown in Fig. 6.

For each graph type, we consider two types of graph signals described as follows:

GS1: The true signals are exactly $\omega$-bandlimited, where $\omega=\theta_{\lfloor N/10\rfloor}$ is the $\lfloor N/10\rfloor$-th eigenvalue of ${\mathbf{L}}$. For example, $\omega=\theta_{50}$ for $N=500$. The non-zero GFT coefficients are randomly generated from ${\cal N}(0,10)$. We then add i.i.d. Gaussian noise generated from ${\cal N}(0,0.1^{2})$ such that the SNR equals to $20$ dB. 2. 2.

GS2: The true signals are generated from multivariate Gaussian distribution ${\cal N}\left(\mathbf{0},({\mathbf{L}}+\delta{\mathbf{I}})^{-1}\right)$, where $\delta=10^{-5}$. Because the power of the generated graph signals is inconsistent, we normalize the signals using ${\mathbf{x}}^{\prime}=\frac{{\mathbf{x}}-mean({\mathbf{x}})}{std({\mathbf{x}})}$, where $mean({\mathbf{x}})=\sum_{i}x_{i}/N$ and $std({\mathbf{x}})=\sqrt{\frac{\sum_{i}\left(x_{i}-mean({\mathbf{x}})\right)^{2}}{N}}$, respectively. We then add i.i.d. Gaussian noise generated from ${\cal N}(0,0.1^{2})$ such that the SNR equals to $20$ dB.

We compare the performance of our proposed algorithm with six recent graph sampling set selection methods, of which implementations are available:

Random graph sampling with non-uniform probability distribution (referred to as Random) [27]. 2. 2.

Deterministic sampling set selection methods, including E-optimal [24], Spectral proxies (SP) [16], MFN [22], MIA [20] and Eigendecomposition-free (Ed-free) [9], respectively.

Our proposed algorithm BS-GDA avoids expensive eigen-decomposition and matrix computations, thus we build MEX functions for the iterative BFS and the greedy set covering, leading to a MATLAB+MEX implementation. For all experiments, we set the target numerical precision $\epsilon=10^{-5}$ for binary search and $\mu=0.01$. We empirically set hop constraint $p=12$ to consider no more than $12$ hops for each subset $\Omega_{i}$ estimation. It is analogus to the neighborhood considered by Chebyshev polynomial approximation with order $12$ in [9].

VI-B Running Time Comparisons

We compare the running time of our proposed BS-GDA with other competing algorithms on random sensor graphs and community graphs with increasing graph size $N$. We select sampling budget $K=N/10$ nodes for both graphs. Fig. 7 shows the curves of running time versus graph size $N$. Specifically, for $N=3000$, we further compute the speedup factor (SF) [9] of our algorithm with respect to the competing algorithms using

[TABLE]

as summarized in Table. II.

We observe that BS-GDA is the fastest deterministic sampling set selection algorithm for both graphs. As mentioned, E-optimal, SP and MFN all require computation of extreme eigen-pairs of the graph Laplacian matrix (or submatrix), which is computationally demanding. Though there is no explicit computation of eigen-pairs in MIA, the large matrix series computation also results in large computation cost. Thus, BS-GDA is hundreds to thousands times faster than E-optimal, SP, MFN and MIA.

The recent eigen-decomposition-free method, Ed-free [9], utilizes the localized operator to design an intuitive sampling strategy, whose complexity is ${\mathcal{O}}(NK\bar{Q})$ in Table I. Our BS-GDA is faster than Ed-free for both graphs, in particular on the community graph. The reason is that edges within each community are relatively dense, and thus $\bar{Q}$ for community graph is much larger than $\bar{Q}$ for random sensor graph, leading to a slower speed for Ed-free. In contrast, the computation efficiency of BS-GDA is similar on both graphs.

Random [27] is faster than our proposed BS-GDA, but its reconstruction error performance is noticeably worse than BS-GDA, as shown in the following section.

VI-C Reconstruction Error Comparisons

We also compare the resulting reconstruction MSEs between our proposed algorithm and the competing ones. We choose $N=500$ for random sensor graphs, community graphs and Barabási-Albert graphs. Minnesota road graph is of fixed size $N=2642$. All the graph sampling set selection algorithms are compared on these four graphs with increasing sampling budgets.

For each graph, we consider two types of graph signals, as described in Sec. VI-A. We randomly generate $50$ graph signals for each type, and for each signal we add $50$ different i.i.d. Gaussian noise, resulting in $2500$ noisy signals for each experiment. For signal type GS2, because the generated signals are not strictly bandlimited, we approximately set the bandwidth to $\omega=\theta_{\lfloor N/10\rfloor}$ for sampling algorithms that require this parameter.

We use our proposed BS-GDA and the competing algorithms to sample the noisy signals. All the competing algorithms are run with the default or the recommended settings in the referred publications. We then reconstruct original graph signals by solving GLR-based reconstruction (8) with $\mu=0.01$.

The average MSEs in terms of sampling budgets are illustrated in Fig. 8. As shown, our proposed BS-GDA achieves comparable or smaller MSE values to the competing deterministic schemes on all the four graphs and achieves much better performance than the random sampling scheme Random.

For eigen-decomposition based algorithms, the performances of E-optimal, MFN and SP are unsatisfactory when the sampling budget is much smaller than the bandwidth. Further, these eigen-decomposition based algorithms are computationally too expensive for large graphs, such as the Minnesota road graph. For eigen-decomposition-free algorithms, MIA and Ed-free work well on the random sensor graphs, the community graphs and the Barabási-Albert graphs, but their performance is poor on the large Minnesota road graph. To avoid explicit eigen-decomposition, symmetrically normalized ${\mathbf{L}}_{n}$ is designated for MIA and Ed-free, whose eigenvalues are theoretically bounded in $[0,2]$, to utilize Chebyshev polynomials of ${\mathbf{L}}_{n}$ to approximate the low-pass graph filter operator. These two schemes are not robust enough on the graph signals generated on the Minnesota road graph. It has been observed that, in general, the performance of Random is not comparable to deterministic graph sampling algorithms.

Since the previous deterministic graph sampling algorithms such as E-optimal and MFN are designed assuming a standard least-squares reconstruction, we further test all the sampling algorithms with least-squares reconstruction on random sensor graphs, as shown in Fig. 9. We set the sampling budget to be no smaller than the bandwidth for GS1 or GS2. As observed in Fig. 9, E-optimal, MFN and SP achieve robust performance. For MIA, Ed-free and the proposed BS-GDA, their reconstruction MSEs are somewhat larger when the sampling budget is close to the bandwidth. As the sampling budget increases, all deterministic sampling algorithms have similar reconstruction MSEs. In contrast, the performance of Random is noticeably worse than the deterministic sampling algorithms.

The proposed BS-GDA solves the dual (21) via greedy set-covering and performs a binary search to find the optimal lower bound of the primal (20), which reduces the global upper bound of the MSE value. Hence, BS-GDA guarantees robust performance for different graphs in the experiments. We further show an illustrative example of sampling $K=11$ nodes on the community graph of size $N=500$ with 11 communities in Fig. 10. This illustration clearly shows that BS-GDA captures the graph structure very well and evenly samples one node in each community. This matches our sampling analogy to the set cover (SC) problem, where each sampled node covers many uncovered nodes.

VII Conclusion

To address the “graph sampling with noise” problem, in this paper, we first assume a biased graph Laplacian regularization (GLR) based signal reconstruction scheme from samples, and then propose a fast graph sampling set selection algorithm to maximize the stability of the resulting linear system. In particular, our proposed BS-GDA maximizes the lower-bound of $\lambda_{\min}$ of the coefficient matrix via Gershgorin disc alignment, which is proven to mean also minimizing the upper-bound of reconstruction MSE. Unlike most previous graph sampling schemes, our proposed BS-GDA does not require computation of any eigen-pairs. Extensive simulation experiments show that our disc-based sampling algorithm runs substantially faster than existing eigen-decomposition-based schemes (hundreds to thousands times faster) for medium-sized graphs and outperforms other eigen-decomposition-free schemes in reconstruction error for large-sized graphs.

Appendix A Proof of Proposition 4

Proof.

Since matrix ${\mathbf{C}}$ is similar to ${\mathbf{B}}$, we have $\lambda_{\min}({\mathbf{B}})=\lambda_{\min}({\mathbf{C}})$. Hence (22) is equivalent to

[TABLE]

From GCT, we see that $\lambda_{\min}({\mathbf{C}})$ is bounded from below by the smallest left-end $\ell_{i}$, i.e.

[TABLE]

We next prove that the largest left-end $\ell_{i}$ is the upper bound of $\lambda_{\min}({\mathbf{C}})$, i.e.,

[TABLE]

Because ${\mathbf{B}}={\mathbf{A}}+\mu{\mathbf{L}}$ and ${\mathbf{S}}=\operatorname*{diag}({\mathbf{s}})$, ${\mathbf{s}}>\mathbf{0}$, the diagonal elements of ${\mathbf{C}}$ are non-negative $c_{ii}\geq 0$ and the off-diagonal elements of ${\mathbf{C}}$ are non-positive $c_{ij}\leq 0$ for $i\neq j$. We consider the matrix $\widetilde{\mathbf{C}}:=-{\mathbf{C}}+t{\mathbf{I}}$, where $t$ is chosen large enough so that the diagonal elements of $\widetilde{\mathbf{C}}$ are all non-negative, and hence the whole matrix $\widetilde{\mathbf{C}}$ is non-negative. Then by the Perron-Frobenius Theorem [35], $\widetilde{\mathbf{C}}$ has a positive eigenvalue $\phi_{\max}$ (with a non-negative eigenvector), which is also the eigenvalue of largest absolute value. Since the eigenvalues $\lambda$’s of ${\mathbf{C}}$ and those $\phi$’s of $\widetilde{\mathbf{C}}$ are related simply by $\phi=-\lambda+t$, it follows that

[TABLE]

is the minimum eigenvalue of ${\mathbf{C}}$.

Let ${\mathbf{v}}=[1,1,\ldots,1]^{\top}$ and consider the vector (recalling that $\widetilde{c}_{ij}=-c_{ij}=|c_{ij}|$ for $i\neq j$)

[TABLE]

Now by Collatz-Wielandt’ bound [35], we have

[TABLE]

Together with (36), we have

[TABLE]

Thus, we conclude that

[TABLE]

as required. ∎

Appendix B Proof of NP-Completeness

We prove that the decision version of the disc coverage (DC) problem (27) is NP-complete by transforming the decision version of the set cover (SC) problem [15] to a special case of DC. The SC decision problem asks the following binary decision question: given a family of subsets ${\cal F}$ of a finite element set ${\cal X}$ and positive integer $K\leq|{\cal F}|$, does there exist a cover ${\cal C}$ of size $K$ or less, i.e., ${\cal C}\subseteq{\cal F}$ and $|{\cal C}|\leq K$, where each element $x\in{\cal X}$ is in at least one member of ${\cal C}$?

For each instance of SC, we perform the following transformation. First, for each element $j$ in set, $x_{j}\in{\cal X}$, we construct one element node $j$ in a graph ${\cal G}$. Next, for each subset in collection $c_{i}\in{\cal F}$, we construct one subset node $|{\cal X}|+i$ and connect it to element nodes $j$ that correspond to elements $x_{j}\in c_{i}$. We construct one super node $|{\cal X}|+|{\cal F}|+1$ and connect it to $|{\cal F}|$ subset nodes. Finally, we construct $K+1$ extra nodes and connect them to the super node. Coverage subset $\Omega_{i}$ for each node $i$ is set to be all 1-hop neighbors of node $i$ in the constructed graph.

As an example, consider the following SC instance:

[TABLE]

Fig. 11 illustrates the corresponding constructed graph. The corresponding DC decision problem is: does there exist a sampling vector ${\mathbf{a}}$ with objective of (27) $K+1$ or less, while having each graph node $i$ included in at least one sample $i$’s coverage subset $\Omega_{i}$?

We prove that the SC decision problem is true iff the corresponding DC decision problem is true. Suppose first that there exists a cover ${\cal C}$ of size $|{\cal C}|\leq K$. We observe that a corresponding sample selection of super node $|{\cal X}|+|{\cal F}|+1$ and subset nodes representing members $c_{i}\in{\cal C}$ will mean that each graph node is included in at least one sample coverage subset $\Omega_{i}$, i.e., all nodes are within one hop of selected samples.

Now suppose that there exists a sample selection of size $K+1$ such that all nodes are within one hop of selected samples. We first observe that the super node must be selected. Otherwise, just covering the $K+1$ extra nodes would require $K+1$ node selections, leaving no leftover budget to select and cover element and subset nodes.

Given that the super node is selected, there remain $K$ sample budget to cover $|{\cal X}|$ element nodes. We can identify a cover of $K^{\prime}$ subsets in SC, $K^{\prime}\leq K$, from the chosen $K$ node samples as follows. Suppose there are $K_{1}$, $K_{2}$ and $K_{3}$ selected element, subset and extra nodes respectively, where $K_{1}+K_{2}+K_{3}=K$. For each one of $K_{2}$ chosen subset node, we select the corresponding subset in SC. For each one of $K_{1}$ chosen element node, we randomly select one of its connected subset node and choose the corresponding subset in SC (if not already chosen in previous step). We ignore all $K_{3}$ chosen extra nodes. We see that the number of chosen subsets in SC is no larger than $K_{1}+K_{2}\leq K$, but covers all elements in set ${\cal X}$. Thus a feasible sample set for the DC decision problem maps to a feasible cover ${\cal C}$ for the SC decision problem. We have thus proven that an instance of the SC problem maps to our construction of the DC problem. Therefore, we can conclude that the DC problem is at least as hard as the SC problem. Since the SC decision problem is known to be NP-complete, the DC decision problem is also NP-complete.

Given the DC decision problem is NP-complete, the corresponding DC optimization problem in (27) is NP-hard. $\Box$