Stationary signal processing on graphs

Introduction

Stationarity is a traditional hypothesis in signal processing used to represent a special type of statistical relationship between samples of a temporal signal. The most commonly used is wide-sense stationarity, which assumes that the first two statistical moments are invariant under translation, or equivalently that the correlation between two samples depends only on their time difference. Stationarity is a corner stone of many signal analysis methods. The expected frequency content of stationary signals, called Power Spectral Density (PSD), provides an essential source of information used to build signal models, generate realistic surrogate data or perform predictions. In Figure 1, we present an example of a stationary process (blue curve) and two predictions (red and green curves). As the blue signal is a realization of a stationary process, the red curve is more probable than the green one because it respects the frequency content of the observed signal.

Classical stationarity is a statement of statistical regularity under arbitrary translations and thus requires a regular structure (often "time"). However many signals do not live on such a regular structure. For instance, imagine that instead of having one sensor returning a temporal signal, we have multiple sensors living in a two-dimensional space, each of which delivers only one value. In this case (see Figure 2 left), the signal support is no longer regular. Since there exists an underlying continuum in this example (2D space), one could assume the existence of a 2D stationary field and use Kriging to interpolate observations to arbitrary locations, thus generalizing stationarity for a regular domain but irregularly spaced samples.

On the contrary, the goal of this contribution is to generalize stationarity for an irregular domain that is represented by a graph, without resorting to any underlying regular continuum. Graphs are convenient for this task as they are able to capture complicated relations between variables. In this work, a graph is composed of vertices connected by weighted undirected edges and signals are now scalar values observed at the vertices of the graph. Our approach is to use a weak notion of translation invariance, define on a graph, that captures the structure (if any) of the data. Whereas classical stationarity means correlations are computed by translating the auto-correlation function, here correlations are given by localizing a common graph kernel, which is a generalized notion of translation as detailed in Section 2.2.

Figure 2 (left) presents an example of random multivariate variable living in a 2-dimensional space. Seen as scattered samples of an underlying 2D stochastic function, one would (rightly) conclude it is not stationary. However, under closer inspection, the observed values look stationary within the spiral-like structure depicted by the graph in Figure 2 (right). The traditional Kriging interpolation technique would ignore this underlying structure and conclude that there are always rapid two dimensional variations in the underlying continuum space. This problem does not occur in the graph case, where the statistical relationships inside the data follow the graph edges resulting in this example in signals oscillating smoothly over the graph.

A typical example of a stationary signal on a graph would be the result of a survey performed by the users of a social network. If there is a relationship between a user’s answers and those of his neighbors, this relationship is expected to be constant among all users. Using stationarity on the graph, we could predict the most probable answer for users that never took the survey.

We use spectral graph theory to extend the notion of stationarity to a broader class of signals. Leveraging the graph localization operator, we establish the theoretical basis of this extension in Section 3. We show that the resulting notion of stationarity is equivalent to the proposition of Girault [2, Definition 16], although the latter is not defined in terms of a localization operator. Localization is a very desirable feature, since it naturally expresses the scale at which samples are strongly correlated.

Since our framework depends on the power spectral density (PSD), we generalize the Welch method in Section 4 and obtain a scalable and robust way to estimate the PSD. It improves largely the covariance estimation when the number of signals is limited.

Based on the generalization of Wiener filters, we propose a new regularization term for graph signal optimization instead of the traditional Dirichlet prior, that depends on the noise level and on the PSD of the signal. The new optimization scheme presented in Section 5 has three main advantages: 1) it allows to deal with an arbitrary regularization parameter, 2) it adapts to the data optimally as we prove that the optimization model is a Maximum A Posteriori (MAP) estimator, and 3) it is more scalable and robust than a traditional Gaussian estimator.

Finally, in Section 6, we show experimentally that common datasets such as USPS follow our stationarity assumption. In section 7, we exploit this fact to perform missing data imputation and we show how stationarity improves over classical graph models and Gaussian MAP estimator.

2 Related work

Graphs have been used for regularization in data applications for more than a decade and two of the most used models will be presented in Section A. The idea of graph filtering was hinted at by the machine learning community but developed for the spectral graph wavelets proposed by Hammond et al. and extended by Shuman et al. in . While in most cases, graph filtering is based on the graph Laplacian, Moura et al. have suggested to use the adjacency matrix instead.

We note that a probabilistic model using Gaussian random fields has been proposed in . In this model, signals are automatically graph stationary with an imposed covariance matrix. Our model differentiates itself from these contributions because it is based on a much less restrictive hypothesis and uses the point of view of stationarity. A detailed explanation is given at the end of Section 3.

Finally, stationarity on graphs has been recently proposed in by Girault et al. These contributions use a different translation operator, promoting energy preservation over localization. While seemingly different, we show that our approach and Girault’s result in the same graph spectral characterization of stationary signals. Girault et al have also shown that using the Laplacian as a regularizer in a de-noising problem (Tikhonov) is equivalent to applying a Wiener filter adapted to a precise class of graph signals. In [2, pp 100], an expression of graph Wiener filter can be found.

After the publication of the first version of this contribution, additional work was done on the topic. First some PSD estimation methods were proposed in . Then stationarity has been extended to time evolving signals on graphs in .

Background theory

The most fundamental operator in graph signal-processing is the (combinatorial) graph Laplacian, defined as: $L=D-W,$ where $D$ is the diagonal degree matrix ($D[i,i]=d[i]$).

Spectral theory

where $U^{*}$ denotes the transposed conjugate of $U$. The graph Fourier transform is written $\hat{f}=U^{*}f$ and its inverse $f=U\hat{f}$. This Graph Fourier Transform possesses interesting properties further studied in . Note that the graph Fourier transform is equivalent to the Discrete Fourier transform for cyclic graphs. The detailed computation for the "ring" can be found in [22, pp 136-137].

Graph convolutive filters

2 Localization operator

As most graphs do not possess a regular structure, it is not possible to translate a signal around the vertex set with an intuitive shift. As stationarity is an invariance with respect to translation, we need to address this issue first. A solution is present in [10, Equation 26], where Shuman et. al. define the generalized translation for graphs as a convolution with a Kroneker delta. The convolution $\ast$ is defined as the element-wise multiplication in the spectral domain leading to the following generalized translation definition:

Naturally, the generalized translation operator does not perform what we would intuitively expect from it, i.e it does not shift a signal $s$ from node $n$ to node $i$ as this graph may not be shift-invariant. Instead when $\hat{s}$ changes smoothly across the frequencies (more details later on), then $T_{i}s$ is localized around node $i$, while $s$ is in general not localized at a particular node or set of nodes.

In order to avoid this issue, we define the localization operator as follows

Let us now clarify how generalized translation and localization are linked. The main difference between these two operators is the domain on which they are applied. Whereas, the translation operator acts on a discrete signal defined in the time or the vertex domain, the localization operator requires a continuous kernel or alternatively a discrete signal in the spectral domain. Both return a signal in the time or the vertex domain. To summarize, the localization operator can be seen as computing the inverse Fourier transform first and then translating the signal. It is an operator that takes a filter from the spectral domain and localizes it at a given node $i$ while adapting it to the graph structure.

In the classical periodic case ("ring" graph), localization is strongly connected to translation as the localized kernels are translated versions of each other:

3 Stationarity for temporal signals

A signal is Time Wide-Sense Stationary (WSS) if its first two statistical moments are invariant under translation, i.e:

where $\eta_{{\bf x}}$ is called the autocorrelation function of ${\bf x}$.

Note that using (3), the autocorrelation can be written in terms of the localization operator:

where $j=\sqrt{-1}$. As a consequence, when a signal is convolved with a filter $\check{h}$, its PSD is multiplied by the energy of the convolution kernel: for ${\bf y}=\check{h}\ast{\bf x}$, we have

When generalizing these concepts to graphs, the underlying structure for stationarity will no longer be time, but graph vertices.

Stationarity of graph signals

We now generalize stationarity to graph signals. While we define stationarity through the localization operator, Girault uses an isometric translation operator instead. That proposition is briefly described in Section 3.2, where we also show the equivalence of both definitions.

As explained in Section 2.2, the localization operator adapts a kernel to the graph structure. As a result, our idea is to use the localization operator to adapt the correlation between the samples to the graph structure. This results in a localized version of the correlation function, whose properties can then be studied via the associated kernel.

A stochastic graph signal ${\bf x}$ defined on the vertices of a graph $\mathcal{G}$ is called Graph Wide-Sense (or second order) Stationary (GWSS), if and only if it satisfies the following properties:

its covariance is the result of localizing a graph kernel:

The first part of the above definition is equivalent to the first property of time WSS signals. The requirement for the second moment is a natural generalization where we are imposing an invariance with respect to the localization operator instead of the translation. It is a generalization of Definiton 2 using (4). In simple words, the covariance is assumed to be driven by a global kernel (filter) $\gamma_{{\bf x}}$. The localization operator adapts this kernel to the local structure of the graph and provides the correlation between the vertices. Additionally, Definition 3 implies that the spectral components of ${\bf x}$ are uncorrelated.

By Definition 1, the covariance localization operator can be written as:

The choice of the filter $\gamma_{{\bf x}}$ in this result is somewhat arbitrary, but we shall soon see that we are interested in localized kernels. In that case, $\gamma_{{\bf x}}$ will be typically be the lowest degree polynomial satisfying the constraints and can be constructed using Lagrange interpolation for instance.

Definition 3 provides a fundamental property of the covariance. The size of the correlation (distance over the graph) depends on the support of localized the kernel $\mathcal{T}_{i}\gamma_{{\bf x}}$. In [10, Theorem 1 and Corollary 2], it has been proved that the concentration of $\mathcal{T}_{i}\gamma_{{\bf x}}$ around $i$ depends on the regularityA regular kernel can be well approximated by a smooth function, for instance a low order polynomial, over the spectrum of the laplacian. of $\gamma_{{\bf x}}$. For example, if $\gamma_{{\bf x}}$ is a polynomial of degree $K$, it is exactly localized in a ball of radius $K$. Hence we will be mostly interested in such low degree polynomial kernels.

The graph spectral covariance matrix of a stochastic graph signal is given by $\Gamma_{{\bf x}}=U^{*}\Sigma_{{\bf x}}U$. For a GWSS signal this matrix is diagonal and the graph power spectral density (PSD) of ${\bf x}$ becomes:

Table 1 presents the differences and the similarities between the classical and the graph case. For a regular cyclic graph (ring), the localization operator is equivalent to the traditional translation and we recover the classical cyclic-stationarity results by setting $\eta_{{\bf x}}=\mathcal{T}_{0}\gamma_{{\bf x}}$. Our framework is thus a generalization of stationarity to irregular domains.

When $\gamma_{{\bf x}}$ is a bijective function, the covariance matrix contains an important part of the graph structure: the Laplacian eigenvectorsIf the laplacian contains eigenvalues with multiplicity, then the covariance matrix contains all its eigenspaces.. On the contrary, if $\gamma_{{\bf x}}$ is not bijective, some of the graph structure is lost as it is not possible to recover all eigenvectors. This is for instance the case when the covariance matrix is low-rank. As another example, let us consider completely uncorrelated centered samples with variance $1$. In this case, the covariance matrix becomes $\Sigma_{{\bf x}}=I$ and loses all graph information, even if by definition the stochastic signal remains stationary on the graph.

One of the crucial benefits of stationarity is that it is preserved by filtering, while the PSD is simply reshaped by the filter. The same property holds on graphs.

When a graph filter $g$ is applied to a GWSS signal, the result remains GWSS, the mean becomes $m_{{g(L){\bf x}}}=m_{{\bf x}}g(0)$ and the PSD satisfies:

Theorem 2 provides a simple way to artificially produce stationary signals with a prescribed PSD by simply filtering white noise :

The resulting signal will be stationary with PSD $g^{2}$. In the sequel, we assume for simplicity that the signal is centered at , i.e: $m_{\bf x}=0$. Note that the input white noise could well be non-Gaussian.

2 Comparison with the work of B. Girault

Stationarity for graph signals has been defined in the past . The proposed definition is based on an isometric graph translation operator defined for a graph signal $s$ as:

where $b(x)=\exp\left(j2\pi\sqrt{\frac{x}{\rho_{\mathcal{G}}}}\right)$ and $\rho_{\mathcal{G}}$ is an upper bound$\rho_{\mathcal{G}}=\max_{i\in\mathcal{V}}\sqrt{2d[i](d[i]+\bar{d}[i])}$ where $\bar{d}[i]=\frac{\sum_{n}=1^{N}W[i,n]d[n]}{d[i]}$ on $\lambda_{\rm max}$. While this operator conserves the energy of the signal ($\|T_{B}s\|_{2}=\|s\|_{2}$), it does not have localization properties. In a sense, one trades localization for isometry. Using this operator, the stationarity definition of Girault is a natural extension of the classical case (Definition 2).

[2, Definition 16] A stochastic signal ${\bf x}$ on the graph $\mathcal{G}$ is Wide-Sense Stationary (WSS) if and only if

To verify the stationary property of this signal, let us observe this quantity in the spectral domain:

Another key difference is that our definition allows us to generalize the notion of PSD to the graph setting in a simpler manner. To extend the notion of PSD using Girault’s definition, one would have to deal with a block diagonal structure of the covariance matrix in the spectral domain that changes depending on the choice of eigenvectors at eigenvalue multiplicitiesFor a subspace associated with an eigenvalue with multiplicity greater than one, there exist multiple possible sets of eigenvectors..

3 Gaussian random field interpretation

The framework of stationary signals on graphs can be interpreted using Gaussian Markov Random Field (GMRF). Let us assume that the signal ${\bf x}$ is drawn from a distribution

In other words, assuming a GRF probabilistic model with inverse covariance matrix $p(L)$ leads to a stationary graph signal with a $\textrm{PSD}=p^{-1}$. However a stationary graph signal is not necessarily a GRF. Indeed, stationarity assumes statistical properties on the signal that are not necessarily based on Gaussian distribution.

In Section 3 of , Gadde and Ortega have presented a GMRF model for graph signals. But they restrict themselves to the case where $p(L)=L+\delta I$. Following a similar approach Zhang et al. link the inverse covariance matrix of a GMRF with the Laplacian. Our approach is much broader than these two contributions since we do not make any assumption on the function $p(L)$. Finally, we exploit properties of stationary signals, such as the characterization of the PSD, to explicitly solve signal processing problems in Section 5.

Estimation of the signal PSD

As the PSD is central to our method, we need a reliable and scalable way to compute it. Equation (7) suggests a direct estimation method using the Fourier transform of the covariance matrix. We could thus estimate the covariance $\Sigma_{{\bf x}}$ empirically from $N_{s}$ realizations $\{x_{n}\}_{n=1\dots,N_{s}}$ of the stochastic graph signal ${\bf x}$, as

where $\bar{m}_{{\bf x}}[i]=\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}x_{n}[i]$. Then our estimate of the PSD would read

Unfortunately, when the number of nodes is considerable, this method requires the diagonalization of the Laplacian, an operation whose complexity in the general case scales as $O(N^{3})$ for the number of operations and $O(N^{2})$ for memory requirements. Additionally, when the number of available realizations $N_{s}$ is small, it is not possible to obtain a good estimate of the covariance matrix. To overcome these issues, inspired by Bartlett and Welch , we propose to use a graph generalization of the Short Time Fourier transform to construct a scalable estimation method.

Bartlett’s method can be summarized as follows. After removing the mean, the signal is first cut into equally sized segments without overlap. Then, the Fourier transform of each segment is computed. Finally, the PSD is obtained by averaging over segments the squared amplitude of the Fourier coefficients. Welch’s method is a generalization that works with overlapping segments.

where $g$ is the window used for the STFT. This is shown in Figure 4.

Our method is based on this idea, using the windowed graph Fourier transform . Instead of a translated rectangular window in time, we use a kernel $g$ shifted by multiples of a step $\tau$ in the spectral domain, i.e.

We then localize each spectral translation at each individual node of the graph. The coefficients of the graph windowed Fourier transform can be seen as a matrix with elements

where $x$ is a single realization of the stationary stochastic graph signal ${\bf x}$. This estimator provides a discrete approximation of the PSD. Interpolation is used to obtain a continuous estimator. This approach avoids the computation of the eigenvectors and the eigenvalues of the Laplacian.

Finally, the last step consists in computing the ratio between the two quantities and interpolating the discrete points (m\tau,\big{(}g\ast\gamma_{{\bf x}}\big{)}(m\tau)).

Variance of the estimator

Studying the bias of (10) reveals its interest :

where ${\bf x}$ is the stationary stochastic graph signal. For a filter $g$ well concentrated at the origin, (11) gives a smoothed estimate of $\gamma_{{\bf x}}(m\tau)$. This smoothing corresponds to the windowing operation in the vertex domain: the less localized the kernel $g$ in the spectral domain, the more pronounced the smoothing effect in (11) and the more concentrated the window in the vertex domain. It is very interesting to note we recover the traditional trade-off between bias and variance in non-parametric spectral estimation. Indeed, if $g$ is very sharply localized on the spectrum, ultimately a Dirac delta, the estimator (10) is unbiased. Let us now study the variance. Intuitively, if the signal is correlated only over small regions of the vertex set, we could isolate them with localized windows of a small size and averaging those uncorrelated estimates together would reduce the variance. These small size windows on the vertex set correspond to large band-pass kernel $g_{m}$ and therefore large bias. However, if those correlated regions are large, and this happens when the PSD is localized in low-frequencies, we cannot hope to benefit from vertex-domain averaging since the graph is finite. Indeed the corresponding windows $g_{m}$ on the vertex set are so large that a single window spans the whole graph and there is no averaging effect: the variance increases precisely when we try to suppress the bias.

Experimental assessment of the method

Figure 5 shows the results of our PSD-estimation algorithm on a $10$-nearest neighbors graph of $20^{\prime}000$ nodes (random geometric graph, weighted with an exponential kernel) and only $K=1$ realization of the stationary graph signal. We compare the estimation using frames of $M=$ $10$, $30$, $100$ Gaussian filters. The parameters $\sigma$ and $\tau$ are adapted to the number of filters such that the shifted windows have an overlap of approximately $2$ ($\tau=\sigma^{2}=\frac{(M+1)\lambda_{\max}}{M^{2}}$). For this experiment $K_{2}$ is set to $4$ and the Chebysheff polynomial order is $30$ The estimated curves are smoothed versions of the PSD.

Complexity analysis

The approximation scales with the number of edges of the graph $\mathcal{O}(|\mathcal{E}|)$, (which is proportional to $N$ in many graphs). Precisely, our PSD estimation method necessitates $(K+K_{2})M$ filtering operations (with $M$ the number of shifts of $g$). A filtering operation costs approximately $O_{c}|E|$, with $O_{c}$ the order of the Chebysheff polynomial . The final computational cost of the method is thus $\mathcal{O}\left(O_{c}(K+K_{2})M|\mathcal{E}|\right)$.

Error analysis

The difference between the approximation and the exact PSD is caused by three different factors.

The inherent bias of the estimator, which is now directly controlled by the parameter $\sigma$.

We use a fast-filtering method based on a polynomial approximation of the filter. For a rough approximation, $\sigma\gg\frac{\lambda_{\rm max}}{N}$, this error is usually negligible. However, in the other cases, this error may become large.

Graph Wiener filters and optimization framework

To recover ${\bf x}$, Wiener filters can be extended to the graph case:

The expression above can be derived by exactly mimicking the classical case and minimizes the expected quadratic error, which can be written as:

where $\bar{{\bf x}}=g(L){\bf y}$ is the estimator of ${\bf x}$ given ${\bf y}$. Theorem 5 proves the optimality of this filter for the graph case.

Wiener optimization

In this contribution, we would like to address a more general problem. Let us suppose that our measurements are generated as:

Our solution to overcome these issues is to solve the following optimization problem that we suggestively call Wiener optimization

Notice that compared to (15), the parameter $\beta$ is exchanged with the PSD of the noise. As a result, if the noise parameters are unknown, Wiener optimization does not solve completely the issue of finding the regularization parameter. In the noise-less case, one can alternatively solve the following problem

Theoretical motivations for the optimization framework

The second and main motivation is theoretical. If we have a Gaussian Random multivariate signal with i.i.d Gaussian noise, then Problem (16) is a MAP estimator.

If ${\bf x}\sim\mathcal{N}\left(0,s^{2}(L)\right)$ and ${\bf w}_{n}\sim\mathcal{N}\left(0,\sigma^{2}I\right)$, i.e: ${\bf x}$ is GWSS and Gaussian, then problem (16) is a MAP estimator for ${\bf x}|{\bf y}$

with $\Sigma_{{\bf x}{\bf y}}=s^{2}(L)H^{*}$ and $\Sigma_{{\bf y}}=Hs^{2}(L)H^{*}$

Additionally, when $H$ is jointly diagonalizable with $L$, Problem (16) can be solved by a single filtering operation.

If the operator $H$ is diagonalizable with $L$, (i.e: $H=h(L)=Ua(\Lambda)U^{*}$), then problem (16) is optimal with respect to the weighting $w$ in the sense that its solution minimizes the mean square error:

Additionally, the solution can be computed by the application of the corresponding Wiener filter.

Advantage of the Wiener optimization framework over a Gaussian MAP estimator

Theorem 3 shows that the optimization framework is equivalent to a Gaussian MAP estimator. In practice, when the data is only close to stationary, the true MAP estimator will perform better than Wiener optimization. So one could ask why we bother defining stationarity on graphs. Firstly, assuming stationarity allows us for a more robust estimate of the covariance matrix. This is shown is in Figure 5, where only one signal is used to estimate the PSD (and thus the covariance matrix). Another example is the USPS experiment presented in the next section. We estimate the PSD by using only $20$ digits. The final result is much better than a Gaussian MAP based on the empirical covariance. Secondly, we have a scalable solution for Problem (16) (See Algorithm 1 below). On the contrary the classical Gaussian MAP estimator requires the explicit computation of a large part of the covariance matrix and it’s inverse, which are both not scalable operations.

Solving Problem (16)

Note that Problem (16) can be solved with a simple gradient descent. However, for a large number of nodes $N$, the matrix $w(L)$ requires $\mathcal{O}(N^{3})$ operations to be computed and $\mathcal{O}(N^{2})$ bits to be stored. This difficulty can be overcome by applying its corresponding filter operator at each iteration. As already mentioned, the cost of the approximation scale with the number of edges $\mathcal{O}(O_{c}|E|)$ .

Evidence of graph stationarity: illustration with USPS

Stationarity may not be an obvious hypothesis for a general dataset, since our intuition does not allow us to easily capture the kind of shift invariance that is really implied. In this section we give additional insights on stationarity from a more experimental point of view. To do so, we will show that the well-known USPS dataset is close to stationary on a nearest neighbor graph. We show similar results with a dataset of faces.

Images can be considered as signals on the 2-dimensional euclidean plane and, naturally, when the signal is sampled, a grid graph is used as a discretization of this manifold. The corresponding eigenbasis is the 2 dimensional DCTThis is a natural extension of . Many papers have exploited the fact that natural texture images are stationary 2-dimensional signals , i.e stationary signals on the grid graph. In , the authors go one step further and ask the following question: suppose that pixels of images have been permuted, can we recover their relative two-dimensional location? Amazingly, they answer positively adding that only a few thousand images are enough to approximately recover the relative location of the pixels. The grid graph seems naturally encoded within images.

The observation of motivates the following experiment involving stationarity on graphs. Let us select the USPS data set which contains $9298$ digit images of $16\times 16$ pixels. We create 5 classes of data: (a) the circularly shifted digitsWe performed all possible shifts in both directions. Because of this, the covariance matrix becomes Toeplitz, (b) the original digits and (c), (d) and (e) the classes of digit $3$, $7$ and $9$. As a pre-processing step, we remove the mean of each pixel, thus forcing the first moment to be , and focus on the second moment. For those 5 cases, we compute the covariance matrix $\Sigma$ and its "Fourier transform",

for 2 different graphs: (a) the grid and (b) the $20$ nearest neighbors graph. In this latter case, each node is a pixel and is associated to a feature vector containing the corresponding pixel value of all images. We use the squared euclidean distance between feature vectors and an exponential kernel to define edge weights$W[i,n]=e^{\frac{-\|x_{i}-x_{n}\|_{2}^{2}}{\sigma^{2}}}$ if $x_{i}$ is in the $20$ nearest neighbors of $x_{n}$.. We then compute the stationarity level of each class of data with both graphs using the following measure:

The closer $s_{r}(\Gamma)$ is to $1$, the more diagonal the matrix $\Gamma$ is and the more stationary the signal. Table 2 shows the obtained stationarity measures. The less universal the data, the less stationary it is on the grid. Clearly, specificity inside the data requires a finer structure than a grid. This is confirmed by the behavior of the nearest neighbors graph. When only one digit class is selected, the nearest neighbors graph still yields very stationary signals.

To further illustrate this phenomenon on a different dataset, we use the CMUPIE set of cropped faces. With a nearest neighbor graph we obtained a stationarity level of $s_{r}=0.92$. This has already been observed in where the concept of Laplacianfaces is introduced. Finally, in the authors successfully use the graph between features to improve the quality of a low-rank recovery problem. The reason seems to be that the principal components of the data are the lowest eigenvectors of the graph, which is again a stationarity assumption.

To intuitively motivate the effectiveness of nearest neighbors at producing stationary signals, let us define the centering operator $J=I-\boldsymbol{1}\boldsymbol{1}^{\top}/N$. Given $K$ signal $x_{k}$, the matrix of average squared distances between the centered features ($\sum_{i=1}^{N}x_{k}[i]=0$) is directly proportional to the covariance matrix :

where $D[i,n]=\frac{1}{K}\sum_{k=1}^{K}\left(x_{k}[i]-x_{k}[n]\right)^{2}$ and $\bar{\Sigma}_{{\bf x}}[i,n]=\frac{1}{K}\sum_{k=1}^{K}x_{k}[i]x_{k}[n]$. The proof is given in Appendix E. The nearest-neighbors graph can be seen as an approximation of the original distance matrix, which pleads for using it as a good proxy destined to leverage the spectral content of the covariance. Put differently, when using realizations of the signal as features and computing the k-NN graph we are connecting strongly correlated variables via strong edge weights.

Experiments

All experiments were performed with the GSPBox and the UNLocBoX two open-source software library. The code to reproduce all figures of the paper can be downloaded at: https://lts2.epfl.ch/rrp/stationarity/. As the stationary signals are random, the reader may obtain slightly different results. However, conclusions shall remain identical. The models used in our comparisons are detailed in the Appendix A for completeness, where we also detail how the tuning of the parameters is done. All experiments are evaluated with respect to the Signal to Noise Ratio (SNR) measure:

In order to obtain a first insight into applications using stationarity, we begin with some classical problems solved on a synthetic dataset. Compared to real data, this framework allows us to be sure that the signal is stationary on the graph.

We start with a de-convolution example on a random geometric graph. This can model an array of sensors distributed in space or simply a mesh. The signal is chosen with a low frequency band-limited PSD. To produce the measurements, the signal is convolved with the heat kernel $h(\lambda)=e^{-\tau\lambda}$. Additionally, we add some uncorrelated i.i.d Gaussian noise. The heat kernel is chosen because it simulates a heat diffusion process. Using de-convolution we aim at recovering the original signal before diffusion. For this experiment, we put ourselves in an ideal case and suppose that both the PSD of the input signal and the noise level are known.

Figure 8 presents the results. We observe that Wiener filtering is able to de-convolve the measurements. The second plot shows the reconstruction errors for three different methods: Tikhonov presented in problem (23), TV in (25) and Wiener filtering in (13). Wiener filtering performs clearly much better than the other methods because it has a much better prior assumption.

Graph Wiener in-painting

In our second example, we use Wiener optimization to solve an in-painting problem. This time, we suppose that the PSD of the input signal is unknown and we estimate it using $50$ signals. Figure 9 presents quantitative results for the in-painting. Again, we compare three different optimization methods: Tikhonov (22), TV (25) and Wiener (16). Additionally we compute the classical MAP estimator based on the empirical covariance matrix (see 2.23). Wiener optimization performs clearly much better than the other methods because it has a much better prior assumption. Even with $50$ measurements, the MAP estimator performs poorly compared to graph methods. The reason is that the graph contains a lot of the covariance information. Note that the PSD estimated with only one measurement is sufficient to outperform Tikhonov and TV.

2 Meteorological dataset

We apply our methods to a weather measurements dataset, more precisely to the temperature and the humidity. Since intuitively these two quantities are correlated smoothly across space, it suggests that they are more or less stationary on a nearest neighbors geographical graph.

The French national meteorological service has published in open access a datasetAccess to the raw data is possible directly through our code or through the link https://donneespubliques.meteofrance.fr/donnees_libres/Hackathon/RADOMEH.tar.gz with hourly weather observations collected during the Month of January 2014 in the region of Brest (France). From these data, we wish to ascertain that our method still performs better than the two other models (TV and Tikhonov) on real measurements. The graph is built from the coordinates of the weather stations by connecting all the neighbors in a given radius with a weight function $W[i,n]=e^{-d_{in}^{2}\tau}$ where $\tau$ is adjusted to obtain an average degree around $3$ ($\tau$, however, is not a sensitive parameter). For our experiments, we consider every time-step as an independent realization of a GWSS signal. As sole pre-processing, we remove the temperature mean of each station independently. This is equivalent to removing the first moment. Thanks to the $744$ time observation, we can estimate the covariance matrix and check whether the signal is stationary on the graph.

The result of the experiment with temperatures is displayed in Figure 10. The covariance matrix shows a strong correlation between the different weather stations. Diagonalizing it with the Fourier basis of the graph shows that the meteorological instances are not really stationary within the distance graph as the resulting matrix is not really diagonal. However, even in this case, Wiener optimization still outperforms graph TV and Tikhonov models, showing the robustness of the proposed method. In our experiment, we solve a prediction problem with a mask operator covering 50 per cent of measurements and an initial average SNR of $13.4$ dB . We then average the result over $744$ experiments (corresponding to the $744$ observations) to obtain the curves displayed in Figure 10. We observe that Wiener optimization always performs better than the two other methods.

Prediction - Humidity

Using the same graph, we have performed another set of experiments on humidity observations. The results are displayed in Figure 11. In our experiment, we solve a prediction problem with a mask operator covering $50\%$ of measurements and various amount of noise. The rest of the testing framework is identical as for the temperature and the conclusions are similar.

3 USPS dataset

We perform the same kind of in-painting/de-noising experiments with the USPS dataset. For our experiments, we consider every digit as an independent realization of a GWSS signal. As sole pre-processing, we remove the mean of each pixel separately. This ensures that the first moment is . We create the graphThe graph is created using patches of pixels of size $5\times 5$. The pixels’ patches help because we have only a few digits available. When the size of the data increases, a nearest neighbor graph performs even better. and estimate the PSD using only the first $20$ digits and we use $500$ of the remaining ones to test our algorithm. We use a mask covering $50\%$ of the pixel and various amount of noise. We then average the result over $500$ experiments (corresponding to the $500$ digits) to obtain the curves displayed in Figure 12All parameters have been tuned optimally in a probabilistic way. This is possible since the noise is added artificially. The models presented in Appendix A have only one parameter to be tuned: $\epsilon$ which is set to $\epsilon=\sigma\sqrt{\#y}$, where $\sigma$ is the variance of the noise and $\#y$ the number of elements of the vector $y$. In order to be fair with the MAP estimator, we construct the graph with the only $20$ digits used in the PSD estimation. . For this experiment, we also compare with traditional TV de-noising and Tikhonov de-noising. The optimization problems used are similar to (22). Additionally we compute the classical MAP estimator based on the empirical covariance matrix for the solution see ( 2.23). The results presented in Figure 12 show that graph optimization is outperforming classical techniques, meaning that the grid is not the optimal graph for the USPS dataset. Wiener once again outperforms the other graph-based models. Moreover, this experiment shows that our PSD estimation is robust when the number of signals is small. In other words, using the graph allows us for a much better covariance estimation than a simple empirical average. When the number of measurements increases, the MAP estimator improves in performance and eventually outperforms Wiener because the data is close to stationary on the graph.

4 ORL dataset

For this last experiment, we use the ORL face dataset. We have a good indication that this dataset is close to stationary since CMUPIE (a smaller faces dataset) is also close to stationary. Each image has $112\times 92=10304$ pixels making it complicated to estimate the covariance matrix and to use a Gaussian MAP estimator. Wiener optimization, on the other hand, does not necessitate an explicit computation of the covariance matrix. Instead, we estimate the PSD using the algorithm presented in Section 4. A detailed experiment is performed in Figure 13. After adding Gaussian noise to the image, we remove randomly a percentage of the pixels. We consider the obtained image as the measurement and we reconstruct the original image using TV, Tikhonov and Wiener priors. In Figure 14, we display the reconstruction results for various noise levels. We create the graph with $300$ facesWe build a nearest neighbor graph based on the pixels values. and estimate the PSD with $100$ faces. We test the different algorithms on the $100$ remaining faces.

Conclusion

In this contribution, we have extended the common concept of stationarity to graph signals. Using this statistical model, we proposed a new regularization framework that leverages the stationarity hypothesis by using the Power Spectral Density (PSD) of the signal. Since the PSD can be efficiently estimated, even for large graphs, the proposed Wiener regularization framework offers a compelling way to solve traditional problems such as denoising, regression or semi-supervised learning. We believe that stationarity is a natural hypothesis for many signals on graphs and showed experimentally that it is deeply connected with the popular nearest neighbor graph construction. As future work, it would be very interesting to clarify this connection and explore if stationarity could be used to infer the graph structure from training signals, in the spirit of .

Acknowledgments

We thank the anonymous reviewers for their constructive comments that helped us improve the structure of the paper, especially Section III B. We also thank Andreas Loukas, Vassilis Kalofolias and Nauman Shahid for their useful suggestions.

This work has been supported by the Swiss National Science Foundation research project Towards Signal Processing on Graphs, grant number: 2000_21/154350/1.

Appendix A Convex models

Convex optimization has recently become a standard tool for problems such as de-noising, de-convolution or in-painting. Graph priors have been used in this field for more than a decade . The general assumption is that the signal varies smoothly along the edges, which is equivalent to saying that the signal is low-frequency-based. Using this assumption, one way to express mathematically an in-painting problem is the following:

where $M$ is a masking operator and $\epsilon$ a constant computed thanks to the noise level. We could also rewrite the objective function as $x^{T}Lx+\gamma\|Mx-y\|_{2}^{2}$, but this implies a greedy search of the regularization parameter $\gamma$ even when the level of noise is known. For our simulations, we use Gaussian i.i.d. noise of standard deviation $n$. It allows us to optimally set the regularization parameter $\epsilon=n\sqrt{\#y}$, where $\#y$ is the number of elements of the measurement vector.

Graph de-convolution can also be addressed with the same prior assumption leading to

where $h$ is the convolution kernel. To be as generic as possible, we combine problems (22) and (23) together leading to a model capable of performing de-convolution, in-painting and de-noising at the same time:

In order to solve these problems, we use a subset of convex optimization tools called proximal splitting methods. Since we are not going to summarize them here, we encourage a novice reader to consult and the references therein for an introduction to the field.

Appendix B Proof of Theorem 3

The proof is a classic development used in Bayesian machine learning. By assumption ${\bf x}$ is a sample of a Gaussian random multivariate signal ${\bf x}\sim\mathcal{N}\left(m_{{\bf x}},s^{2}(L)\right)$. The measurements are given by

where ${\bf w}_{n}\sim\mathcal{N}\left(0,\sigma^{2}\right)$ and thus have the following first and second moments: ${\bf y}|{\bf x}\sim\mathcal{N}\left(H{\bf x},\sigma^{2}I)\right)$. For simplicity, we assume $s^{2}(L)$ to be invertible. However this assumption is not necessary. We can write the probabilities of ${\bf x}$ and ${\bf y}|{\bf x}$ as:

Appendix C Proof of Theorem 5

Since by hypothesis $H=h(L)=Uh(\Lambda)U^{*}$, we can rewrite the optimization problem (16) in the graph Fourier domain using the Parseval identity $\|{\bf x}\|_{2}=\|U{\bf x}\|_{2}=\|\hat{{\bf x}}\|_{2}$:

As a next step, we use the fact that $\hat{{\bf y}}=h\hat{{\bf x}}+\hat{{\bf w}}_{n}$ to find:

The error performed by the algorithm becomes

The expectation of the error can thus be computed:

which is the Wiener filter associated with the convolution $h(L)=H$. ∎

Appendix D Proof of Theorem 4

Let ${\bf x}$ be GWSS with covariance matrix $\Sigma_{{\bf x}}=s^{2}(L)$ and mean $m_{\bf x}$. The measurements satisfy

where ${\bf w}_{n}$ is i.i.d noise with PSD $\sigma^{2}$. The variable ${\bf y}$ has a covariance matrix $\Sigma_{{\bf y}}=Hs^{2}(L)H^{*}+\sigma^{2}I$ and a mean $m_{{\bf y}}=Hm_{\bf x}$. The covariance between ${\bf x}$ and ${\bf y}$ is $\Sigma_{{\bf x}{\bf y}}=\Sigma_{{\bf y}{\bf x}}^{*}=s^{2}(L)H^{*}$. For simplicity, we assume $s^{2}(L)$ and $Hs^{2}(L)H^{*}+\sigma^{2}I$ to be invertible. However this assumption is not necessary. The Wiener optimization framework reads:

where (D) follows from the Woodbury, Sherman and Morrison formula. The linear estimator of ${\bf x}$ corresponding to Wiener optimization is thus:

We observe that it is equivalent to the solution of the linear minimum mean square error estimator:

with ${\bf y}=H{\bf x}+{\bf w}_{n}$. See [44, Equation 12.6]. ∎

Using similar arguments, we can prove that

and is thus a linear minimum mean square estimator too.

Appendix E Development of equation 21

Let us denote the matrix of squared distances $D[i,j]=\tfrac{1}{K}\sum_{k}|x_{k}[i]-x_{k}[j]|^{2}$ for the samples $\{x_{1},x_{2},\ldots x_{K}\}$ of the random multivariate variable ${\bf x}$ on a $N$ vertex graph. Let us assume further that $m[k]=\sum_{n=1}^{N}x_{k}[n]=0$. We show then that $\bar{\Sigma}_{{\bf x}}=-\tfrac{1}{2}JD_{{\bf x}}J$ where $\bar{\Sigma}_{{\bf x}}$ is the covariance (Gram) matrix defined as $\bar{\Sigma}_{{\bf x}}[i,j]=\tfrac{1}{K}\sum_{k=1}^{K}x_{k}[i]x_{k}[j]$ and $J$ is centering matrix $J[k,l]=\delta_{k}[l]-\tfrac{1}{N}$.

Let us substitute $D[i,j]=\bar{\Sigma}_{{\bf x}}[i,i]+\bar{\Sigma}_{{\bf x}}[j,j]-2\bar{\Sigma}_{{\bf x}}[i,j]$, then we find

Under the assumption $m[k]=\sum_{n=1}^{N}x_{k}[n]=0$, we recover the desired result $\bar{\Sigma}_{{\bf x}}=-\tfrac{1}{2}JD_{{\bf x}}J$. ∎

References