Graph sampling with determinantal processes

TL;DR

This paper introduces a novel graph sampling method using determinantal point processes, enabling efficient and accurate signal recovery on large graphs, especially those with community structures, by leveraging spectral properties or random walks.

Contribution

It proposes a new DPP-based sampling strategy for graph signals that works efficiently on large graphs without requiring spectral access, with theoretical and empirical validation.

Findings

Perfect recovery on small graphs with accessible spectrum.

Effective sampling on large graphs using loop-erased random walks.

Promising results on graphs with strong community structures.

Abstract

We present a new random sampling strategy for k-bandlimited signals defined on graphs, based on determinantal point processes (DPP). For small graphs, ie, in cases where the spectrum of the graph is accessible, we exhibit a DPP sampling scheme that enables perfect recovery of bandlimited signals. For large graphs, ie, in cases where the graph's spectrum is not accessible, we investigate, both theoretically and empirically, a sub-optimal but much faster DPP based on loop-erased random walks on the graph. Preliminary experiments show promising results especially in cases where the number of measurements should stay as small as possible and for graphs that have a strong community structure. Our sampling scheme is efficient and can be applied to graphs with up to $10^6$ nodes.

Figures

Equations

$$\mathbf{y}=\mathsf{M}\mathbf{x}+\mathbf{n}\in {\mathbb{R}}^m,$$

$${\mathbf{x}}_{\text{rec}}=\underset{\mathbf{z}\in \mathrm{span}({\mathsf{U}}_k)}{\mathrm{argmin}}{\Vert \mathsf{M}\mathbf{z}-\mathbf{y}\Vert }^2={\mathsf{U}}_k(\mathsf{M}{\mathsf{U}}_k{)}^{†}\mathbf{y},$$

$${\mathbf{x}}_{\text{rec}}$$

$$=\mathbf{x}+{\mathsf{U}}_k({\mathsf{U}}_k^{⊺}{\mathsf{M}}^{⊺}\mathsf{M}{\mathsf{U}}_k{)}^{-1}{\mathsf{U}}_k^{⊺}{\mathsf{M}}^{⊺}\mathbf{n}.\text{qed}$$

$${\mathcal{A}}^{\text{WCE}}=\arg \underset{\mathcal{A}\text{s.t.}|\mathcal{A}|=k}{\max }{\sigma }_1^2.$$

$${\mathcal{A}}^{\text{MSE}}=\arg \underset{\mathcal{A}\text{s.t.}|\mathcal{A}|=k}{\min }\text{tr}\left[({\mathsf{U}}_k^{⊺}{\mathsf{M}}^{⊺}\mathsf{M}{\mathsf{U}}_k{)}^{-1}\right]=\arg \underset{\mathcal{A}\text{s.t.}|\mathcal{A}|=k}{\min }\sum _{i=1}^k\frac{1}{{\sigma }_i^2},$$

$${\mathcal{A}}^{\text{MV}}=\arg \underset{\mathcal{A}\text{s.t.}|\mathcal{A}|=k}{\max }\prod _{i=1}^k{\sigma }_i^2,$$

$${\mathbf{x}}_{\text{rec}}=\underset{\mathbf{z}\in {\mathbb{R}}^N}{\mathrm{argmin}}{\Vert \mathsf{M}\mathbf{z}-\mathbf{y}\Vert }^2+\gamma {\mathbf{z}}^{⊺}{\mathsf{L}}^r\mathbf{z},$$

$$({\mathsf{M}}^{⊺}\mathsf{M}+\gamma {\mathsf{L}}^r){\mathbf{x}}_{\text{rec}}={\mathsf{M}}^{⊺}\mathbf{y},$$

$$\mathsf{P}=\text{diag}(m{p}_{{\omega }_1}^{\ast }|\dots |m{p}_{{\omega }_m}^{\ast }).$$

$${\mathbf{x}}_{\text{rec}}=\underset{\mathbf{z}\in \mathrm{span}({\mathsf{U}}_k)}{\mathrm{argmin}}{\Vert {\mathsf{P}}^{-1/2}\left(\mathsf{M}\mathbf{z}-\mathbf{y}\right)\Vert }^2;$$

$${\mathbf{x}}_{\text{rec}}=\underset{\mathbf{z}\in {\mathbb{R}}^N}{\mathrm{argmin}}{\Vert {\mathsf{P}}^{-1/2}\left(\mathsf{M}\mathbf{z}-\mathbf{y}\right)\Vert }^2+\gamma {\mathbf{z}}^{⊺}{\mathsf{L}}^r\mathbf{z}.$$

$$\mathbb{P}(\mathcal{S}\subseteq \mathcal{A})=\text{det}({\mathsf{K}}_{\mathcal{S}}),$$

$$\mathbb{E}(|\mathcal{A}|)=\text{tr}(\mathsf{K})=\sum _{i=1}^N{\mu }_i\text{and}\text{Var}(|\mathcal{A}|)=\sum _{i=1}^N{\mu }_i(1-{\mu }_i).$$

$$\mathsf{P}=\text{diag}({\pi }_{{\omega }_1}|\dots |{\pi }_{{\omega }_m})$$

$${\mathbb{E}}_{\mathcal{A}}$$

$$=\sum _{i=1}^N{x}_i^2\sum _{\mathcal{A}\supset i}\frac{\mathbb{P}(\mathcal{A})}{{\pi }_i}=\sum _{i=1}^N{x}_i^2={\Vert \mathbf{x}\Vert }^2\text{qed}$$

$${\mathsf{K}}_k={\mathsf{U}}_k{\mathsf{U}}_k^{⊺}=\mathsf{U}{h}_{{\lambda }_k}(\mathsf{\Lambda }){\mathsf{U}}^{⊺},$$

$$\text{det}({{\mathsf{K}}_k}_{\mathcal{A}})>0,\text{i.e.},\text{det}(\mathsf{M}{\mathsf{U}}_k{\mathsf{U}}_k^{⊺}{\mathsf{M}}^{⊺})=\prod _{i=1}^k{\sigma }_i^2>0,$$

Full text

I Introduction

Graphs are a central modelling tool for network-structured data. Data on a graph, called graph signals [1], such as individual hobbies in social networks, blood flow of brain regions in neuronal networks, or traffic at a station in transportation networks, may be represented by a scalar per node. Studying such signals with respect to the particular graph topology on which it is defined, is the goal of graph signal processing (GSP). One of the top challenges of GSP is not only to adapt classical processing tools to arbitrary topologies defined by the underlying graph, but to do so with efficient algorithms that can handle the large datasets encountered today. A possible answer to this efficiency challenge is dimension reduction, and in particular sampling. Graph sampling [2, 3, 4, 5, 6] consists in measuring a graph signal on a reduced set of nodes carefully chosen to optimize a predefined objective, may it be the estimation of some statistics, signal compression via graph filterbanks, signal recovery, etc.

In an effort to generalize Shannon’s sampling theorem of bandlimited signals, several authors [3, 4, 5] have proposed sampling schemes adapted to the recovery of a priori smooth graph signals, and in particular $k$-bandlimited signals, i.e., linear combinations of only $k$ graph Fourier modes. Classically, there are three main types of sampling schemes that one may try to adapt to graphs: periodic, irregular, and random.

Periodic graph sampling is ill-defined unless the graph is not bipartite (or multipartite) and has only been investigated in the design of graph filterbanks where the objective is to sample “one every two nodes” [7, 8].

Irregular graph sampling of $k$-bandlimited graph signals has been studied by Anis et al. [3] and Chen et al. [4] who, building upon Pensenson’s work [2], look for optimal sampling sets of $k$-bandlimited graphs, in the sense that they are tight (i.e., they contain only $k$ nodes) and that they maximize noise-robustness of the reconstruction. To exhibit such an optimal set, one needs i) to compute the first $k$ eigenvectors of the Laplacian which becomes prohibitive at large scale, and ii) to optimize over all subsets of $k$ nodes, which is a large combinatorial problem. Using greedy algorithms and spectral proxies [3] to sidestep some of the computational bottlenecks, these methods enable to efficiently perform optimal sampling on graphs of size up to a few thousands nodes.

Loosening the tightness constraint on sampling sets, random graph sampling [5] shows that if $O(k\log{k})$ nodes are drawn with replacement according to a particular probability distribution $p^{*}$ that depends on the first $k$ eigenvectors of the Laplacian, then recovery is guaranteed with high probability. $p^{*}$ was shown to be closely related to leverage scores, an important concept in randomized numerical linear algebra [9]. Moreover, an efficient estimation of $p^{*}$, bypassing expensive spectral computations, enables to sample very large graphs; as was shown in [10], where experiments were performed on graphs with $N\simeq 10^{6}$ and $k\simeq 200$.

In this paper, we detail another form of random graph sampling, motivated by a simple experimental context. Consider a graph composed of two almost disconnected communities of equal size. Its first Fourier mode is constant (as always), and its second one is positive in one community, and negative in the other, while its absolute value is more or less constant on the whole graph. Therefore, a 2-bandlimited graph signal will approximately be equal to a constant in one community, and to another constant in the other: any sampling set composed of one node drawn from each community will enable to recover the whole signal. Unfortunately, random graph sampling as proposed in [5] samples independently with replacement: it has a fair chance of sampling only in one community, thereby missing half of the information. To avoid such situations, we propose to take advantage of determinantal point processes (DPP), random processes that introduce negative correlations to force diversity within the sampling set.

Contributions. If the first $k$ eigenvectors of the Laplacian are computable, we exhibit a DPP that always samples sets of size $k$ perfectly embedding $k$-bandlimited graph signals. Moreover, in the case the graph is too large to compute its first $k$ eigenvectors, and building upon recent results linking determinantal processes and random rooted spanning forests [11], we show to what extent a simple algorithm based on loop-erased random walks is an acceptable approximation.

II Background

Consider an undirected weighted graph $\mathcal{G}$ of $N$ nodes represented by its adjacency matrix $\mathsf{W}$, where $\mathsf{W}_{ij}=\mathsf{W}_{ji}\geqslant 0$ is the weight of the connection between nodes $i$ and $j$. Denote by $\mathsf{L}=\mathsf{D}-\mathsf{W}$ its Laplacian matrix, where $\mathsf{D}$ is the diagonal matrix with entries $\mathsf{D}_{ii}=\sum_{j}\mathsf{W}_{ij}$. $\mathsf{L}$ is real and symmetrical, therefore diagonalizable as $\mathsf{U}\mathsf{\Lambda}\mathsf{U}^{\intercal}$ with $\mathsf{U}=(\bm{u}_{1}|\ldots|\bm{u}_{N})\in\mathbb{R}^{N\times N}$ the graph Fourier basis and $\mathsf{\Lambda}=\text{diag}(\lambda_{1},\ldots,\lambda_{N})$ the diagonal matrix of the associated frequencies, that we sort: $0=\lambda_{1}\leqslant\lambda_{2}\leqslant\ldots\leqslant\lambda_{N}$. For any function $h:\mathbb{R}^{+}\rightarrow\mathbb{R}$, we write $h(\mathsf{\Lambda})=\text{diag}(h(\lambda_{1}),\ldots,h(\lambda_{N}))$. The graph filter associated to the frequency response $h$ reads $\mathsf{U}h(\mathsf{\Lambda})\mathsf{U}^{\intercal}$.

Given the frequency interpretation of the graph Fourier modes, one may define “smooth” signals as linear combinations of the first few low-frequency Fourier modes. Writing $\mathsf{U}_{k}=(\bm{u}_{1}|\ldots|\bm{u}_{k})\in\mathbb{R}^{N\times k}$, we have the formal definition:

Definition II.1** ($k$-bandlimited signal on $\mathcal{G}$).**

A signal $\bm{x}\in\mathbb{R}^{N}$ defined on the nodes of the graph $\mathcal{G}$ is $k$-bandlimited with $k\in\mathbb{N}^{*}$ if $\bm{x}\in{\rm span}(\mathsf{U}_{k})$, i.e., $\exists~{}\bm{\alpha}\in\mathbb{R}^{k}$ such that $\bm{x}=\mathsf{U}_{k}\bm{\alpha}$.

Sampling consists in selecting a subset $\mathcal{A}=(\omega_{1},\ldots,\omega_{m})$ of $m$ nodes of the graph and measuring the signal on it. To each possible sampling set, we associate a measurement matrix $\mathsf{M}=(\bm{\delta}_{\omega_{1}}|\bm{\delta}_{\omega_{2}}|\ldots|\bm{\delta}_{\omega_{m}})^{\intercal}\in\mathbb{R}^{m\times N}$ where $\bm{\delta}_{\omega_{i}}(j)=1$ if $j=\omega_{i}$, and 0 otherwise. Given a $k$-bandlimited signal $\bm{x}$, its (possibly noisy) measurement on $\mathcal{A}$ reads:

[TABLE]

where $\bm{n}$ models any kind of measurement noise and/or the fact that $\bm{x}$ may not exactly be in ${\rm span}(\mathsf{U}_{k})$. In this latter case, $\bm{x}=\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\bm{x}+\bm{e}$, where $\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\bm{x}$ is the closest signal to $\bm{x}$ that is in ${\rm span}(\mathsf{U}_{k})$; and $\mathsf{M}\bm{x}=\mathsf{M}\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\bm{x}+\bm{n}$ with $\bm{n}=\mathsf{M}\bm{e}$.

II-A If $\mathsf{U}_{k}$ is known

If $\mathsf{U}_{k}$ is known, the measurement on $\mathcal{A}$ reads $\bm{y}=\mathsf{M}\mathsf{U}_{k}\bm{\alpha}+\bm{n}$ and one recovers the signal by solving:

[TABLE]

where $(\mathsf{M}\mathsf{U}_{k})^{\dagger}\in\mathbb{R}^{k\times m}$ is the Moore-Penrose pseudo-inverse of $\mathsf{M}\mathsf{U}_{k}\in\mathbb{R}^{m\times k}$ and $\left\|.\right\|$ the Euclidean norm. Of course, perfect reconstruction is impossible if $m<k$, as the system would be underdetermined. Hereafter, we thus suppose $m\geqslant k$. Denote by $\sigma_{1}\leqslant\ldots\leqslant\sigma_{k}$ the singular values of $\mathsf{M}\mathsf{U}_{k}$.

Proposition II.1** (Classical result).**

If $\sigma_{1}>0$, then perfect reconstruction (up to the noise level) is possible.

Proof.

If $\sigma_{1}>0$, then $\mathsf{U}_{k}^{\intercal}\mathsf{M}^{\intercal}\mathsf{M}\mathsf{U}_{k}$ is invertible; and since $\bm{x}=\mathsf{U}_{k}\bm{\alpha}$, we have:

[TABLE]

In the following, in the cases where $\mathsf{U}_{k}$ is known, we fix $m=k$. Chen et al. [4] showed that a sample of size $k$ always exists s.t. $\sigma_{1}>0$. In fact, in general, many possible such subsets exist. Finding the optimal one is a matter of definition. Authors in [4, 3] propose two optimality definitions via noise robustness. A first option is to minimize the worst-case error, which translates to finding the subset that maximizes $\sigma_{1}^{2}$:

[TABLE]

A second option is to find the subset that minimizes the mean square error (assuming $\mathbb{E}(\bm{n}\bm{n}^{\intercal})$ is the identity) :

[TABLE]

where tr is the trace operator.

We argue that another possible definition of optimal subset is such that the determinant of $\mathsf{U}_{k}^{\intercal}\mathsf{M}^{\intercal}\mathsf{M}\mathsf{U}_{k}$ is maximal, i.e.:

[TABLE]

where MV stands for Maximum Volume, as the determinant is proportional to the volume spanned by the $k$ sampled lines of $\mathsf{U}_{k}$. This definition of optimality does not have a straight-forward noise robustness interpretation, even though it is related to the previous definitions (in particular, increasing $\sigma_{1}$ necessarily increases the determinant). Nevertheless, such volume sampling as it is sometimes called [12] has been proved optimal in the related problem of low-rank approximation [13], and is largely studied in diverse contexts [14].

In all three cases, finding the optimal subset is a very large combinatorial problem (it implies a search over all combinations of $k$ nodes): in practice, greedy algorithms such as Alg. 1 are used to find approximate solutions. In the following, we refer to $\hat{\mathcal{A}}^{\text{WCE}}$, $\hat{\mathcal{A}}^{\text{MSE}}$ and $\hat{\mathcal{A}}^{\text{MV}}$ the (deterministic) solutions obtained by greedy optimization of the above 3 objectives. We refer to $\hat{\mathcal{A}}^{\text{maxvol}}$ the solution to the maxvol algorithm proposed in [15], which is another approximation of $\mathcal{A}^{\text{MV}}$.

II-B If $\mathsf{U}_{k}$ is unknown

If $\mathsf{U}_{k}$ is unknown, one may use a proxy that penalizes high frequencies in the solution, and thus recover an approximation of the original signal by solving the regularized problem:

[TABLE]

where $\gamma>0$ is the regularization parameter and $r$ the power of the Laplacian that controls the strength of the high frequency penalization. By differentiating with respect to $\bm{z}$, one has:

[TABLE]

which may be solved by direct inversion if $N$ is not too large, or –thanks to the objective’s convexity– by iterative methods such as gradient descent.

II-C Reweighting in the case of random sampling

In random sampling, the sample $\mathcal{A}=(\omega_{1},\ldots,\omega_{m})$ is a random variable, and therefore so is the measurement $\bm{y}=\mathsf{M}\bm{x}$. Consider two signals $\bm{x}_{1}$ and $\bm{x}_{2}\neq\bm{x}_{1}$, and their measurements $\bm{y}_{1}$ and $\bm{y}_{2}$. If $\bm{y}_{1}=\bm{y}_{2}$, perfect reconstruction is impossible. To prevent this, one possible solution is to reweight the measurement such that the expected norm of the reweighted measurement equals the signal’s norm. For instance, in [5], $m$ nodes are drawn independently with replacement. At each draw, the probability to sample node $i$ is $p^{*}_{i}=\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|^{2}/k$. One has $\mathbb{E}\left(\left\|\mathsf{P}^{-1/2}\mathsf{M}\bm{x}\right\|^{2}\right)=\left\|\bm{x}\right\|^{2}$ if

[TABLE]

This guarantees that if $\left\|\bm{x}_{1}-\bm{x}_{2}\right\|>0$, then, for a large enough $m$, the reweighted measures will necessarily be distinct in the measurement space, i.e., $\left\|\mathsf{P}^{-1/2}(\bm{y}_{1}-\bm{y}_{2})\right\|>0$, ensuring that $\bm{x}_{1}$ and $\bm{x}_{2}$ have a chance of being recovered.

Recovery with known $\mathsf{U}_{k}$ (Eq. (2)) thus transforms into

[TABLE]

and recovery with unknown $\mathsf{U}_{k}$ (Eq. (6)) transforms into

[TABLE]

III Determinantal processes

Denote by $[N]$ the set of all subsets of $\{1,2,\ldots,N\}$.

Definition III.1** (Determinantal Point Process).**

Consider a point process, i.e., a process that randomly draws an element $\mathcal{A}\in[N]$. It is determinantal if, for every $\mathcal{S}\subseteq\mathcal{A}$,

[TABLE]

where $\mathsf{K}\in\mathbb{R}^{N\times N}$, a semi-definite positive matrix s.t. $0\preceq\mathsf{K}\preceq 1$, is called the marginal kernel; and $\mathsf{K}_{\mathcal{S}}$ is the restriction of $\mathsf{K}$ to the rows and columns indexed by the elements of $\mathcal{S}$.

III-A Sampling from a DPP and signal recovery

Given the constraints on $\mathsf{K}$, its eigendecomposition $\mathsf{K}=\sum_{i=1}^{N}\mu_{i}\bm{v}_{i}\bm{v}_{i}^{\intercal}$ always exists, and Alg. 2 provides a sample from the associated DPP [14]. Note that $|\mathcal{A}|$, the number of elements in $\mathcal{A}$, is distributed as the sum of $N$ Bernoulli trials of probability $\mu_{i}$ (see the “for loop” of Alg. 2). In particular:

[TABLE]

Consider $\mathcal{A}$, the random variable of a DPP with kernel $\mathsf{K}$; and $\mathsf{M}$ its associated measurement matrix. Denote by $\pi_{i}=\mathsf{K}_{ii}$ the marginal probability that $\{i\}\subseteq\mathcal{A}$ and write:

[TABLE]

Proposition III.1**.**

$\forall\bm{x}\in\mathbb{R}^{N}~{}~{}\mathbb{E}_{\mathcal{A}}\left(\left\|\mathsf{P}^{-1/2}\mathsf{M}\bm{x}\right\|^{2}\right)=\left\|\bm{x}\right\|^{2}.$**

Proof.

Let $\bm{x}\in\mathbb{R}^{N}$. In the following, we need the probability $\mathbb{P}(\mathcal{A})$ that $\mathcal{A}$ is sampled (instead of marginal probabilities). One can show that $\mathbb{P}(\mathcal{A})$ is proportional to the determinant of the restriction to $\mathcal{A}$ of a matrix called $L$-ensemble [14]. Note that: $\forall i~{}\sum_{\mathcal{A}\supset i}\mathbb{P}(\mathcal{A})=\mathbb{P}(\{i\}~{}\text{is sampled})=\pi_{i}$. Thus:

[TABLE]

This result is the first step to proove a restricted isometry property (RIP) for $\mathsf{P}^{-1/2}\mathsf{M}$ – currently work in progress. In this paper’s context, it shows how to reweight the measurement in the case of sampling from a DPP: recovery should be performed via optimization of (9) or (10) with $\mathsf{P}$ as in (12).

III-B If $\mathsf{U}_{k}$ is known: the ideal low-pass marginal kernel

Consider the DPP defined by the following marginal kernel:

[TABLE]

with $h_{\lambda_{k}}$ s.t. $h_{\lambda_{k}}(\lambda)=1$ if $\lambda\leqslant\lambda_{k}$ and [math] otherwise. In GSP words, $\mathsf{K}_{k}$ is the ideal low-pass with cutting frequency $\lambda_{k}$.

Proposition III.2**.**

A sample from the DPP with marginal kernel $\mathsf{K}_{k}$ defined as in (13) is of size $k$ and $\forall\bm{x}\in{\rm span}{(\mathsf{U}_{k})}$, the measurement $\bm{y}=\mathsf{M}\bm{x}+\bm{n}\in\mathbb{R}^{k}$ enables perfect reconstruction up to the noise level.

Proof.

$\mathsf{K}_{k}=\sum_{i=1}^{N}\mu_{i}\bm{u}_{i}\bm{u}_{i}^{\intercal}$, with $\mu_{1}=\ldots=\mu_{k}=1$ and $\mu_{k+1}=\ldots=\mu_{N}=0$. According to (11), $\mathbb{E}(|\mathcal{A}|)=k$ and $\text{Var}(|\mathcal{A}|)=0$, implying that any sample $\mathcal{A}$ has size $k$. Also, by definition of a DPP, the fact that $\mathcal{A}$ is sampled implies:

[TABLE]

implying $\sigma_{1}>0$. Using Prop. II.1 completes the proof. ∎

The most probable sample from this DPP is the sample $\mathcal{A}$ of size $k$ that maximizes $\text{det}(\mathsf{M}\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\mathsf{M}^{\intercal})$ i.e., $\mathcal{A}^{\text{MV}}$, the solution to the maximal volume optimisation problem (5).

III-C If $\mathsf{U}_{k}$ is unknown: Wilson’s marginal kernel

First, add a node $\Delta$ to the graph, that we will call “absorbing state” and that we connect to all existing $N$ nodes with a weight $q$. We will propagate random walks on this extended graph. Note that the probability to jump from any node $i\neq\Delta$:

• •

to any other node $j\neq\Delta$ is equal to $\mathsf{W}_{ij}/(\mathsf{D}_{ii}+q)$;

• •

to $\Delta$ is $q/(\mathsf{D}_{ii}+q)$.

Once the random walk reaches $\Delta$, it cannot escape from it. Wilson and Propp [16, 17] introduced Algorithm 3 in the context of spanning tree sampling. This Algorithm outputs a set $\mathcal{Y}$ of nodes by propagating loop-erased random walks.

Proposition III.3** (Corollary of Thm 3.4 in [11]).**

$\mathcal{Y}$* is a sample from a DPP with marginal kernel $\mathsf{K}_{q}=\mathsf{U}g_{q}(\Lambda)\mathsf{U}^{\intercal}$, with $g_{q}(\lambda)=\frac{q}{q+\lambda}$.*

Proof.

Given that the weight $q$ we added is finite and constant over all nodes, Thm 3.4 and Lemma 3.3 in [11] state that $\mathcal{Y}$ is distributed as a DPP with marginal kernel equal to $(\mathsf{L}+q\mathsf{I})^{-1}q\mathsf{I}$ where $\mathsf{I}$ is the identity matrix (note that $L$ in [11] refers to $-\mathsf{L}$ in our paper). Decomposing $\mathsf{L}$ as $\mathsf{U}\Lambda\mathsf{U}^{\intercal}$ finishes the proof. ∎

Prop. III.2 shows that sampling with $\mathsf{K}_{k}$ is optimal in the sense it enables perfect reconstruction with only $m=k$ measurements. Unfortunately, as $N$ or $k$ increase, computing $\mathsf{K}_{k}$ becomes prohibitive. Considering $g_{q}$ as an approximation of $h_{\lambda_{k}}$, one may see sampling with $\mathsf{K}_{q}$ as an approximation of optimal sampling with $\mathsf{K}_{k}$. Prop. III.3 has a beautiful consequence: one does not need to explicitly compute $\mathsf{K}_{q}$ (which would also be too expensive) to sample from it. In fact, Alg. 3 efficiently samples from $\mathsf{K}_{q}$. Thus, Alg. 3 approximates sampling from $\mathsf{K}_{k}$ without any spectral computation of $\mathsf{L}$.

IV Experiments

The Stochastic Block Model (SBM). We consider random community-structured graphs drawn from the SBM. We specifically look at graphs with $k$ communities of same size $N/k$. In the SBM, the probability of connection between any two nodes $i$ and $j$ is $q_{1}$ if they are in the same community, and $q_{2}$ otherwise. One can show that the average degree reads $c=q_{1}\left(\frac{N}{k}-1\right)+q_{2}\left(N-\frac{N}{k}\right)$. Thus, instead of providing the probabilities $(q_{1},q_{2})$, one may characterize a SBM by considering $(\epsilon=\frac{q_{2}}{q_{1}},c)$. The larger $\epsilon$, the fuzzier the community structure. In fact, authors in [18] show that above the critical value $\epsilon_{c}=(c-\sqrt{c})/(c+\sqrt{c}(k-1))$, community structure becomes undetectable in the large $N$ limit. In the following, $N=100$, $k=2$, $c=16$ and we let $\epsilon$ vary.

To generate a $k$-bandlimited signal, we draw $k$ realisations of a normal distribution with mean [math] and variance $1$, to build the vector $\bm{\alpha}\in\mathbb{R}^{k}$. Then, renormalize to make sure that $\left\|\bm{\alpha}\right\|^{2}=1$, and write $\bm{x}=\mathsf{U}_{k}\bm{\alpha}$. In our experiments, the measurement noise is Gaussian with mean 0 and $\sigma_{n}=10^{-4}$.

First series: $\mathsf{U}_{k}$ known. We compare DPP sampling from $\mathsf{K}_{k}$ followed by recovery using (9); to deterministic sampling with $\hat{\mathcal{A}}^{\text{WCE}}$, $\hat{\mathcal{A}}^{\text{MSE}}$, $\hat{\mathcal{A}}^{\text{MV}}$ and $\hat{\mathcal{A}}^{\text{maxvol}}$ presented in Sec. II-A followed by recovery using (2). We show in Fig. 1a) the recovery performance with respect to $\epsilon/\epsilon_{c}$: all greedy methods perform similarly, and on average slightly outperform the DPP-based method, especially as the graphs become less structured.

Second series: $\mathsf{U}_{k}$ unknown. We compare negatively correlated sampling from $\mathsf{K}_{q}$ (using Wilson’s algorithm) followed by recovery using (10) with $\mathsf{P}$ as in (12); to the uncorrelated random sampling [5] where nodes are sampled independently with replacement from $p^{*}$ with $p^{*}_{i}=\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|^{2}/k$, followed by recovery using (10) with $\mathsf{P}$ as in (8). Several comments are in order. i) the number of samples from Wilson’s algorithm is not known in advance, but we know its expected value $\mathbb{E}(|\mathcal{A}|)=\sum_{i}q/(q+\lambda_{i})$. In order to explore the behavior around the critical number of measurements $m=k$, one needs to choose $q$ s.t. $\sum_{i}q/(q+\lambda_{i})\simeq k$, but without computing the $\lambda_{i}$! To do so, we follow the proposition of [11] (see discussion around Eq. (4.14)) based on a few runs of Alg. 3. Once a sample of size approximately equal to $k$ is exhibited, and for fair comparison, one uses the same number of nodes for the uncorrelated sampling. ii) exact computation of $p^{*}$ for the uncorrelated sampling requires $\mathsf{U}_{k}$. We follow the efficient Alg. 1 of [5] to approximate $p^{*}$. iii) recovery using $\mathsf{P}$ as in (12) requires to know $\pi_{i}=\sum_{j=1}^{N}g_{q}(\lambda_{j})\bm{u}_{j}(i)^{2}$ for each node $i$. To estimate $\pi_{i}$ without spectral decomposition of $\mathsf{L}$, an elegant solution is via fast graph filtering of random signals, as in Alg. 4. Building upon the Johnson-Lindenstrauss lemma, one can show that $\hat{\pi}_{i}$ concentrates around $\pi_{i}$ for $n=O(\log{N})$ (see a similar proof in [19]). We fix $n=20\log{N}$ and $d=30$.

The reconstruction’s parameters are $\gamma$ and $r$. Following [5], $r$ is fixed to $4$. Fig. 1b) shows that the reconstruction performance saturates at $\gamma=10^{-5}$. With these values of $r$ and $\gamma$, we compare DPP sampling with $\mathsf{K}_{q}$ vs uncorrelated sampling in Fig. 1c) with respect to the number of measurements, for different values of $\epsilon$. The lower $\epsilon$, i.e. the stronger the community structure, the better DPP sampling is compared to uncorrelated sampling. In terms of computation time, a precise comparison of both methods versus the different parameters is out of this paper’s scope. Nevertheless, to give an idea of the method’s scalability, for a SBM with $N=10^{5}$ (resp. $10^{6}$), $k=2$, $c=16$ and $\epsilon=\epsilon_{c}/5$, and with $q=5\cdot 10^{-4}$, Wilson’s algorithm outputs in average 5 (resp. 36) samples in a mean time of 7 (resp. 90) seconds, using Python on a laptop.

V Conclusion

We first show that sampling from a DPP with marginal kernel $\mathsf{K}_{k}=\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}$ i.e. the projector onto the first $k$ graph Fourier modes, outputs a subset of size $k$ that enables perfect reconstruction. The robustness to noise of the reconstruction is comparable to state-of-the-art greedy sampling algorithms. Moreover, in the important and fairly common case where the first $k$ eigenvectors of the Laplacian are not computable, we show that Wilson’s algorithm (Alg. 3) may be leveraged as an efficient graph sampling scheme that approximates sampling from $\mathsf{K}_{k}$. Preliminary experiments on the SBM suggest that in the cases where the number of samples $m$ needs to stay close to $k$ and where the community structure of the graph is strong, sampling with Wilson’s algorithm enables a better reconstruction than the state-of-the-art uncorrelated random sampling. Further investigation, both theoretical and experimental, is necessary to better specify which families of graphs are favorable to DPP sampling, and which are not.