\'Echantillonnage de signaux sur graphes via des processus d\'eterminantaux

TL;DR

This paper introduces a novel sampling strategy for k-bandlimited graph signals using determinantal point processes, enabling efficient reconstruction with fewer samples without full Laplacian diagonalization.

Contribution

It proposes a new method based on determinantal point processes for sampling graph signals, avoiding costly Laplacian diagonalization and improving sampling efficiency.

Findings

Reconstruction with only O(k) samples is possible.

The proposed sampling algorithm is faster than existing methods by an order of k.

The method avoids partial diagonalization of the graph Laplacian.

Abstract

We consider the problem of sampling k-bandlimited graph signals, ie, linear combinations of the first k graph Fourier modes. We know that a set of k nodes embedding all k-bandlimited signals always exists, thereby enabling their perfect reconstruction after sampling. Unfortunately, to exhibit such a set, one needs to partially diagonalize the graph Laplacian, which becomes prohibitive at large scale. We propose a novel strategy based on determinantal point processes that side-steps partial diagonalisation and enables reconstruction with only O(k) samples. While doing so, we exhibit a new general algorithm to sample determinantal process, faster than the state-of-the-art algorithm by an order k.

Figures

Equations

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

$${\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}.$$

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

$$\text{det}\left({\mathsf{K}}_{{\mathcal{S}}_{n-1}\cup \{i\}}\right)$$

$$=p_{n-1}(i)\text{det}\left({\mathsf{K}}_{{\mathcal{S}}_{n-1}}\right).\text{vspace}-0.3cm$$

$$\mathbb{P}(\mathcal{A})$$

$$\forall l\in [1,m],\sum _{i=1}^N{p}_{l-1}(i)=\sum _{i=1}^N{p}_0(i)-\sum _{i=1}^N{\mathsf{K}}_{{\mathcal{S}}_{l-1},i}^{⊺}{\mathsf{K}}_{{\mathcal{S}}_{l-1}}^{-1}{\mathsf{K}}_{{\mathcal{S}}_{l-1},i}\text{vspace}-0.3cm$$

$$\sum _{i=1}^N$$

$$=\text{Tr}\left(({\mathsf{MXX}}^{⊺}{\mathsf{M}}^{⊺}{)}^{-1}{\mathsf{MXX}}^{⊺}{\mathsf{XX}}^{⊺}\mathsf{M}\right)=\text{Tr}({\mathsf{I}}_{l-1})=l-1,\text{vspace}-0.3cm$$

$$\mathbb{P}(\mathcal{A})=\frac{1}{Z}\text{det}({\mathsf{K}}_{\mathcal{A}})\text{leavevmodenobreak avecleavevmodenobreak }Z=\prod _{l=1}^{m}d-l+1,\text{vspace}-0.3cm$$

$$\forall n\forall i\sum _{l=1}^n{f}_l(i{)}^2={\mathsf{K}}_{{\mathcal{S}}_n,i}^{⊺}{\mathsf{K}}_{{\mathcal{S}}_n}^{-1}{\mathsf{K}}_{{\mathcal{S}}_n,i}.\text{vspace}-0.3cm$$

$$\forall n\forall i,j\sum _{l=1}^n{f}_l(i)f_l(j)={\mathsf{K}}_{{\mathcal{S}}_n,i}^{⊺}{\mathsf{K}}_{{\mathcal{S}}_n}^{-1}{\mathsf{K}}_{{\mathcal{S}}_n,j}.\text{vspace}-0.3cm$$

$$\forall i,j{\mathsf{K}}_{{\mathcal{S}}_n,i}^{⊺}{\mathsf{K}}_{{\mathcal{S}}_n}^{-1}{\mathsf{K}}_{{\mathcal{S}}_n,j}={\mathsf{K}}_{s_1,i}{\mathsf{K}}_{s_1,j}/{\mathsf{K}}_{s_1,s_1}=f_1(i)f_1(j).\text{vspace}-0.3cm$$

$${\mathsf{K}}_{{\mathcal{S}}_n,i}^{⊺}{\mathsf{K}}_{{\mathcal{S}}_n}^{-1}{\mathsf{K}}_{{\mathcal{S}}_n,j}$$

$$\forall i,j{f}_n(i)f_n(j)=\frac{z_n(i)z_n(j)}{z_n(s_n)}.\text{vspace}-0.3cm$$

$$\forall i{f}_n(i)=\frac{{\mathsf{K}}_{s_n,i}-{\sum }_{l=1}^{n-1}{f}_l(i)f_l(s_n)}{\sqrt{{\mathsf{K}}_{s_n}-{\sum }_{l=1}^{n-1}{f}_l(s_n{)}^2}}.\text{vspace}-0.3cm$$

$$\forall i{f}_n(i)=\frac{z_n(i)}{\sqrt{z_n(s_n)}},\text{vspace}-0.3cm$$

$$\forall \lambda \in [0,{\lambda }_N]{\tilde{h}}_k(\lambda )=\sum _{l=1}^r{\alpha }_l{\lambda }^l\simeq h_k(\lambda ).$$

$$\text{On a : }{\mathsf{K}}_k\simeq \mathsf{U}{\tilde{h}}_k(\mathsf{\Lambda }){\mathsf{U}}^{⊺}=\sum _{l=1}^r{\alpha }_l\mathsf{U}{\Lambda }^l{\mathsf{U}}^{⊺}=\sum _{l=1}^r{\alpha }_l{\mathsf{L}}^l={\tilde{h}}_k(\mathsf{L}).\text{vspace}-0.1cm$$

$$\mathbb{E}$$

$$={\mathbf{\delta }}_i^{⊺}({\tilde{h}}_k(\mathsf{L}){)}^2{\mathbf{\delta }}_i\simeq {\mathbf{\delta }}_i^{⊺}{\mathsf{K}}_k^2{\mathbf{\delta }}_i={\mathsf{K}}_k(i,i)\text{vspace}-0.3cm$$

$$p_0(i)\simeq {\Vert {\mathbf{\delta }}_i^{⊺}{\tilde{h}}_k(\mathsf{L})\mathsf{R}\Vert }^2,\text{vspace}-0.3cm$$

Full text

[TABLE]

À partir de (Démonstration.), et sachant que i) $\mathsf{K}$ est SDP: $\forall\mathcal{S},\,\text{det}(\mathsf{K}_{\mathcal{S}})\geqslant 0$, ii) $d\geqslant m$, on peut montrer que $p_{n}(i)\geqslant 0$ et $\sum_{i}{p_{n}}(i)\neq 0$, c’est-à-dire qu’à chaque itération de la boucle, la probabilité $\mathbb{P}(s)$ est bien définie. La boucle étant répétée $m$ fois, la sortie de l’algorithme, notée $\mathcal{A}$, est nécessairement de taille $m$, en accord avec le point i) de la Déf. 2. Enfin, montrons que $\mathbb{P}(\mathcal{A})$ est conforme au point ii). Par construction de $\mathcal{A}$ :

[TABLE]

Or, en écrivant (Démonstration.) pour $i=s_{n}$, et en itérant, on obtient : $\prod_{l=1}^{m}p_{l-1}(s_{l})=\text{det}(\mathsf{K}_{\mathcal{A}})$. Reste à montrer que le dénominateur de (6) ne dépend pas des échantillons choisis. C’est là où l’hypothèse d’un noyau de projection est essentielle. On a :

[TABLE]

On a $\sum_{i=1}^{N}p_{0}(i)=\text{Tr}(\mathsf{XX}^{\intercal})=\text{Tr}(\mathsf{X}^{\intercal}\mathsf{X})=d$. De plus, soit $\mathsf{M}$ la matrice de mesure associée à $\mathcal{S}_{l-1}$:

[TABLE]

par invariance de la trace aux permutations circulaires. Ainsi :

[TABLE]

ce qui termine la preuve. ∎

$\mathcal{A}^{\text{\footnotesize MV}}$, l’ensemble à approcher, est l’ensemble le plus probable du $m$-PPD (avec $m=k$) à noyau $\mathsf{K}_{k}=\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}$. Si nous avons les ressources pour calculer $\mathsf{U}_{k}$, une première stratégie d’échantillonnage est donc l’Alg. 1 appliqué à $\mathsf{K}_{k}$. Dans le cas contraire, nous proposons dans la suite un nouvel algorithme d’échantillonnage de $m$-PPD, de complexité inférieure d’un facteur $m$, qui permet d’appliquer des techniques d’approximation polynomiale, évitant ainsi toute étape de diagonalisation.

4  Approximation via des filtres sur graphe

4.1  Réécriture de l’algorithme d’échantillonnage

Proposition 2.

• L’algorithme 2 est équivalent à l’algorithme 1 : il échantillonne aussi un $m$-PPD à noyau de projection $\mathsf{K}$.*

[TABLE]

Nous allons montrer plus généralement que :

[TABLE]

Pour ce faire, nous proposons une récurrence.

Initialisation. C’est vrai pour $n=1$, où $\mathcal{S}_{n}$ est réduit à $\{s_{1}\}$ :

[TABLE]

Hypothèse. Supposons que (9) est vraie à l’étape $n-1$.

Récurrence. Montrons qu’elle est également vraie à l’étape $n$. En utilisant l’identité de Woodbury sur $\mathsf{K}_{\mathcal{S}_{n}}^{-1}$, on montre que:

[TABLE]

où $z_{n}(i)=\mathsf{K}_{s_{n},i}-\mathsf{K}_{\mathcal{S}_{n-1},s_{n}}^{\intercal}\mathsf{K}_{\mathcal{S}_{n-1}}^{-1}\mathsf{K}_{\mathcal{S}_{n-1},i}$. En remplaçant $\mathsf{K}_{\mathcal{S}_{n-1},i}^{\intercal}\mathsf{K}_{\mathcal{S}_{n-1}}^{-1}\mathsf{K}_{\mathcal{S}_{n-1},j}$ par $\sum_{l=1}^{n-1}f_{l}(i)f_{l}(j)$ grâce à l’hypothèse, il nous reste à montrer que :

[TABLE]

Or, par construction dans l’Algorithme 2, $f_{n}(i)$ s’écrit :

[TABLE]

En utilisant une seconde fois l’hypothèse (2), on montre que :

[TABLE]

ce qui prouve (10) et termine la preuve. ∎

L’algorithme 2 est un algorithme général pour échantillonner un $m$-PPD à noyau de projection. Sa complexité est en $O(Nm^{2})$, alors que la complexité de l’algorithme 1 est en $O(Nm^{3})$. Notons que des idées similaires pour réduire la complexité d’un facteur $m$ existent dans la littérature (par exemple dans [9]) mais sous des formes un peu cachées, et, à notre connaissance, n’ont jamais été vraiment explicitées dans le cas discret.

4.2  Approcher l’échantillonnage d’un PPD

Nous allons voir que la manière dont l’algorithme 2 fait appel à $\mathsf{K}$ permet de tirer profit d’approximations polynomiales usuellement utilisées lors d’opérations de filtrage sur graphe. Rappelons que nous considérons le $m$-PPD associé au noyau $\mathsf{K}_{k}$. Or : $\mathsf{K}_{k}=\mathsf{U}h_{k}(\mathsf{\Lambda})\mathsf{U}^{\intercal},$ où $h_{k}(\mathsf{\Lambda})=\text{diag}(h_{k}(\lambda_{1}),\ldots,h_{k}(\lambda_{N}))$ et $h_{k}(\lambda)$ est tel que $h_{k}(\lambda)=1$ si $\lambda\leqslant\lambda_{k}$ et $h_{k}(\lambda)=0$ sinon. En traitement du signal sur graphes, $\mathsf{K}_{k}$ est un filtre passe-bas idéal de fréquence de coupure $\lambda_{k}$.

Approximation polynomiale. [7] Considérons le polynôme de Tchebychev $\tilde{h}_{k}$ de degré $r$ qui approche au mieux $h_{k}$:

[TABLE]

[TABLE]

Filtrage rapide sur graphe. On ne calcule jamais explicitement $\tilde{h}_{k}(\mathsf{L})$ qui est en général dense de taille $N\times N$. En revanche, étant donné un signal $\bm{x}$ défini sur le graphe, le signal filtré par $h_{k}$, $\mathsf{K}_{k}\bm{x}$, est approché par $\tilde{h}_{k}(\mathsf{L})\bm{x}=\sum_{l=1}^{r}\alpha_{l}\mathsf{L}^{l}\bm{x}$, qui se calcule via $r$ multiplications matrice-vecteur si bien que le nombre d’opérations nécessaires pour filtrer un signal est de l’ordre de $r|E|$ où $|E|$ est le nombre de liens du graphe.

Estimation de la diagonale de $\mathsf{K}_{k}$. Soit $\mathsf{R}\in\mathbb{R}^{N\times n}$ une matrice contenant $n$ signaux aléatoires gaussiens de moyenne nulle et de variance $1/n$. Notons que:

[TABLE]

La $i^{\text{\footnotesize\`{e}me}}$ valeur de la diagonale de $\mathsf{K}_{k}$, $p_{0}(i)$, est donc approchée par la norme de la $i^{\text{\footnotesize\`{e}me}}$ ligne de $\tilde{h}_{k}(\mathsf{L})\mathsf{R}$, i.e. :

[TABLE]

et, via le lemme de Johnson-Lindenstrauss, on montre qu’un nombre de signaux aléatoires $n=O(\log{N})$ est suffisant pour avoir une estimation convenable avec haute probabilité [10].

Algorithme d’approximation. On souhaite approcher l’algorithme 2 avec entrée $\mathsf{K}_{k}$ sans avoir à calculer $\mathsf{U}_{k}$. Pour ce faire, on suit l’algorithme 3 en lui donnant comme entrée111pour avoir $h_{k}(\lambda)$, il faut $\lambda_{k}$, que nous estimons en suivant [4, Sec. 4.2]. $\mathsf{L}$, $h_{k}(\lambda)$, $r=50$ et $m$ le nombre d’échantillons souhaité. Au lieu d’accéder directement à la diagonale d’un noyau $\mathsf{K}$ connu, on l’estime à l’aide de l’équation (15). Puis, au lieu d’accéder aux colonnes $\bm{k}_{s}$ directement, on les estime via le filtrage rapide associé: $\bm{k}_{s}\simeq\tilde{h}_{k}(\mathsf{L})\bm{\delta}_{s}$. Au vu des erreurs d’approximation successives, l’écart entre les $\{\bm{f}_{l}\}$ calculés au cours des algorithmes 2 et 3 ne cesse de s’accroître au fur et à mesure de la boucle. Pour stabiliser l’algorithme 3 : i) le nouvel échantillon est celui qui maximise $\bm{p}$ (mise à part l’estimation de la diagonale de $\mathsf{K}$, l’algorithme n’est donc plus aléatoire); ii) et on s’autorise $m\geqslant k$ (ce qui n’a pas de sens pour l’algorithme 2 appliqué à $\mathsf{K}_{k}$ qui est de rang $k$). Aussi, en pratique, on fait quelques modifications mineures : i) après la mise à jour de $\bm{p}$, on force $p(s_{i})=0$ pour tous les échantillons $s_{i}$ déjà choisis; ii) lors de la normalisation de $\bm{f}_{n}$, si $\bm{f}_{n}(s_{n})\leqslant 0$ à cause de l’approximation, alors on normalise $\bm{f}_{n}$ par $\sqrt{\left\|\bm{f}_{n}\right\|/N}$. Ces choix sont cependant arbitraires : une comparaison rigoureuse des différentes possibilités de stabilisation fera l’objet de travaux futurs.

5  Simulations

Trois graphes. Nous considérons le graphe routier du Minnesota ($N=2642$), le graphe correspondant au maillage 3D d’un objet (graphe “Bunny” avec $N=2503$), et une réalisation d’un graphe aléatoire avec communautés issue du modèle stochastique par blocs (MSB) avec $N=1000$, 10 communautés de même taille et de paramètre de structure $\epsilon=\epsilon_{c}/4$, où $\epsilon_{c}$ correspond à la limite où les structures en communautés ne sont plus détectables; et plus $\epsilon$ est petit, plus la structure en communautés est forte. Voir par exemple[5, sec. 4] pour des détails sur $\epsilon$; [4, Fig. 1] pour des illustrations des deux premiers graphes.

Pour générer les signaux à bande limitée. Dans tous ces exemples, nous choisissons $k=10$. Chaque signal $\bm{x}$ à bande limitée $k$ est généré en calculant $\mathsf{U}_{k}$, puis en calculant $\bm{x}=\mathsf{U}_{k}\bm{\alpha}$ où $\bm{\alpha}\in\mathbb{R}^{k}$ est tiré aléatoirement selon une gaussienne $\mathcal{G}(0,1)$, puis en normalisant $\bm{x}$ pour qu’il soit de norme 1.

Les différentes méthodes d’échantillonnage comparées. On compare 5 méthodes d’échantillonnage : i) l’échantillonnage iid uniforme sans remise, ii) l’échantillonnage iid selon $\bm{p}=\text{diag}(\mathsf{K}_{k})$ sans remise (comme dans [4]), iii) la même stratégie mais avec une distribution $\bm{p}\simeq\text{diag}(\mathsf{K}_{k})$ estimée via (15), iv) l’échantillonnage dépendant en suivant l’algorithme 3, et v) l’échantillonnage selon le PPD idéal, c’est-à-dire en suivant l’algorithme 1 avec la légère modification qu’à chaque itération de la boucle on tire l’échantillon qui maximise la probabilité (cela donne de meilleures performances de reconstruction dans les cas testés). Les quatre premières méthodes peuvent échantillonner un nombre quelconque de nœuds $m$; la dernière échantillonne nécessairement $k$ nœuds.

Le bruit de mesure $\bm{n}\in\mathbb{R}^{m}$ est tiré selon une Gaussienne de moyenne nulle et d’écart-type $\sigma$, et est ajouté à la mesure du signal, comme dans (1). On fixe $\sigma=10^{-3}$.

Reconstruction des signaux échantillonnés. Quelle que soit la méthode choisie, on reconstruit les signaux avec (2), c’est-à-dire avec la connaissance de $\mathsf{U}_{k}$. Il existe des moyens de reconstruction qui ne nécessitent pas $\mathsf{U}_{k}$ (voir [4]); mais nous souhaitons ici comparer exclusivement l’échantillonnage.

Discussion des résultats. La Fig. 1 montre que notre méthode approche la performance du PPD à noyau $\mathsf{K}_{k}$ plus rapidement que les autres, dès $m=O(k)$. Pour le graphe du Minnesota, à petit $m$, notre performance est moindre. Ceci dit, nous retrouvons une meilleure performance pour d’autres choix de $k$. De nombreuses expériences ont été menées sur d’autres réalisations du MSB, avec d’autres valeurs de $k$ (différentes de $k=10$ choisie ici), $\epsilon$ et $\sigma$, et des communautés de taille hétérogène: nous trouvons toujours des comportements similaires.

6  Conclusion

Nous proposons l’algorithme 2, une réécriture moins coûteuse de l’algorithme 1 d’échantillonnage de $m$-DPP à noyau de projection. Dans les cas où le noyau est de la forme $\mathsf{K}=h_{k}(\mathsf{L})$ où $\mathsf{L}$ est une matrice diagonalisable (de préférence parcimonieuse), nous proposons l’algorithme 3 qui, via des approximations polynomiales, permet d’approcher l’échantillonnage du $m$-PPD associé à $\mathsf{K}$ sans jamais calculer explicitement $\mathsf{K}$, potentiellement dense, et en évitant tout calcul de diagonalisation de $\mathsf{L}$. Nous illustrons avec succès l’intérêt de cet algorithme pour l’échantillonnage de signaux sur graphe à bande limitée.

Remerciements. Ce travail a été en partie soutenu par le LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01), l’ANR GenGP (ANR-16-CE23-0008), le LIA CNRS/Melbourne Univ Geodesic.

Références

• [1]

A. Çivril, M. Magdon-Ismail. On selecting a maximum volume sub-matrix of a matrix and related problems. Theoretical Computer Science, no. 47, vol. 410, p. 4801–4811, 2009.

• [2]

S. Chen, R. Varma, A. Sandryhaila, J. Kovacevic. Discrete Signal Processing on Graphs: Sampling Theory. IEEE Transactions on Signal Processing, no. 24, vol. 63, p. 6510–6523, 2015.

• [3]

A. Anis, A. Gadde, A. Ortega. Efficient Sampling Set Selection for Bandlimited Graph Signals Using Graph Spectral Proxies. IEEE Transactions on Signal Processing, no. 14, vol. 64, p. 3775–3789, 2016.

• [4]

G. Puy, N. Tremblay, R. Gribonval, P. Vandergheynst. Random sampling of bandlimited signals on graphs. Applied and Computational Harmonic Analysis, in press, 2016.

• [5]

N. Tremblay, P-0. Amblard, S. Barthelmé. Graph sampling with determinantal processes. arXiv preprint, 1703.01594, 2017.

• [6]

A. Kulesza, B. Taskar, Determinantal point processes for machine learning. Found. and Trends in Mach. Learn., no. 2, vol. 5, p. 123–286, 2012.

• [7]

D. Shuman, S. Narang, P. Frossard, A. Ortega, P. Vandergheynst, The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains. IEEE Signal Processing Magazine, no. 3, vol. 30, p. 83–98, 2013.

• [8]

R. Bardenet, A. Hardy, Monte Carlo with determinantal point processes. arXiv preprint, 1605.00361, 2016.

• [9]

F. Lavancier, J. Moller, E. Rubak, Determinantal point process models and statistical inference. Journal of the Royal Statistical Society: Series B (Statistical Methodology), no. 4, vol. 77, p. 853–877, 2015

• [10]

N. Tremblay, G. Puy, R. Gribonval, P. Vandergheynst. Accelerated spectral clustering using graph filtering of random signals. ICASSP, 2016.