Efficient Sampling Set Selection for Bandlimited Graph Signals Using Graph Spectral Proxies

I Introduction

Graphs provide a natural representation for data in many applications, such as social networks, web information analysis, sensor networks and machine learning . They can also be used to represent conventional data, such as images and videos . A graph signal is a function defined over the nodes of a graph. Graph signal processing aims to extend the well-developed tools for analysis of conventional signals to signals on graphs while exploiting the underlying connectivity information . In this paper, we extend the theory of sampling for graph signals by developing fast and efficient algorithms for sampling set selection.

Sampling theory of graph signals similarly deals with the problem of recovering a signal from its samples on a subset of nodes of the graph. The smoothness assumption on a graph signal is formalized in terms of bandlimitedness in a graph Fourier basis. The graph Fourier basis is given by the eigenvectors and eigenvalues of certain variation operators (e.g., graph Laplacian) that measure the variation in a graph signal while taking into account the underlying connectivity. To formulate a sampling theory for graph signals we need to consider the following questions: 1. Given a subset of nodes to be sampled, what is the maximum bandwidth (in the appropriate graph Fourier domain) that a signal can have so that it can be uniquely and stably reconstructed from those samples? 2. Given the signal bandwidth, what is the best subset of nodes to be sampled for a unique and stable reconstruction? Stability is an important issue in the choice of sampling set. In practice, signals are only approximately bandlimited and/or samples are noisy. A poor choice of sampling set can result in a very ill-conditioned reconstruction operator which amplifies the sample perturbations caused by noise and model mismatch and thus, lead to large reconstruction errors. Hence, selecting a sampling set that gives stable reconstructions is vital.

The problem of selecting sampling sets for recovery of smooth, bandlimited graph signals, arises in many applications. A prominent example is active semi-supervised learning , where a learner is allowed to specify a set of points to be labeled, given a budget, before predicting the unknown labels. In this setting, class indicator vectors can be considered as smooth or bandlimited graph signals and the set of points to be labeled as the sampling set. Therefore the task of actively choosing the training set in this scenario is equivalent to finding the best possible sampling set, under a given sampling budget. Other applications of sampling set selection include selective activation of sensors in sensor networks, and design of graph based lifting transforms in image compression . Signals of interest in these applications are also smooth with respect to the graph and the goal is to find the best sampling locations that minimize the reconstruction error.

Most recent approaches for formulating a sampling theory for graph signals involve two steps – first, computing a portion of the graph Fourier basis, and second, using the basis elements either to check if a unique and stable reconstruction is possible with the given samples or to search for the best subset for sampling. However, when the graphs of interest are large, computing and storing multiple eigenvectors of their variation operators increases the numerical complexity and memory requirement significantly. Therefore, we propose a technique that achieves comparable results by using the variation operator directly and skipping the intermediate step of eigen-decomposition.

Sampling theory for graph signals was first studied in , where a sufficient condition for unique recovery of signals is stated for a given sampling set. Using this condition, gives a bound on the maximum bandwidth that a signal can have, so that it can be uniquely reconstructed from its samples on a given subset of nodes. We refine this bound in our previous work by considering a necessary and sufficient condition for sampling. Using this condition, we also propose a direct sampling set selection method that finds a set approximately maximizing this bound so that a larger space of graph signals can be uniquely reconstructed. However, does not explain why maximizing this bound leads to stable reconstructions. Moreover, the results are specific to undirected graphs. Thus, the main contributions of this paper are to extend our prior work and propose an efficient sampling set selection algorithm that generalizes easily for different graphs, while also considering the issue of stability.

Previous methods for sampling set selection in graphs can be classified into two types, namely spectral-domain methods and vertex-domain methods, which are summarized below.

Most of the recent work on sampling theory of graph signals assumes that a portion of the graph Fourier basis is explicitly known. We classify these methods as spectral-domain approaches since they involve computing the spectrum of the variation operator. For example, the work of requires computation and processing of the first $r$ eigenvectors of the graph Laplacian to construct a sampling set that guarantees unique (but not necessarily stable) reconstruction for a signal spanned by those eigenvectors. Similarly, a greedy algorithm for selecting stable sampling sets for a given bandlimited space is proposed in . It considers a spectral-domain criterion, using minimum singular values of submatrices of the graph Fourier transform matrix, to minimize the effect of sample noise in the worst case. The work of creates a link between the uncertainty principle for graph signals and sampling theory to arrive at similar criteria in the presence of sample noise. It is also possible to generalize this approach using ideas from the theory of optimal experiment design and define other spectral-domain optimality criteria for selecting sampling sets that minimize different measures of reconstruction error when the samples are noisy (for example, the mean squared error). Greedy algorithms can then be used to find sets which are approximately optimal with respect to these criteria.

I-A2 Vertex-domain approaches

There exist alternative approaches to sampling set selection that do not consider graph spectral information and instead rely on vertex-domain characteristics. Examples include and , which select sampling sets based on maximum graph cuts and spanning trees, respectively. However, these methods are better suited for designing downsampling operators required in bipartite graph multiresolution transforms . Specifically, they do not consider the issue of optimality of sampling sets in terms of quality of bandlimited reconstruction. Further, it can be shown that the maximum graph-cut based sampling set selection criterion is closely related to a special case of our proposed approach. There exists an alternate vertex-domain sampling approach, described in the work of , that involves successively shifting a signal using the adjacency matrix and aggregating the values of these signals on a given node. However, sampling using this strategy requires aggregating the sample values for a neighborhood size equal to the dimension of the bandlimited space, which can cover a large portion of the graph.

The sampling strategies described so far involve deterministic methods of approximating optimal sampling sets. There also exists a randomized sampling strategy that guarantees a bound on the worst case reconstruction error in the presence of noise by sampling nodes independently based on a judiciously designed distribution over the nodes. However, one needs to sample much more nodes than the dimension of the bandlimited space to achieve the error bound.

I-B Contributions of this work

It is possible to extract and process useful spectral information about a graph signal even when the graph Fourier basis is not known. For example, spectral filters in the form of polynomials of the variation operator are used in the design of wavelet filterbanks for graph signals to offer a trade-off between frequency-domain and vertex-domain localization. In our work, we use a similar technique of extracting spectral information from signals using $k$-hop localized operations, without explicitly computing the graph Fourier basis elements. Our main contributions can be summarized as follows:

Motivated by spectral filters localized in the vertex domain, we define graph spectral proxies based on powers of the variation operator to approximate the bandwidth of graph signals. These proxies can be computed using localized operations in a distributed fashion with minimal storage cost, thus forming the key ingredient of our approach. These proxies have a tunable parameter $k$ (equal to the number of hops), that provides a trade-off between accuracy of the approximation versus the cost of computation.

Using these proxies, we give an approximate bound on the maximum bandwidth of graph signals (cutoff frequency) that guarantees unique reconstruction with the given samples. We show that this bound also gives us a measure of reconstruction stability for a given sampling set.

We finally introduce a greedy, iterative gradient based algorithm that aims to maximize the bound, in order to select an approximately optimal sampling set of given size.

A specific case of these spectral proxies based on the undirected graph Laplacian, and a sampling set selection algorithm, has been introduced earlier in our previous work . With respect to the key new contributions are as follows. We generalize the framework to a variety of variation operators, thereby making it applicable for both undirected and directed graphs. We provide a novel interpretation for the cutoff frequency function in as a stability measure for a given sampling set that one can maximize. We show that the spectral proxies arise naturally in the expression for the bound on the reconstruction error when the samples are noisy or the signals are only approximately bandlimited. Thus, an optimal sampling set in our formulation minimizes this error bound. We also show that our algorithm is equivalent to performing Gaussian elimination on the graph Fourier transform matrix in a certain limiting sense, and is therefore closely related to spectral-domain approaches. Numerical complexity of the proposed algorithm is evaluated and compared to other state of the art methods that were introduced after was published. Finally, we evaluate the performance and complexity of the proposed algorithm through extensive experiments using different graphs and signal models.

The rest of the paper is organized as follows. Section II defines the notation used in the paper. This is followed by the concepts of frequency and bandlimitedness for signals, on both undirected and directed graphs, based on different variation operators. In Section III, we consider the problems of bandlimited reconstruction, uniqueness and stable sampling set selection, assuming that the graph Fourier basis is known. Section IV addresses these problems using graph spectral proxies. The effectiveness of our approach is demonstrated in Section V through numerical experiments. We conclude in Section VI with some directions for future work.

II Background

II-B Notions of Frequency for Graph Signals

In order to formulate a sampling theory for graph signals, we need a notion of frequency that enables us to characterize the level of smoothness of the graph signal with respect to the graph. The key idea, which is used in practice, is to define analogs of operators such as shift or variation from traditional signal processing, that allow one to transform a signal or measure its properties while taking into account the underlying connectivity over the graph. Let ${\bf L}$ be such an operator in the form of an $N\times N$ matrixAlthough ${\bf L}$ has been extensively used to denote the combinatorial Laplacian in graph theory, we overload this notation to make the point that any such variation operator can be defined to characterize signals of interest in the application at hand.. A variation operator creates a notion of smoothness for graph signals through its spectrum. Specifically, assume that ${\bf L}$ has eigenvalues $|\lambda_{1}|\leq\ldots\leq|\lambda_{N}|$ and corresponding eigenvectors $\{{\bf u}^{1},\ldots,{\bf u}^{N}\}$. Then, these eigenvectors provide a Fourier-like basis for graph signals with the frequencies given by the corresponding eigenvalues. For each ${\bf L}$, one can also define a variation functional $\text{Var}({\bf L},{\bf f})$ that measures the variation in any signal ${\bf f}$ with respect to ${\bf L}$. Such a definition should induce an ordering of the eigenvectors which is consistent with the ordering of eigenvalues. More formally, if $|\lambda_{i}|\leq|\lambda_{j}|$, then $\text{Var}({\bf L},{\bf u}^{i})\leq\text{Var}({\bf L},{\bf u}^{j})$.

where ${\cal R}=\{1,\ldots,r\}$. The space of $\omega$-bandlimited signals is called Paley-Wiener space and is denoted by $PW_{\omega}(G)$ . Note that $PW_{\omega}(G)=\text{range}({\bf U}_{{\cal V}{\cal R}})$ (i.e., the span of columns of ${\bf U}_{{\cal V}{\cal R}}$). Bandwidth of a signal ${\bf f}$ is defined as the largest among absolute values of eigenvalues corresponding to non-zero GFT coefficients of ${\bf f}$, i.e.,

A key ingredient in our theory is an approximation of the bandwidth of a signal using powers of the variation operator ${\bf L}$, as explained in Section IV. Since this approximation holds for any variation operator, the proposed theory remains valid for all of the choices of GFT in Table I.

II-C Examples of variation operators

In undirected graphs, the most commonly used variation operator is the combinatorial Laplacian given by:

where ${\bf D}$ is the diagonal degree matrix $\text{diag}\{d_{1},\dots,d_{N}\}$ with $d_{i}=\sum_{j}w_{ij}$. Since, $w_{ij}=w_{ji}$ for undirected graphs, this matrix is symmetric. As a result, it has real non-negative eigenvalues $\lambda_{i}\geq 0$ and an orthogonal set of eigenvectors. The variation functional associated with this operator is known as the graph Laplacian quadratic form and is given by

One can normalize the combinatorial Laplacian to obtain the symmetric normalized Laplacian and the (asymmetric) random walk Laplacian given as

Both ${\bf L}_{\textit{sym}}$ and ${\bf L}_{\textit{rw}}$ have non-negative eigenvalues. However the eigenvectors of ${\bf L}_{\textit{rw}}$ are not orthogonal as it is asymmetric. The eigenvectors of ${\bf L}_{\textit{sym}}$, on the other hand, are orthogonal. The variation functional associated with ${\bf L}_{\textit{sym}}$ has a nice interpretation as it normalizes the signal values on the nodes by the degree:

II-C2 Variation on directed graphs

Note that variation operators defined for directed graphs can also be used for undirected graphs since each undirected edge can be thought of as two oppositely pointing directed edges.

where $p=1,2$ and $\mu_{\text{max}}$ denotes the eigenvalue of ${\bf W}$ with the largest magnitude. It can be shown that for two eigenvalues $|\mu_{i}|<|\mu_{j}|$ of ${\bf W}$, the corresponding eigenvectors ${\bf v}_{i}$ and ${\bf v}_{j}$ satisfy $\text{Var}^{p}_{TV}({\bf v}^{i})<\text{Var}^{p}_{TV}({\bf v}^{j})$. In order to be consistent with our convention, one can define the variation operator as ${\bf L}=\mathbf{I}-{\bf W}/|\mu_{\text{max}}|$ which has the same eigenvectors as ${\bf W}$ with eigenvalues $\lambda_{i}=1-\mu_{i}/|\mu_{\text{max}}|$. This allows us to have the same ordering for the graph frequencies and the variations in the basis vectors. Note that for directed graphs, where ${\bf W}$ is not symmetric, the GFT basis vectors will not be orthogonal. Further, for some adjacency matrices, there may not exist a complete set of linearly independent eigenvectors. In such cases, one can use generalized eigenvectors in the Jordan normal form of ${\bf W}$ as stated before .

This notion of variation is based on the hub-authority model for specific directed graphs such as a hyperlinked environment (e.g., the web). This model distinguishes between two types of nodes. Hub nodes ${\cal H}$ are the subset of nodes which point to other nodes, whereas authority nodes ${\cal A}$ are the nodes to which other nodes point. Note that a node can be both a hub and an authority simultaneously. In a directed network, we need to define two kinds of degrees for each node $i\in{\cal V}$, namely the in-degree $p_{i}=\sum_{j}w_{ji}$ and the out-degree $q_{i}=\sum_{j}w_{ij}$. The co-linkage between two authorities $i,j\in{\cal A}$ or two hubs $i,j\in{\cal H}$ is defined as

respectively, and can be thought of as a cumulative link weight between two authorities (or hubs). Based on this, one can define a variation functional for a signal ${\bf f}$ on the authority nodes as

In order to write the above functional in a matrix form, define ${\bf T}={\bf D}_{q}^{-1/2}{\bf W}{\bf D}_{p}^{-1/2}$, where ${\bf D}_{p}^{-1/2}$ and ${\bf D}_{q}^{-1/2}$ are diagonal matrices with

It is possible to show that $\text{Var}_{\cal A}({\bf f})={\bf f}^{\top}{\bf L}_{\cal A}{\bf f}$, where ${\bf L}_{\cal A}=\mathbf{I}-{\bf T}^{\top}{\bf T}$. A variation functional for a signal ${\bf f}$ on the hub nodes can be defined in the same way as (10) and can be written in a matrix form as $\text{Var}_{\cal H}({\bf f})={\bf f}^{\top}{\bf L}_{\cal H}{\bf f}$, where ${\bf L}_{\cal H}=\mathbf{I}-{\bf T}{\bf T}^{\top}$. A convex combination $\text{Var}_{\gamma}({\bf f})=\gamma\text{Var}_{\cal A}({\bf f})+(1-\gamma)\text{Var}_{\cal H}({\bf f})$, with $\gamma\in$, can be used to define a variation functional for ${\bf f}$ on the whole vertex set ${\cal V}$. Note that the corresponding variation operator ${\bf L}_{\gamma}=\gamma{\bf L}_{\cal A}+(1-\gamma){\bf L}_{\cal H}$ is symmetric and positive semi-definite. Hence, eigenvectors and eigenvalues of ${\bf L}_{\gamma}$ can be used to define an orthogonal GFT similar to the undirected case, where the variation in the eigenvector increases as the corresponding eigenvalue increases.

Every directed graph has an associated random walk with a probability transition matrix ${\bf P}$ given by

By the Perron-Frobenius theorem, if ${\bf P}$ is irreducible then it has a stationary distribution $\pi$ which satisfies \hbox{\boldmath\pi}{\bf P}=\hbox{\boldmath\pi} . One can then define the following variation functional for signals on directed graphs :

Note that if the graph is undirected, the above expression reduces to (7) since, in that case, \hbox{\boldmath\pi}_{i}=d_{i}/\sum_{j}d_{j}. Intuitively, \hbox{\boldmath\pi}_{i}{\bf P}_{ij} can be thought of as the probability of transition from node $i$ to $j$ in the steady state. We expect it to be large if $i$ is similar to $j$. Thus, a big difference in signal values on nodes similar to each other contributes more to the variation. A justification for the above functional in terms of generalization of normalized cut to directed graphs is given in . Let \hbox{\boldmath\Pi}={\hbox{diag}}\{\hbox{\boldmath\pi}_{1},\ldots,\hbox{\boldmath\pi}_{n}\}. Then $\text{Var}_{\text{rw}}({\bf f})$ can be written as ${\bf f}^{\top}{\bf L}{\bf f}$, where

It is easy to see that the above ${\bf L}$ is a symmetric positive semi-definite matrix. Therefore, its eigenvectors can be used to define an orthonormal GFT, where the variation in the eigenvector increases as the corresponding eigenvalue increases.

Table I summarizes different choices of GFT bases based on the above variation operators. Our theory applies to all of these choices of GFT (with the caveat that diagonalizability is assumed in the definition of adjacency-based GFT).

III Sampling theory for graph signals

In this section, we address the issue of uniqueness and stability of bandlimited graph signal reconstruction and discuss different optimality criteria for sampling set selection assuming that the graph Fourier basis (i.e., the spectrum of the corresponding variation operator) is known. The uniqueness conditions in this section are equivalent to the ones in . However, the specific form in which we present these conditions lets us give a GFT-free definition of cutoff frequency. This together with the spectral proxies defined later in Section IV allows us to circumvent the explicit computation of the graph Fourier basis to ensure uniqueness and find a good sampling set.

The results in this section are useful when the graphs under consideration are small and thus, computing the spectrum of their variation operators is computationally feasible. They also serve as a guideline for tackling the aforementioned questions when the graphs are large and computation and storage of the graph Fourier basis is impractical.

In order to give a necessary and sufficient condition for unique identifiability of any signal ${\bf f}\in PW_{\omega}(G)$ from its samples ${\bf f}_{\cal S}$ on the sampling set ${\cal S}$, we first state the concept of uniqueness set .

A subset of nodes ${\cal S}$ is a uniqueness set for the space $PW_{\omega}(G)$ iff ${\bf x}_{\cal S}={\bf y}_{\cal S}$ implies ${\bf x}={\bf y}$ for all ${\bf x},{\bf y}\in PW_{\omega}(G)$.

Unique identifiability requires that no two bandlimited signals have the same samples on the sampling set as ensured by the following theorem in our previous work .

${\cal S}$ is a uniqueness set for $PW_{\omega}(G)$ if and only if $PW_{\omega}(G)\cap L_{2}({\cal S}^{c})=\{{\bf 0}\}$.

Let ${\cal R}=\{1,\dots,r\}$, where $\lambda_{r}$ is the largest graph frequency less than $\omega$. Then ${\cal S}$ is a uniqueness set for $PW_{\omega}(G)$ if and only if ${\bf U}_{{\cal S}{\cal R}}$ has full column rank.

If ${\bf U}_{{\cal S}{\cal R}}$ has a full column rank, then a unique reconstruction $\hat{{\bf f}}\in PW_{\omega}(G)$ can be obtained by finding the unique least squares solution to ${\bf f}_{\cal S}={\bf U}_{{\cal S}{\cal R}}{\bf c}$:

III-B Issue of Stability and Choice of Sampling set

Note that selecting a sampling set ${\cal S}$ for $PW_{\omega}(G)$ amounts to selecting a set of rows of ${\bf U}_{{\cal V}{\cal R}}$. It is always possible to find a sampling set of size $r=\dim{PW_{\omega}(G)}$ that uniquely determines signals in $PW_{\omega}(G)$ as proven below.

For any $PW_{\omega}(G)$, there always exists a uniqueness set ${\cal S}$ of size $|{\cal S}|=r$.

Since $\{{\bf u}^{i}\}_{i=1}^{r}$ are linearly independent, the matrix ${\bf U}_{{\cal V}{\cal R}}$ has full column rank equal to $r$. Further, since the row rank of a matrix equals its column rank, we can always find a linearly independent set ${\cal S}$ of $r$ rows such that ${\bf U}_{{\cal S}{\cal R}}$ has full rank that equals $r$, thus proving our claim. ∎

In most cases picking $r$ nodes randomly gives a full rank ${\bf U}_{{\cal S}{\cal R}}$. However, all sampling sets of given size are not equally good. A bad choice of ${\cal S}$ can give an ill-conditioned ${\bf U}_{{\cal S}{\cal R}}$ which in turn leads to an unstable reconstruction $\hat{{\bf f}}$. Stability of reconstruction is important when the true signal ${\bf f}$ is only approximately bandlimited (which is the case for most signals in practice) or when the samples are noisy. The reconstruction error in this case depends not only on noise and model mismatch but also on the choice of sampling set. The best sampling set achieves the smallest reconstruction error.

Since ${\bf f}\in PW_{\omega}(G)$, ${\bf U}_{{\cal V}{\cal R}}{\bf U}_{{\cal S}{\cal R}}^{+}{\bf f}_{\cal S}={\bf f}$. The reconstruction error equals ${\bf e}=\hat{{\bf f}}-{\bf f}={\bf U}_{{\cal V}{\cal R}}{\bf U}_{{\cal S}{\cal R}}^{+}{\bf n}$. If we assume that the entries of ${\bf n}$ are iid with zero mean and unit variance, then the covariance matrix of the reconstruction error is given by

Different costs can be defined to measure the reconstruction error as a function of the error covariance matrix. These cost functions are based on optimal design of experiments . If we define the optimal sampling set ${\cal S}^{\text{opt}}$ of size $m$, as the set which minimizes the mean squared error, then assuming ${\bf U}_{{\cal V}{\cal R}}$ has orthonormal columns, we have

This is analogous to the so-called $A$-optimal design. Similarly, minimizing the maximum eigenvalue of the error covariance matrix leads to $E$-optimal design. For an orthonormal ${\bf U}_{{\cal V}{\cal R}}$, the optimal sampling set with this criterion is given by

where $\sigma_{\min}(.)$ denotes the smallest singular value of a matrix. It can be thought of as a sampling set which minimizes the worst case reconstruction error. The above criterion is equivalent to the one proposed in . Further, one can show that when ${\bf U}_{{\cal V}{\cal R}}$ does not have orthonormal columns, (17) and (18) produce sampling sets that minimize upper bounds on the mean squared and worst case reconstruction errors respectively. Note that both $A$ and $E$-optimality criteria lead to combinatorial problems, but it is possible to develop greedy approximate solutions to these problems.

So far we assumed that the true signal ${\bf f}\in PW_{\omega}(G)$ and hence, ${\bf U}_{{\cal V}{\cal R}}{\bf U}_{{\cal S}{\cal R}}^{+}{\bf f}_{\cal S}={\bf f}$. However, in most applications, the signals are only approximately bandlimited. The reconstruction error in such a case is analyzed next.

III-B2 Effect of model mismatch

Let ${\bf P}={\bf U}_{{\cal V}{\cal R}}{\bf U}_{{\cal V}{\cal R}}^{\top}$ be the projector for $PW_{\omega}(G)$ and ${\bf Q}={\bf S}{\bf S}^{\top}$ be the projector for $L_{2}({\cal S})$. Assume that the true signal is given by ${\bf f}={\bf f}^{*}+\Delta{\bf f}$, where ${\bf f}^{*}={\bf P}{\bf f}$ is the bandlimited component of the signal and $\Delta{\bf f}={\bf P}^{\perp}{\bf f}$ captures the “high-pass component” (i.e., the model mismatch). If we use (14) for reconstructing ${\bf f}$, then a tight upper bound on the reconstruction error is given by

where $\theta_{\max}$ is the maximum angle between subspaces $PW_{\omega}(G)$ and $L_{2}({\cal S})$ defined as

$\cos(\theta_{\max})>0$ when the uniqueness condition in Theorem 1 is satisfied and the error is bounded. Intuitively, the above equation says that for the worst case error to be minimum, the sampling and reconstruction subspaces should be as aligned as possible.

We define an optimal sampling set ${\cal S}^{\text{opt}}$ of size $m$ for $PW_{\omega}(G)$ as the set which minimizes the worst case reconstruction error. Therefore, $L_{2}({\cal S}^{\text{opt}})$ makes the smallest maximum angle with $PW_{\omega}(G)$. It is easy to show that $\cos(\theta_{\max})=\sigma_{\min}({\bf U}_{{\cal S}{\cal R}})$. Thus, to find this set we need to solve a similar problem as (18). As stated before, this problem is combinatorial. It is possible to give a greedy algorithm to get an approximate solution. A simple greedy heuristic to approximate ${\cal S}^{\text{opt}}$ is to perform column-wise Gaussian elimination over ${\bf U}_{{\cal V}{\cal R}}$ with partial row pivoting. The indices of the pivot rows in that case form a good estimate of ${\cal S}^{\text{opt}}$ in practice.

Table II summarizes the different set selection criteria and corresponding search algorithms under various assumptions about the signal. However, the methods described above require computation of many eigenvectors of the variation operator ${\bf L}$. We circumvent this issue in the next section, by defining graph spectral proxies based on powers of ${\bf L}$. These spectral proxies do not require eigen-decomposition of ${\bf L}$ and still allow us to define a measure of quality of sampling sets. As we will show, these proxies arise naturally in the expression for the bound on the reconstruction error. Thus, a sampling set optimal with respect to these spectral proxies ensures a small reconstruction error bound.

IV Sampling Set Selection using Graph Spectral Proxies

As discussed earlier, graphs considered in most real applications are very large. Hence, computing and storing the graph Fourier basis explicitly is often practically infeasible. We now present techniques that allow us to express the condition for unique bandlimited reconstruction and methods for sampling set selection via simple operations using the variation operator. The following discussion holds for any choice of the variation operator ${\bf L}$ in Table I.

In order to obtain a measure of quality for a sampling set ${\cal S}$, we first find the cutoff frequency associated with it, which can be defined as the largest frequency $\omega$ such that ${\cal S}$ is a uniqueness set for $PW_{\omega}(G)$. It follows from Theorem 1 that, for ${\cal S}$ to be a uniqueness set of $PW_{\omega}(G)$, $\omega$ needs to be less than the minimum possible bandwidth that a signal in $L_{2}({\cal S}^{c})$ can have. This would ensure that no signal from $L_{2}({\cal S}^{c})$ can be a part of $PW_{\omega}(G)$. Thus, the cutoff frequency $\omega_{c}({\cal S})$ for a sampling set ${\cal S}$ can be expressed as:

To use the equation above, we first need a tool to approximately compute the bandwidth $\omega(\phi)$ of any given signal $\phi$ without computing the Fourier coefficients explicitly. Our proposed method for bandwidth estimation is based on the following definition:

For an operator ${\bf L}$ with real eigenvalues and eigenvectors, $\omega_{k}({\bf f})$ can be shown to increase monotonically with $k$:

These quantities are bounded from above, as a result, $\lim_{k\rightarrow\infty}\omega_{k}({\bf f})$ exists for all ${\bf f}$. Consequently, it is easy to prove that if $\omega({\bf f})$ denotes the bandwidth of a signal ${\bf f}$, then

Note that (24) also holds for an asymmetric ${\bf L}$ that has complex eigenvalues and eigenvectors. The proofs of (23) and (24) are provided in the Appendix. These properties give us an important insight: as we increase the value of $k$, the spectral proxies tend to have a value close to the actual bandwidth of the signal, i.e., they essentially indicate the frequency localization of the signal energy. Therefore, using $\omega_{k}(\phi)$ as a proxy for $\omega(\phi)$ (i.e. bandwidth of $\phi$) is justified and this leads us to define the cut-off frequency estimate of order k as

Using the definitions of $\Omega_{k}({\cal S})$ and $\omega_{c}({\cal S})$ along with (23) and (24), we conclude that for any $k_{1}<k_{2}$:

Using (26) and (21), we now state the following proposition:

For any $k$, ${\cal S}$ is a uniqueness set for $PW_{\omega}(G)$ if, $\omega<\Omega_{k}({\cal S})$. $\Omega_{k}({\cal S})$ can be computed from (25) as

where $\sigma_{1,k}$ denotes the smallest eigenvalue of the reduced matrix $(({\bf L}^{\top})^{k}{\bf L}^{k})_{{\cal S}^{c}}$. Further, if $\psi_{1,k}$ is the corresponding eigenvector, and $\phi_{k}^{*}$ minimizes $\omega_{k}(\phi)$ in (25) (i.e. it approximates the smoothest possible signal in $L_{2}({\cal S}^{c})$), then

We note from (26) that to get a better estimate of the true cut-off frequency, one simply needs a higher $k$. Therefore, there is a trade-off between accuracy of the estimate on the one hand, and complexity and numerical stability on the other (that arise by taking higher powers of ${\bf L}$).

IV-B Best Sampling Set of Given Size

As shown in Proposition 2, $\Omega_{k}({\cal S})$ is an estimate of the smallest bandwidth that a signal in $L_{2}({\cal S}^{c})$ can have and any signal in $PW_{\omega}(G)$ is uniquely sampled on ${\cal S}$ if $\omega<\Omega_{k}({\cal S})$. Intuitively, we would like the projection of $L_{2}({\cal S}^{c})$ along $PW_{\omega}(G)$ to be as small as possible. Based on this intuition, we propose the following optimality criterion for selecting the best sampling set of size $m$:

To motivate the above criterion more formally, let ${\bf P}$ denote the projector for $PW_{\omega}(G)$. The minimum gap between the two subspaces $L_{2}({\cal S}^{c})$ and $PW_{\omega}(G)$ is given by:

We now show that $\Omega_{k}({\cal S})$ also arises in the bound on the reconstruction error when the reconstruction is obtained by variational energy minimization:

It was shown in that if ${\bf f}\in PW_{\omega}(G)$, then the reconstruction error $\|\hat{{\bf f}}_{m}-{\bf f}\|/\|{\bf f}\|$, for a given $m$, is upper-bounded by $2(\omega/\Omega_{1}({\cal S}))^{m}$. This bound is suboptimal and can be improved by replacing $\Omega_{1}({\cal S})$ with $\Omega_{k}({\cal S})$ (which, from (26), is at least as large as $\Omega_{1}({\cal S})$) for any $k\leq m$, as shown in the following theorem:

Let $\hat{{\bf f}}_{m}$ be the solution to (31) for a signal ${\bf f}\in PW_{\omega}(G)$. Then, for any $k\leq m$,

Note that $(\hat{{\bf f}}_{m}-{\bf f})\in L_{2}({{\cal S}^{c}})$. Therefore, from (25)

(33) follows from triangle inequality. (34) holds because $\hat{{\bf f}}_{m}$ minimizes $\|{\bf L}^{m}\hat{{\bf f}}_{m}\|$ over all sample consistent signals. (35) follows from the definition of $\omega_{m}({\bf f})$ and the last step follows from (24) and (26). ∎

Note that for the error bound in (32) to go to zero as $m\to\infty$, $\omega$ must be less than $\Omega_{k}({\cal S})$. Thus, increasing $\Omega_{k}({\cal S})$ allows us to reconstruct signals in a larger bandlimited space using the variational method. Moreover, for a fixed $m$ and $k$, a higher value of $\Omega_{k}({\cal S})$ leads to a lower reconstruction error bound. The optimal sampling set ${\cal S}^{\rm opt}_{k}$ in (29) essentially minimizes this error bound.

IV-C Finding the Best Sampling Set

The problem posed in (29) is a combinatorial problem because we need to compute $\Omega_{k}({\cal S})$ for every possible subset ${\cal S}$ of size $m$. We therefore formulate a greedy heuristic to get an estimate of the optimal sampling set. Starting with an empty sampling set ${\cal S}$ ($\Omega_{k}({\cal S})=0$) we keep adding nodes (from ${\cal S}^{c}$) one-by-one while trying to ensure maximum increase in $\Omega_{k}({\cal S})$ at each step. To achieve this, we first consider the following quantity:

where ${\bf 1}_{\cal S}:{\cal V}\rightarrow\{0,1\}$ denotes the indicator function for the subset ${\cal S}$ (i.e. ${\bf 1}({\cal S})={\bf 1}$ and ${\bf 1}({\cal S}^{c})={\bf 0}$). Note that the right hand side of (36) is simply a relaxation of the constraint in (25). When $\alpha\gg 1$, the components ${\bf x}({\cal S})$ are highly penalized during minimization, hence, forcing values of ${\bf x}$ on ${\cal S}$ to be vanishingly small. Thus, if ${\bf x}^{\alpha}_{k}({\bf 1}_{\cal S})$ is the minimizer in (36), then $[{\bf x}^{\alpha}_{k}({\bf 1}_{\cal S})]({\cal S})\rightarrow{\bf 0}$. Therefore, for $\alpha\gg 1$,

Now, to tackle the combinatorial nature of our problem, we allow a binary relaxation of the indicator ${\bf 1}_{S}$ in (36), to define the following quantities

When ${\bf t}={\bf 1}_{\cal S}$, we know that the minimizer of (39) with respect to ${\bf x}$ for large $\alpha$ is $\phi_{k}^{*}$. Hence,

The equation above gives us the desired greedy heuristic - starting with an empty ${\cal S}$ (i.e., ${\bf 1}_{\cal S}={\bf 0}$), if at each step, we include the node on which the smoothest signal $\phi_{k}^{*}\in L_{2}({\cal S}^{c})$ has maximum energy (i.e. ${\bf 1}_{\cal S}(i)\leftarrow 1,i=\text{arg max}_{j}(\phi_{k}^{*}(j))^{2}$), then $\lambda^{\alpha}_{k}({\bf t})$ and in effect, the cut-off estimate $\Omega_{k}({\cal S})$, tend to increase maximally. We summarize the method for estimating ${\cal S}^{\text{opt}}_{k}$ in Algorithm 1.

One can show that the cutoff frequency estimate $\Omega_{k}({\cal S})$ associated with a sampling set can only increase (or remain unchanged) when a node is added to it. This is stated more formally in the following proposition.

Let ${\cal S}_{1}$ and ${\cal S}_{2}$ be two subsets of nodes of $G$ with ${\cal S}_{1}\subseteq{\cal S}_{2}$. Then $\Omega_{k}({\cal S}_{1})\leq\Omega_{k}({\cal S}_{2})$.

This turns out to be a straightforward consequence of the eigenvalue interlacing property for symmetric matrices.

Let ${\bf B}$ be a symmetric $n\times n$ matrix. Let ${\cal R}=\{1,2,\ldots,r\}$, for $1\leq r\leq n-1$ and ${\bf B}_{r}={\bf B}_{{\cal R}}$. Let $\lambda_{k}({\bf B}_{r})$ be the $k$-th largest eigenvalue of ${\bf B}_{r}$. Then the following interlacing property holds:

The above theorem implies that if ${\cal S}_{1}\subseteq{\cal S}_{2}$, then ${\cal S}_{2}^{c}\subseteq{\cal S}_{1}^{c}$ and thus, $\lambda_{\min}\left[\left(({\bf L}^{\top})^{k}{\bf L}^{k}\right)_{{\cal S}^{c}_{1}}\right]\leq\lambda_{\min}\left[\left(({\bf L}^{\top})^{k}{\bf L}^{k}\right)_{{\cal S}^{c}_{2}}\right]$.

From Section III, we know that the optimal sampling set can be obtained by maximizing $\sigma_{\text{min}}\left({\bf U}_{{\cal S}{\cal R}}\right)$ with respect to ${\cal S}$. A heuristic to obtain good sampling sets is to perform a column-wise Gaussian elimination with pivoting on the eigenvector matrix ${\bf U}$. Then, a sampling set of size $i$ is given by the indices of zeros in the $(i+1)^{\text{th}}$ column of the echelon form. We now show that the greedy heuristic proposed in Algorithm 1 is closely related to this rank-revealing Gaussian elimination procedure through the following observation:

Let $\Phi$ be the matrix whose columns are given by the smoothest signals $\phi_{\infty}^{*}$ obtained sequentially after each iteration of Algorithm 1 with $k=\infty$, (i.e., $\Phi=\left[\phi^{*}_{\infty}|_{|{\cal S}|=0}\;\phi^{*}_{\infty}|_{|{\cal S}|=1},\;\dots\right]$). Further, let ${\bf T}$ be the matrix obtained by performing column-wise Gaussian elimination on ${\bf U}$ with partial pivoting. Then, the columns of ${\bf T}$ are equal to the columns of $\Phi_{\infty}^{*}$ within a scaling factor.

If ${\cal S}$ is the smallest sampling set for uniquely representing signals in $PW_{\omega}(G)$ and $r={\rm dim}\;PW_{\omega}(G)$, then we have the following:

The smoothest signal $\phi^{*}_{\infty}\in L_{2}({\cal S}^{c})$ has bandwidth $\lambda_{r+1}$.

Therefore, $\phi^{*}_{\infty}|_{|{\cal S}|=r}$ is spanned by the first $r+1$ frequency basis elements $\{{\bf u}_{1},\dots,{\bf u}_{r+1}\}$. Further, since $\phi^{*}_{\infty}|_{|{\cal S}|=r}$ has zeroes on exactly $r$ locations, it can be obtained by performing Gaussian elimination on ${\bf u}_{r+1}$ using ${\bf u}_{1},{\bf u}_{2},\dots,{\bf u}_{r}$. Hence the $(r+1)^{\text{th}}$ column of $\Phi$ is equal (within a scaling factor) to the $(r+1)^{\text{th}}$ column of ${\bf T}$. Pivoting comes from the fact that the $(i+1)^{\text{th}}$ sampled node is given by the index of the element with maximum magnitude in $\phi_{\infty}^{*}|_{|{\cal S}|=i}$, and is used as the pivot to zeros out elements with same index in subsequent columns. ∎

The above result illustrates that Algorithm 1 is an iterative procedure that approximates a rank-revealing Gaussian elimination procedure on ${\bf U}_{{\cal V}{\cal R}}$. For the known-spectrum case, this is a good heuristic for maximizing $\sigma_{min}\left({\bf U}_{{\cal S}{\cal R}}\right)$. In other words, our method directly maximizes $\sigma_{min}\left({\bf U}_{{\cal S}{\cal R}}\right)$ without going through the intermediate step of computing ${\bf U}_{{\cal V}{\cal R}}$. As we shall see in the next subsection, this results in significant savings in both time and space complexity.

IV-D Complexity and implementation issues

We note that in the algorithm, computing the first eigen-pair of $(({\bf L}^{\top})^{k}{\bf L}^{k})_{{\cal S}^{c}}$ is the major step for each iteration. There are many efficient iterative methods, such as those based on Rayleigh quotient minimization, for computing the smallest eigen-pair of a matrix . The atomic step in all of these methods consists of matrix-vector products. Specifically, in our case, this step involves evaluating the expression $(({\bf L}^{\top})^{k}{\bf L}^{k})_{{\cal S}^{c}}{\bf x}$. Note that we do not actually need to compute the matrix $(({\bf L}^{\top})^{k}{\bf L}^{k})_{{\cal S}^{c}}$ explicitly, since the expression can be implemented as a sequence of matrix-vector products as

Evaluating the expression involves $2k$ matrix-vector products and has a complexity of $O(k|E|)$, where $|E|$ is the number of edges in the graph. Moreover, a localized and parallel implementation of this step is possible in the case of sparse graphs. The number of iterations required for convergence of the eigen-pair computation iterations is a function of the eigen-value gaps and hence dependent on the graph structure and edge-weights.

For the methods of and , one needs to compute a portion of the eigenvector matrix, i.e., ${\bf U}_{{\cal V}{\cal S}}$ (assuming $|{\cal R}|=|{\cal S}|$). This can be done using block-based Rayleigh quotient minimization methods , block-based Kryolov subspace methods such as Arnoldi/Lanczos iterations or deflation methods in conjunction with single eigen-pair solvers . The complexity of these methods increases considerably as the number of requested eigen-pairs increases, making them impractical. On the other hand, our method requires computing a single eigen-pair at each iteration, making it viable for cases when a large number of samples are required. Moreover, the sample search steps in the methods of and require an SVD solver and a linear system solver, respectively, thus making them much more complex in comparison to our method, where we only require finding the maximum element of a vector. Our algorithm is also efficient in terms of space complexity, since at any point we just need to store ${\bf L}$ and one vector. On the other hand, require storage of at least $|{\cal S}|$ eigenvectors.

A summary of the complexities of all the methods is given in Table III. The eigen-pair computations for are assumed to be performed using a block version of the Rayleigh quotient minimization method, which has a complexity of $O((|E||{\cal S}|+C|{\cal S}|^{3})T_{1})$, where $T_{1}$ denotes the number of iterations for convergence, and $C$ is a constant. The complexity of computing one eigen-pair in our method is $O(k|E||S|T_{2}(k))$, where $T_{2}(k)$ denotes the average number of iterations required for convergence of a single eigen-pair. $T_{1}$ and $T_{2}(k)$ required to achieve a desired error tolerance are functions of the eigen-gaps of ${\bf L}$ and ${\bf L}^{k}$ respectively. In general, $T_{2}(k)>T_{1}$, since ${\bf L}^{k}$ has lower eigengaps near the smallest eigenvalue. Increasing the parameter $k$ further flattens the spectrum of ${\bf L}^{k}$ near the smallest eigenvalue leading to an increase in $T_{2}(k)$, since one has to solve a more ill-conditioned problem. We illustrate this in the next section through experiments that compare the running times of all the methods.

The choice of the parameter $k$ depends on the desired accuracy – a larger value of $k$ gives a better sampling set, but increases the complexity proportionally, thus providing a trade-off. Through experiments, we show in the next section that the quality of the sampling set is more sensitive to choice of $k$ for sparser graphs. This is because increasing $k$ results in the consideration of more global information while selecting samples. On the other hand, dense graphs have a lower diameter and there is relatively little information to be gained by increasing $k$.

V Experiments

We now numerically evaluate the performance of the proposed work. The experiments involve comparing the reconstruction errors and running times of different sampling set selection algorithms in conjunction with consistent bandlimited reconstruction (14)Although reconstruction using (14) requires explicit computation of ${\bf U}_{{\cal V}{\cal R}}$, there exist efficient localized reconstruction algorithms that circumvent this . However, in the current work, we restrict our attention to the problem of sampling set selection.. We compare our approach with the following methods:

This method uses a greedy algorithm to approximate the ${\cal S}$ that maximizes $\sigma_{\min}({\bf U}_{{\cal S}{\cal R}})$. Consistent bandlimited reconstruction (14) is then used to estimate the unknown samples.

At each iteration $i$, this method finds the representation of ${\bf u}_{i}$ as $\sum_{j<i}\beta_{j}{\bf u}_{j}+\sum_{u\notin{\cal S}}\alpha_{u}{\bf 1}_{u}$, where ${\bf 1}_{u}$ is the delta function on $u$. The node $v$ with maximum $|\alpha_{v}|$ is sampled. Reconstruction is done using (14).

Both the above methods assume that a portion of the frequency basis is known and the signal to be recovered is exactly bandlimited. As a baseline, we also compare all sampling set selection methods against uniform random sampling.

We first give some simple examples on the following simulated undirected graphs:

Erdös-Renyi random graph (unweighted) with $1000$ nodes and connection probability $0.01$.

Small world graph (unweighted) with $1000$ nodes. The underlying regular graph with degree $8$ is rewired with probability $0.1$.

Barabási-Albert random network with $1000$ nodes. The seed network is a fully connected graph with $m_{0}=4$ vertices, and each new vertex is connected to $m=4$ existing vertices randomly. This model, as opposed to $G1$ and $G2$, is a scale-free network, i.e., its degrees follow a power law $P(k)\sim k^{-3}$.

The performance of the sampling methods depends on the assumptions about the true signal and sampling noise. For each of the above graphs, we consider the problem in the following scenarios:

The true signal is noise-free and exactly bandlimited with $r=\dim{PW_{\omega}(G)}=50$. The non-zero GFT coefficients are randomly generated from ${\cal N}(1,0.5^{2})$.

The true signal is exactly bandlimited with $r=50$ and non-zero GFT coefficients are generated from ${\cal N}(1,0.5^{2})$. The samples are noisy with additive iid Gaussian noise such that the SNR equals 20dB.

The true signal is approximately bandlimited with an exponentially decaying spectrum. Specifically, the GFT coefficients are generated from ${\cal N}(1,0.5^{2})$, followed by rescaling with the following filter (where $r=50$):

We generate 50 signals from each of the three signal models on each of the graphs, use the sampling sets obtained from the all the methods to perform reconstruction and plot the mean of the mean squared error (MSE) for different sizes of sampling sets. For our algorithm, we set the value of $k$ to 2, 8 and 14. The result is illustrated in Figure 1. Note that when the size of the sampling set is less than $r=50$, the results are quite unstable. This is expected, because the uniqueness condition is not satisfied by the sampling set. Beyond $|{\cal S}|=r$, we make the following observations:

For the noise-free, bandlimited signal model F1, all methods lead to zero reconstruction error as soon as the size of the sampling set exceeds the signal cutoff $r=50$ (error plots for this signal model are not shown). This is expected from the sampling theorem. It is interesting to note that in most cases, uniform random sampling does equally well, since the signal is noise-free and perfectly bandlimited.

For the noisy signal model F2 and the approximately bandlimited model F3, our method has better or comparable performance in most cases. This indicates that our method is fairly robust to noise and model mismatch. Uniform random sampling performs very badly as expected, because of lack of stability considerations.

Parameter $k$ in the definition of spectral proxies controls how closely we estimate the bandwidth of any signal ${\bf f}$. Spectral proxies with higher values of $k$ give a better approximation of the bandwidth. Our sampling set selection algorithm tries to maximize the smallest bandwidth that a signal in $L_{2}({{\cal S}^{c}})$ can have. Using higher values of $k$ allows us to estimate this smallest bandwidth more closely, thereby leading to better sampling sets as demonstrated in Figure 2. Intuitively, maximizing $\Omega_{k}({\cal S})$ with $k=1$ ensures that the sampled nodes are well connected to the unsampled nodes and thus, allows better propagation of the observed signal information. Using $k>1$ takes into account multi-hop paths while ensuring better connectedness between ${\cal S}$ and ${{\cal S}^{c}}$. This effect is especially important in sparsely connected graphs and the benefit of increasing $k$ becomes less noticeable when the graphs are dense as seen in Figure 2. However, this improvement in performance in the case of sparse graphs comes at the cost of increased numerical complexity.

Running time

We also compare the running times of the sampling set selection methods for different sizes of the graph. For our experiments, we generate symmetrized Erdös-Renyi random graphs of different sizes with parameter $0.01$, and measure the average running time of selecting $5\%$ of the samples in MATLAB. For computing the eigen-pairs, we use the code for the Locally Optimal Block Prec-conditioned Conjugate Gradient (LOBPCG) method available online (this was observed to be faster than MATLAB’s inbuilt sparse eigensolver eigs, which is based on Lanczos iterations). The results of the experiments are shown in Table IV. We observe that the rate of increase of running time as the graph size increases is slower for our method compared to other methods, thus making it more practical. Note that the increase with respect to $k$ is nonlinear since the eigengaps are a function of $k$ and lead to different number of iterations required for convergence of the eigenvectors.

V-B A Real World Example

We first compare the performance of the proposed method against M1 and M2 using the normalized adjacency matrix based GFT with the variation operator ${\bf L}=\mathbf{I}-{\bf D}^{-1}{\bf W}$. The bandwidth parameter $r$ is set to $50$. The plot of classification error averaged over the 10 dataset instances vs. number of labels is presented in Figure 3(a). It shows that the proposed method has comparable performance despite being localized. The performance is also affected by the choice of the variation operators (or, the GFT bases). Figure 3(b) shows that the variation operators based on the hub-authority model and random walk offer higher classification accuracy and thus, are more suited for this particular application. Their superior performance can be explained by looking at the signal representation in the respective GFT domains. Figure 3(c) shows the fraction of signal energy captured in increasing number of GFT coefficients starting from low frequency. Since the hub-authority model based GFT and random walk based GFT offer more energy compaction than adjacency based GFT, the signal reconstruction quality using these bases is naturally better.

VI Conclusion

We studied the problem of selecting an optimal sampling set for reconstruction of bandlimited graph signals. The starting point of our framework is the notion of the Graph Fourier Transform (GFT) which is defined via an appropriate variation operator. Our goal is to find good sampling sets for reconstructing signals which are bandlimited in the above frequency domain. We showed that when the samples are noisy or the true signal is only approximately bandlimited, the reconstruction error depends not only on the model mismatch but also on the choice of sampling set. We proposed a measure of quality for the sampling sets, namely the cutoff frequency, that can be computed without finding the GFT basis explicitly. A sampling set that maximizes the cutoff frequency is shown to minimize the reconstruction error. We also proposed a greedy algorithm which finds an approximately optimal set. The proposed algorithm can be efficiently implemented in a distributed and parallel fashion. Together with localized signal reconstruction methods, it gives an effective method for sampling and reconstruction of smooth graph signals on large graphs.

The present work opens up some new questions for future research. The problem of finding a sampling set with maximum cutoff frequency is combinatorial. The proposed greedy algorithm gives only an approximate solution to this problem. It would be useful to find a polynomial time algorithm with theoretical guarantees on the quality of approximation. Further, the proposed set selection method is not adaptive, i.e., the choice of sampling locations does not depend on previously observed samples. This can be a limitation in applications that require batch sampling. In such cases, it would be desirable to have an adaptive sampling set selection scheme which takes into account the previously observed samples to refine the choice of nodes to be sampled in the future.

In this section, we prove the monotonicity and convergence properties of $\omega_{k}({\bf f})$.

If ${\bf L}$ has real eigenvalues and eigenvectors, then for any $k_{1}<k_{2}$, we have $\omega_{k_{1}}({\bf f})\leq\omega_{k_{2}}({\bf f}),\forall{\bf f}$.

We first expand $\omega_{k_{1}}({\bf f})$ as follows:

If ${\bf L}$ has real entries, but complex eigenvalues and eigenvectors, then these occur in conjugate pairs, hence, the above summation is real. However, in that case, $\omega_{k}({\bf f})$ is not guaranteed to increase in a monotonous fashion, since $c_{ij}$’s are not real and Jensen’s inequality breaks down. ∎

Let $\omega({\bf f})$ be the bandwidth of any signal ${\bf f}$. Then, the following holds:

We first consider the case when ${\bf L}$ has real eigenvalues and eigenvectors. Let $\omega({\bf f})=\lambda_{p}$, then we have:

Now, if ${\bf L}$ has complex eigenvalues and eigenvectors, then these have to occur in conjugate pairs since ${\bf L}$ has real entries. Hence, for this case, we do a similar expansion as above and take $|\lambda_{p}|$ out of the expression. Then, the limit of the remaining term is once again equal to 1. ∎

References