Stationary Graph Signals

• Stationary graph signals are random signals indexed by graph vertices, where stationarity is defined relative to a graph operator rather than traditional translations.
• Their analysis relies on spectral diagonalization through a graph Fourier transform, enabling covariance to function as a localized graph filter for robust PSD estimation.
• These signals underpin advanced applications including signal filtering, graph learning, and optimal estimation techniques in network inference and prediction.

Stationary graph signals are random signals indexed by the vertices of a graph whose first- and second-order structure is defined relative to a graph operator rather than an ordinary Euclidean translation. In the foundational graph signal processing formulations, stationarity is tied to a graph Fourier basis induced by a graph Laplacian or a more general graph shift operator, so that covariance becomes a graph filter or is diagonalized by the graph Fourier transform (GFT); in constructive formulations, the process is generated by filtering white noise on the graph. Later work extended the concept to diffusion-based means, joint time-vertex processes, Hilbert-space-valued graph processes, local stationarity on irregular domains, and graph learning problems in which stationarity is used to infer the underlying network itself (Perraudin et al., 2016, Marques et al., 2016, Gama et al., 2018).

1. Foundational definitions and shift dependence

A stationary graph-signal model begins with a choice of graph operator. In the graph-Laplacian formulation, the combinatorial Laplacian is

$L=D-W,$

with eigendecomposition

$L = U\Lambda U^*,$

and the GFT is

$\hat f = U^* f.$

The localization operator

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$

plays the role of graph translation. A stochastic graph signal is graph wide-sense stationary when its mean is constant and its covariance is a localized graph kernel,

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$

This makes stationarity a statement about graph-adapted localization rather than ordinary lag invariance (Perraudin et al., 2016).

A second foundational line defines weak stationarity with respect to a normal graph shift operator

$S=V\Lambda V^H$

by requiring that the process be the output of a linear shift-invariant graph filter applied to white noise: $$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$ Under distinct eigenvalues, this is equivalent to requiring that the covariance and the graph shift be simultaneously diagonalizable, so that graph stationarity is the graph-domain analogue of the classical equivalence between filtered white noise and Fourier-diagonal covariance (Marques et al., 2016).

A third line, developed for diffusion-based graph processes, defines graph wide-sense stationarity relative to a normal, nonnegative shift $S$ by two conditions: $\mu_x = \mathbb E[x] = \mu\, v_1$ for the Perron eigenvector $v_1$, and

$L = U\Lambda U^*,$0

Here the mean is not necessarily constant across nodes; it is proportional to the dominant eigenvector of the shift. This definition is explicitly different from the constant-mean definition above and makes the graph analogue of the DC mode the Perron mode rather than the all-ones vector unless $L = U\Lambda U^*,$1 (Gama et al., 2018).

A further strengthening appears in work on symmetric shifts and graph-polynomial dynamics, where a stationary graph signal is defined to have zero mean and covariance equal to a polynomial graph filter of the shift,

$L = U\Lambda U^*,$2

This definition is stronger than mere commutation, because a covariance may commute with $L = U\Lambda U^*,$3 without being a polynomial in $L = U\Lambda U^*,$4, especially when $L = U\Lambda U^*,$5 has repeated eigenvalues (Chen et al., 16 Sep 2025).

A recurrent source of confusion is therefore definitional rather than terminological. The literature does not use a single universal mean condition. Some formulations impose a constant mean or a mean in the null space of the Laplacian; some impose zero mean; and diffusion-based formulations impose alignment with the Perron eigenvector. What remains common is that stationarity is always shift-dependent: it is defined relative to a chosen graph operator and its spectral basis, not independently of that choice (Perraudin et al., 2016, Gama et al., 2018).

2. Spectral characterization and power spectral density

Across the main formulations, the decisive second-order property is spectral diagonalization. For graph-stationary signals, the covariance takes the form

$L = U\Lambda U^*,$6

with $L = U\Lambda U^*,$7 diagonal, or equivalently

$L = U\Lambda U^*,$8

The graph power spectral density (PSD) is then the diagonal of the covariance in the graph Fourier basis,

$L = U\Lambda U^*,$9

This is the graph counterpart of the Wiener–Khintchine characterization (Perraudin et al., 2016, Marques et al., 2016).

The constructive filtering model makes the PSD interpretation immediate. If

$\hat f = U^* f.$0

then the graph-filter frequency response is

$\hat f = U^* f.$1

and for white input the output covariance is diagonalized by the same graph Fourier basis. In the weak-stationary graph-process formulation,

$\hat f = U^* f.$2

and filtering a stationary input reshapes the PSD according to

$\hat f = U^* f.$3

The graph Fourier coefficients are therefore uncorrelated, and the GFT basis is also the Karhunen–Loève basis for a stationary graph process (Marques et al., 2016).

The same spectral logic underlies Gaussian stationary-signal graph learning. If

$\hat f = U^* f.$4

then

$\hat f = U^* f.$5

for symmetric $\hat f = U^* f.$6, so $\hat f = U^* f.$7, $\hat f = U^* f.$8, and $\hat f = U^* f.$9 share eigenvectors. Provided $$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$0 is full rank, the precision $$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$1 also commutes with $$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$2,

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$3

This commutativity is the operational form used to couple Gaussian likelihoods to graph stationarity in later graph-learning methods (Buciulea et al., 2023, Buciulea et al., 2024).

PSD estimation on graphs follows both nonparametric and parametric routes. In the weak-stationary graph-process framework, the periodogram and correlogram coincide, and for $$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$4 independent realizations the periodogram is

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$5

It is unbiased, and for Gaussian signals on symmetric shifts its covariance is

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$6

with mean-squared error

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$7

The same paper develops window-based average periodograms, filter banks, and parametric MA, AR, and ARMA graph-process estimation (Marques et al., 2016).

A scalable alternative is the graph Welch/Bartlett-type estimator based on localized spectral windows. For windows $$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$8, the estimator is

$$T_i g[n] := \sum_{\ell=0}^{N-1} g(\lambda_\ell) u_\ell^*[i] u_\ell[n] = (g(L)\delta_i)[n] = g(L)[i,n]$$9

with bias

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$0

This returns a smoothed PSD and makes the bias–variance tradeoff explicit (Perraudin et al., 2016).

3. Means, realization averages, and graph ergodicity

In graph settings, the relationship between ensemble means and averages from a single realization is subtler than in ordinary time series. For diffusion-based graph stationarity, the ensemble mean is

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$1

so it is generally node-varying. The paper therefore defines a graph-domain realization average not as an arithmetic mean but as a graph shift average,

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$2

Unbiasedness requires

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$3

which yields

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$4

Operationally, each power $\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$5 diffuses information $\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$6 hops through the graph, so each node forms a weighted combination of local and multi-hop observations (Gama et al., 2018).

This estimator is itself wide-sense stationary, with output PSD

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$7

The construction is therefore a low-pass graph filter: it preserves the Perron/DC component,

$\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$8

and attenuates non-DC frequencies. At node $\Sigma_x[i,n] = T_i\gamma_x[n] = \gamma_x(L)[i,n].$9, the deviation bound is

$S=V\Lambda V^H$0

Unlike the classical variance-of-sample-mean bound, this depends explicitly on eigenvector localization at the node, so convergence is spatially nonuniform (Gama et al., 2018).

Under spectral separation conditions, the graph weak law of large numbers states that if $S=V\Lambda V^H$1 and

$S=V\Lambda V^H$2

or if $S=V\Lambda V^H$3, then

$S=V\Lambda V^H$4

This reproduces the classical order $S=V\Lambda V^H$5, but only “in at least some nodes” rather than uniformly over all nodes. The paper also gives

$S=V\Lambda V^H$6

which makes explicit that worst-node behavior can be qualitatively different (Gama et al., 2018).

A second major result is that the graph shift average is generally not mean-squared-error optimal. For unbiased graph-filter estimators

$S=V\Lambda V^H$7

with PSD

$S=V\Lambda V^H$8

the MSE criterion

$S=V\Lambda V^H$9

is minimized by the ideal low-pass graph filter

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$0

The same solution minimizes a D-optimality criterion. In words, the optimal estimator keeps only the Perron/DC mode and suppresses every other graph frequency exactly (Gama et al., 2018).

4. Joint, generalized, and locally stationary graph processes

Stationary graph-signal theory extends in several orthogonal directions. For time-varying graph signals, joint wide-sense stationarity (JWSS) is defined on the joint Fourier basis

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$1

with joint Fourier transform

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$2

A process is JWSS when

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$3

and

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$4

with $$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$5 diagonal. The diagonal entries

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$6

define the joint power spectral density (JPSD). This is the time-vertex analogue of both temporal and graph stationarity and allows non-separable graph-frequency-dependent temporal dynamics (Perraudin et al., 2016, Loukas et al., 2016).

For prediction, the joint causal model

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$7

decouples in the graph Fourier domain into $$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$8 scalar ARMA models,

$$x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l, \qquad \mathbb E\{ww^H\}=I.$$9

This decoupling theorem is the main computational consequence of joint stationarity for forecasting. Under invertibility, the one-step prediction error equals the innovation,

$S$0

so the predictor is MSE-optimal for the fitted causal model (Loukas et al., 2016).

A separate operator-theoretic development defines time-vertex stationarity by invariance under a bivariate joint translation operator

$S$1

In that formulation, a zero-mean time-vertex process is JWSS if

$S$2

for all $S$3, and this is equivalent to covariance diagonalization by the joint Fourier transform. The associated JPSD is genuinely bivariate in temporal and graph frequencies rather than only a function of a product-graph eigenvalue (Jalili et al., 2020).

Generalization in another direction replaces scalar values at each vertex by elements of a separable Hilbert space $S$4. In that setting, a generalized graph random process $S$5 is jointly wide-sense stationary when

$S$6

which implies the covariance expansion

$S$7

The coefficients

$S$8

are then uncorrelated, with variances $S$9. Standard scalar graph stationarity is recovered by taking $\mu_x = \mathbb E[x] = \mu\, v_1$0; time-vertex stationarity is recovered by taking $\mu_x = \mathbb E[x] = \mu\, v_1$1 or a corresponding temporal operator (Jian et al., 2021).

Local stationarity weakens the global model. A locally stationary graph process is defined as

$\mu_x = \mathbb E[x] = \mu\, v_1$2

with smooth membership vectors $\mu_x = \mathbb E[x] = \mu\, v_1$3 satisfying

$\mu_x = \mathbb E[x] = \mu\, v_1$4

Its vertex-frequency spectrum is

$\mu_x = \mathbb E[x] = \mu\, v_1$5

For a globally stationary process, this matrix degenerates to a rank-1 form with identical rows. Under localization and spectral-separation assumptions, the full process can be approximated locally by WSS processes on subgraphs (Canbolat et al., 2023).

These extensions suggest a common theme: graph stationarity is best viewed as a spectral symmetry condition that survives substantial changes in domain, but the symmetry may live on a graph, on a time-vertex product, on a graph–Hilbert-space product, or only locally over graph regions (Perraudin et al., 2016, Jian et al., 2021, Canbolat et al., 2023).

5. Graph learning from stationary signals

Stationarity is also a structural prior for inferring unknown graphs from data. In diffusion-based graph inference, if observations are generated by

$\mu_x = \mathbb E[x] = \mu\, v_1$6

from i.i.d. latent signals $\mu_x = \mathbb E[x] = \mu\, v_1$7, then the covariance shares eigenvectors with the unknown diffusion operator $\mu_x = \mathbb E[x] = \mu\, v_1$8. Fixing the covariance eigenvectors reduces inference to selecting admissible eigenvalues

$\mu_x = \mathbb E[x] = \mu\, v_1$9

subject to linear entrywise nonnegativity constraints, spectral bounds $v_1$0, and normalization $v_1$1. The admissible set is therefore a convex polytope in eigenvalue space (Pasdeloup et al., 2016).

This spectral-template viewpoint reappears in robust graph learning from stationary signals. The classical robust spectral-template model

$v_1$2

with a hard normalization can be infeasible. A log-barrier reformulation replaces the normalization by

$v_1$3

yielding

$v_1$4

This model is always feasible, and the associated finite-sample analysis gives non-asymptotic bounds on objective values, solution sets, and degree vectors (Liu et al., 2023).

When Gaussianity is added, graph learning can be written as a joint estimation of precision and graph operator: $v_1$5 This formulation, called GGSR, treats Graphical Lasso as the restrictive special case in which graph and precision are essentially identified, while pure stationarity-based methods correspond to covariance-eigenspace matching without Gaussian likelihood. The paper develops an alternating convex BSUM scheme and reports that the Gaussian-plus-stationary model can require about ten times fewer samples than a stationarity-only approach in a polynomial covariance setting (Buciulea et al., 2023).

A related development, Polynomial Graphical Lasso, keeps the Gaussian likelihood

$v_1$6

but allows the precision to be any polynomial of the sought graph by imposing approximate commutativity,

$v_1$7

together with sparse-graph penalties on $v_1$8. This yields a nonconvex but biconvex problem solved by alternating graph and precision updates. The model explicitly generalizes Graphical Lasso from “sparse precision = graph” to “sparse graph, precision is a graph polynomial” (Buciulea et al., 2024).

Online graph learning uses the same commutativity principle. With streaming stationary graph signals, one may update the empirical covariance recursively and solve

$v_1$9

by a proximal-gradient step at each time. Under a strong-convexity condition expressed through

$L = U\Lambda U^*,$00

the online iterates track the time-varying batch optimum within a neighborhood whose size depends on the variability of the optimum itself (Shafipour et al., 2020).

Hidden nodes modify the stationarity relation rather than invalidating it. If the full covariance and graph shift satisfy

$L = U\Lambda U^*,$01

then the observed block obeys

$L = U\Lambda U^*,$02

This observation supports both joint inference of multiple graphs with hidden variables and online graph estimation from incomplete graph signals, using convex objectives that combine sparse graph penalties with $L = U\Lambda U^*,$03-regularization on the hidden-node correction matrix $L = U\Lambda U^*,$04 (Rey et al., 2021, Buciulea et al., 2024).

6. Inference tasks, filtering, and applications

Once a stationary model is specified, it supports classical statistical operations in graph form. In the graph-only case, a noisy linear inverse problem

$L = U\Lambda U^*,$05

with stationary signal PSD $L = U\Lambda U^*,$06 and noise PSD $L = U\Lambda U^*,$07 leads to the graph Wiener filter

$L = U\Lambda U^*,$08

when $L = U\Lambda U^*,$09. More generally, the stationary prior yields the optimization

$L = U\Lambda U^*,$10

which has both MAP and LMMSE interpretations under the stated Gaussian assumptions (Perraudin et al., 2016).

In the generalized Hilbert-space setting, stationarity again yields Wiener filters acting diagonally in the joint spectral basis. For denoising with $L = U\Lambda U^*,$11, the modewise gain is

$L = U\Lambda U^*,$12

and for signal completion the estimator has an explicit form involving the projected covariance operator. The framework supports multichannel, time-vertex, and continuous-time recovery tasks (Jian et al., 2021).

For dynamical systems with polynomial state and observation matrices,

$L = U\Lambda U^*,$13

with $L = U\Lambda U^*,$14 and $L = U\Lambda U^*,$15, Kalman filtering preserves stationarity when the initial state and initial estimation error are stationary. The gain and error covariance remain polynomial graph filters,

$L = U\Lambda U^*,$16

and modewise graph-spectral recursions follow from simultaneous diagonalization. In the simulations reported in that work, Kalman filtering yields lower reconstruction error than static inverse filtering and the zero-signal baseline (Chen et al., 16 Sep 2025).

Prediction, segmentation, and mean estimation also benefit directly from stationarity. Jointly stationary graph processes admit causal predictors that outperform disjoint per-node temporal models when the JPSD is non-separable (Loukas et al., 2016). Streams of graph signals with piecewise-constant mean can be segmented in the graph spectral domain because stationarity diagonalizes the residual covariance; the resulting change-point method uses sparse GFT representations of segment means and comes with a non-asymptotic oracle inequality (Concha et al., 2020). Single-realization mean estimation can be performed by graph shift averaging or optimal graph filtering, including Gaussian-Markov random field examples in sensor networks (Gama et al., 2018).

Applications in the cited literature are correspondingly broad. They include denoising and regression on arbitrary graphs (Perraudin et al., 2016), spectral estimation on synthetic and real-world graphs (Marques et al., 2016), short-term prediction of multivariate graph processes (Loukas et al., 2016), denoising and signal completion on epilepsy, air-quality, weather, and continuous-time synthetic data (Jian et al., 2021), graph learning from synthetic networks and financial or sectoral data (Buciulea et al., 2023, Buciulea et al., 2024), robustness studies on Protein and Reddit graphs (Liu et al., 2023), topology tracking from streaming graph signals (Shafipour et al., 2020), and mean or state estimation in graph-based sensor and dynamical systems (Gama et al., 2018, Chen et al., 16 Sep 2025).

A plausible implication is that stationary graph signals now function less as a single model class than as a family of graph-spectral statistical assumptions. Their common core is covariance structure aligned with a graph operator; their main differences concern what counts as the appropriate mean, whether covariance merely commutes with the shift or must be a polynomial of it, whether the stationarity is global or local, and whether inference proceeds from covariance alone or jointly with Gaussian likelihoods and latent-variable corrections.