---
title: Stationary Graph Signals
url: https://www.emergentmind.com/topics/stationary-graph-signals
type: topic
---

# Stationary Graph Signals

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 [1601.02522] [1603.04667] [1803.04550].

## 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 [1601.02522].

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 [1603.04667].

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
\[
C_x = V \operatorname{diag}(p)\,V^H.
\]
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 \(v_1\propto \mathbf 1\) [1803.04550].

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,
\[
\operatorname{cov}(x)=h(S).
\]
This definition is stronger than mere commutation, because a covariance may commute with \(S\) without being a polynomial in \(S\), especially when \(S\) has repeated eigenvalues [2509.12605].

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 [1601.02522] [1803.04550].

## 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
\[
\Sigma_x = U\Gamma_x U^*,
\]
with \(\Gamma_x\) diagonal, or equivalently
\[
C_x = V\operatorname{diag}(p)V^H.
\]
The graph power spectral density (PSD) is then the diagonal of the covariance in the graph Fourier basis,
\[
\gamma_x(\lambda_\ell)=\left(U^*\Sigma_x U\right)_{\ell,\ell},
\qquad
p:=\operatorname{diag}(V^H C_x V).
\]
This is the graph counterpart of the Wiener–Khintchine characterization [1601.02522] [1603.04667].

The constructive filtering model makes the PSD interpretation immediate. If
\[
x = Hw, \qquad H=\sum_{l=0}^{L-1} h_l S^l,
\]
then the graph-filter frequency response is
\[
\tilde h = \Psi h,
\]
and for white input the output covariance is diagonalized by the same graph Fourier basis. In the weak-stationary graph-process formulation,
\[
p = |\tilde h|^2,
\]
and filtering a stationary input reshapes the PSD according to
\[
p_y = |\tilde h|^2 \circ p_x.
\]
The graph Fourier coefficients are therefore uncorrelated, and the GFT basis is also the Karhunen–Loève basis for a stationary graph process [1603.04667].

The same spectral logic underlies Gaussian stationary-signal graph learning. If
\[
x = Hw, \qquad \mathbb E[ww^\top]=I,
\]
then
\[
C_x = HH^\top = H^2
\]
for symmetric \(S\), so \(C_x\), \(H\), and \(S\) share eigenvectors. Provided \(C_x\) is full rank, the precision \(\Theta=C_x^{-1}\) also commutes with \(S\),
\[
S\Theta=\Theta S.
\]
This commutativity is the operational form used to couple Gaussian likelihoods to graph stationarity in later graph-learning methods [2303.07041] [2404.02621].

PSD estimation on graphs follows both nonparametric and parametric routes. In the weak-stationary graph-process framework, the periodogram and correlogram coincide, and for \(R\) independent realizations the periodogram is
\[
\hat p_{\mathrm{pg}}=\frac{1}{R}\sum_{r=1}^R |\tilde x_r|^2.
\]
It is unbiased, and for Gaussian signals on symmetric shifts its covariance is
\[
\Sigma_{\mathrm{pg}}=(2/R)\operatorname{diag}^2(p),
\]
with mean-squared error
\[
\mathrm{MSE}(\hat p_{\mathrm{pg}})=(2/R)\|p\|_2^2.
\]
The same paper develops window-based average periodograms, filter banks, and parametric MA, AR, and ARMA graph-process estimation [1603.04667].

A scalable alternative is the graph Welch/Bartlett-type estimator based on localized spectral windows. For windows \(g_m\), the estimator is
\[
\bar{\gamma}_x(m\tau)=\frac{\|g_m(L)x\|_2^2}{\|g_m(L)\|_F^2},
\]
with bias
\[
\frac{E\{\|g_m(L)\tilde{x}\|_2^2\}}{\|g_m(L)\|_F^2}
=
\frac{\sum_{\ell=0}^{N-1} g_m(\lambda_\ell)^2\,\gamma_x(\lambda_\ell)}
{\sum_{\ell=0}^{N-1} g_m(\lambda_\ell)^2}.
\]
This returns a smoothed PSD and makes the bias–variance tradeoff explicit [1601.02522].

## 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
\[
\mu_x=\mu v_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,
\[
\hat{\mu}_N := \frac{1}{\alpha(S)} \sum_{\ell=0}^{N-1} S^\ell x.
\]
Unbiasedness requires
\[
\alpha(S)=\sum_{\ell=0}^{N-1}\lambda_1^\ell,
\]
which yields
\[
\hat{\mu}_N =
\frac{1}{\sum_{\ell=0}^{N-1}\lambda_1^\ell}\sum_{\ell=0}^{N-1}S^\ell x,
\qquad
\mathbb E[\hat{\mu}_N]=\mu_x.
\]
Operationally, each power \(S^\ell x\) diffuses information \(\ell\) hops through the graph, so each node forms a weighted combination of local and multi-hop observations [1803.04550].

This estimator is itself wide-sense stationary, with output PSD
\[
q_n = p_n \left|\sum_{\ell=0}^{N-1}\lambda_n^\ell\right|^2
\Big/
\left|\sum_{\ell=0}^{N-1}\lambda_1^\ell\right|^2.
\]
The construction is therefore a low-pass graph filter: it preserves the Perron/DC component,
\[
q_1=p_1,
\]
and attenuates non-DC frequencies. At node \(k\), the deviation bound is
\[
\Pr\!\left(\left|[\hat{\mu}_N-\mu_x]_k\right|>\epsilon\right)
\le
\frac{1}{\epsilon^2}\sum_{n=1}^N q_n |v_{k,n}|^2.
\]
Unlike the classical variance-of-sample-mean bound, this depends explicitly on eigenvector localization at the node, so convergence is spatially nonuniform [1803.04550].

Under spectral separation conditions, the graph weak law of large numbers states that if \(\lambda_1>1\) and
\[
\frac{|\lambda_n|}{\lambda_1}=o\!\left(N^{-\delta/2(N-1)}\right)
\quad \text{for some } \delta>0 \text{ and all } n\ge 2,
\]
or if \(\lambda_1=1\), then
\[
\min_{k=1,\ldots,N}\Pr\!\left(\left|[\hat{\mu}_N-\mu_x]_k\right|>\epsilon\right)
\le
\frac{p_1}{N\epsilon^2}+o(N^{-\delta}).
\]
This reproduces the classical order \(O(1/(N\epsilon^2))\), but only “in at least some nodes” rather than uniformly over all nodes. The paper also gives
\[
\max_{k=1,\ldots,N}\Pr\!\left(\left|[\hat{\mu}_N-\mu_x]_k\right|>\epsilon\right)
\le
\frac{p_1}{\epsilon^2}+o(N^{-\delta}),
\]
which makes explicit that worst-node behavior can be qualitatively different [1803.04550].

A second major result is that the graph shift average is generally not mean-squared-error optimal. For unbiased graph-filter estimators
\[
z_N =
\frac{1}{\sum_{\ell=0}^{N-1} h_\ell \lambda_1^\ell}
\sum_{\ell=0}^{N-1} h_\ell S^\ell x,
\qquad
\mathbb E[z_N]=\mu_x,
\]
with PSD
\[
r_n = p_n \frac{|\tilde h_n|^2}{|\tilde h_1|^2},
\]
the MSE criterion
\[
\operatorname{tr}[C_z] =
\sum_{n=1}^N p_n \frac{|\tilde h_n|^2}{|\tilde h_1|^2}
\]
is minimized by the ideal low-pass graph filter
\[
\tilde h_n=0,\quad n=2,\ldots,N,
\qquad
\tilde h_1\neq 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 [1803.04550].

## 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
\[
U_J = U_T \otimes U_G,
\]
with joint Fourier transform
\[
\operatorname{JFT}\{x\}=U_J^*x.
\]
A process is JWSS when
\[
\mathbb E[x]=0_{NT}
\]
and
\[
\Sigma_x = U_J H U_J^*,
\]
with \(H\) diagonal. The diagonal entries
\[
H_{k,k}=h(\lambda_n,\omega_\tau)
\]
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 [1606.06962] [1607.03313].

For prediction, the joint causal model
\[
\sum_{p=0}^{P} a_p(L_G)\, x_{t-p} = \sum_{q=0}^{Q} b_q(L_G)\, \varepsilon_{t-q}
\]
decouples in the graph Fourier domain into \(N\) scalar ARMA models,
\[
\sum_{p=0}^{P} a_p(n)\,\widehat x_{t-p}(n)
=
\sum_{q=0}^{Q} b_q(n)\,\widehat \varepsilon_{t-q}(n).
\]
This decoupling theorem is the main computational consequence of joint stationarity for forecasting. Under invertibility, the one-step prediction error equals the innovation,
\[
e_t=\varepsilon_t,
\]
so the predictor is MSE-optimal for the fitted causal model [1607.03313].

A separate operator-theoretic development defines time-vertex stationarity by invariance under a bivariate joint translation operator
\[
T_J^{(\upsilon,\vartheta)} = T_D^\upsilon \otimes T_G^\vartheta.
\]
In that formulation, a zero-mean time-vertex process is JWSS if
\[
R_x
=
\mathbb E\!\left[
\big(T_J^{(\upsilon,\vartheta)}x\big)
\big(T_J^{(\upsilon,\vartheta)}x\big)^*
\right]
\]
for all \((\upsilon,\vartheta)\), 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 [2004.00298].

Generalization in another direction replaces scalar values at each vertex by elements of a separable Hilbert space \(H\). In that setting, a generalized graph random process \(X\) is jointly wide-sense stationary when
\[
C_X\circ (A_G\otimes A_H)=(A_G\otimes A_H)\circ C_X,
\]
which implies the covariance expansion
\[
C_X = \sum_{k,\tau} p_X(k,\tau)\, P_{\phi_k\otimes \psi_\tau}.
\]
The coefficients
\[
\widehat X_{k,\tau}=\langle X,\phi_k\otimes \psi_\tau\rangle
\]
are then uncorrelated, with variances \(p_X(k,\tau)\). Standard scalar graph stationarity is recovered by taking \(H=\mathbb C\); time-vertex stationarity is recovered by taking \(H=\mathbb C^T\) or a corresponding temporal operator [2112.01127].

Local stationarity weakens the global model. A locally stationary graph process is defined as
\[
x=\sum_{k=1}^K M_k\,U\,h_k(\Lambda)\,U^T w
\]
with smooth membership vectors \(m_k\) satisfying
\[
m_k^T L m_k \le C
\quad\text{and}\quad
\|h_k(\Lambda)\|_F^2=1.
\]
Its vertex-frequency spectrum is
\[
M=\sum_{k=1}^K m_k h_k^T.
\]
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 [2309.01657].

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 [1606.06962] [2112.01127] [2309.01657].

## 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
\[
x_i = T^{k(i)} y_i
\]
from i.i.d. latent signals \(y_i\), then the covariance shares eigenvectors with the unknown diffusion operator \(T\). Fixing the covariance eigenvectors reduces inference to selecting admissible eigenvalues
\[
\widetilde T = \mathcal X \operatorname{diag}(\widetilde\lambda)\mathcal X^\top
\]
subject to linear entrywise nonnegativity constraints, spectral bounds \(\widetilde\lambda_i\in[-1,1]\), and normalization \(\widetilde\lambda_1=1\). The admissible set is therefore a convex polytope in eigenvalue space [1605.02569].

This spectral-template viewpoint reappears in robust graph learning from stationary signals. The classical robust spectral-template model
\[
\|\!C_n S - S C_n\!\|_F \le \delta
\]
with a hard normalization can be infeasible. A log-barrier reformulation replaces the normalization by
\[
-\alpha \mathbf 1^\top \log(S\mathbf 1),
\]
yielding
\[
\min_S \|S\|_{1,1} - \alpha \mathbf 1^\top \log(S\mathbf 1)
\quad
\text{s.t.}\quad
S\in\mathcal S,\ 
\|SC_n-C_nS\|_F^2\le \delta_n^2.
\]
This model is always feasible, and the associated finite-sample analysis gives non-asymptotic bounds on objective values, solution sets, and degree vectors [2305.01379].

When Gaussianity is added, graph learning can be written as a joint estimation of precision and graph operator:
\[
\hat \Theta,\hat S
=
\arg\min_{\Theta \succeq 0,\; S\in\mathcal S}
-\log\det(\Theta)+\operatorname{tr}(\Theta \hat C)+\rho\|S\|_0
\quad
\text{s.t.}\quad
S\Theta=\Theta S.
\]
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 [2303.07041].

A related development, Polynomial Graphical Lasso, keeps the Gaussian likelihood
\[
-\log(\det(\Theta))+\operatorname{tr}(\hat C\Theta)
\]
but allows the precision to be any polynomial of the sought graph by imposing approximate commutativity,
\[
\|\Theta S-S\Theta\|_F\le \delta,
\]
together with sparse-graph penalties on \(S\). 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” [2404.02621].

Online graph learning uses the same commutativity principle. With streaming stationary graph signals, one may update the empirical covariance recursively and solve
\[
S^\star \in \arg\min_{S\in\mathcal S}
\|S\|_1 + \frac{\mu}{2}\|\hat C_t S-S\hat C_t\|_F^2
\]
by a proximal-gradient step at each time. Under a strong-convexity condition expressed through
\[
Q_t = \hat C_t\otimes I_N - I_N\otimes \hat C_t,
\]
the online iterates track the time-varying batch optimum within a neighborhood whose size depends on the variability of the optimum itself [2007.03653].

Hidden nodes modify the stationarity relation rather than invalidating it. If the full covariance and graph shift satisfy
\[
CS=SC,
\]
then the observed block obeys
\[
C_{\mathcal O} S_{\mathcal O} + P = S_{\mathcal O} C_{\mathcal O} + P^\top,
\qquad
P := C_{\mathcal O\mathcal H} S_{\mathcal H\mathcal O}.
\]
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 \(\ell_{2,1}\)-regularization on the hidden-node correction matrix \(P\) [2110.03666] [2409.08760].

## 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
\[
y = Hx + w_n
\]
with stationary signal PSD \(s^2(\lambda)\) and noise PSD \(n(\lambda)\) leads to the graph Wiener filter
\[
g(\lambda)=\frac{h(\lambda)s^2(\lambda)}
{h^2(\lambda)s^2(\lambda)+n(\lambda)}
\]
when \(H=h(L)\). More generally, the stationary prior yields the optimization
\[
\bar{x}|y = \argmin_x \|Hx-y\|_2^2 + \|w(L)(x-m_x)\|_2^2,
\qquad
w(\lambda)=\frac{\sqrt{n(\lambda)}}{s(\lambda)},
\]
which has both MAP and LMMSE interpretations under the stated Gaussian assumptions [1601.02522].

In the generalized Hilbert-space setting, stationarity again yields Wiener filters acting diagonally in the joint spectral basis. For denoising with \(Y=X+E\), the modewise gain is
\[
g_{k,\tau}=
\frac{p_X(k,\tau)}{p_X(k,\tau)+p_E(k,\tau)},
\]
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 [2112.01127].

For dynamical systems with polynomial state and observation matrices,
\[
x_k=A_k x_{k-1}+\sigma_k e_k,\qquad
z_k=B_k x_k+\tilde\sigma_k \tilde e_k,
\]
with \(A_k=a_k(S)\) and \(B_k=b_k(S)\), Kalman filtering preserves stationarity when the initial state and initial estimation error are stationary. The gain and error covariance remain polynomial graph filters,
\[
K_k=g_k(S),\qquad P_k=p_k(S),
\]
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 [2509.12605].

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 [1607.03313]. 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 [2006.10628]. Single-realization mean estimation can be performed by graph shift averaging or optimal graph filtering, including Gaussian-Markov random field examples in sensor networks [1803.04550].

Applications in the cited literature are correspondingly broad. They include denoising and regression on arbitrary graphs [1601.02522], spectral estimation on synthetic and real-world graphs [1603.04667], short-term prediction of multivariate graph processes [1607.03313], denoising and signal completion on epilepsy, air-quality, weather, and continuous-time synthetic data [2112.01127], graph learning from synthetic networks and financial or sectoral data [2303.07041] [2404.02621], robustness studies on Protein and Reddit graphs [2305.01379], topology tracking from streaming graph signals [2007.03653], and mean or state estimation in graph-based sensor and dynamical systems [1803.04550] [2509.12605].

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.

Source: https://www.emergentmind.com/topics/stationary-graph-signals