Time-Vertex Graph Signal Analysis

• Time-vertex graph signals are defined on graph vertices evolving over time, coupling spatial structure with temporal dynamics for localized and effective signal processing.
• Joint spectral frameworks, including the time-vertex Fourier transform and its fractional generalizations, enable precise denoising, reconstruction, and adaptive filtering on irregular domains.
• Sampling and reconstruction theories for TVGS emphasize minimal critical sampling and structured recovery methods, with practical applications in traffic analysis, climate data, and network dynamics.

Time-vertex graph signals (TVGS) are signals defined on graph vertices and evolving over time, typically represented either as functions $f \in L^2(V \times T)$ or as matrices $X \in \mathbb{C}^{N \times T}$ whose rows are temporal signals on graph vertices and whose columns are graph signals at individual time instants. The framework couples graph structure with temporal dynamics through joint transforms, joint Laplacians, stochastic models, and sampling operators, making it suitable for localization, denoising, interpolation, propagation analysis, and reconstruction on irregular spatiotemporal domains (Perraudin et al., 2016).

1. Formal signal model and graph-time representation

A standard TVGS model treats the signal as $X \in \mathbb{C}^{N \times T}$, where each row is a temporal signal on a graph vertex and each column is a graph signal at a time instant. This representation appears in finite, infinite-length, and continuous-time settings. In the finite case, the signal is indexed on a graph and a finite set of time instants; in the continuous-time case, the temporal transform is an FT; in the infinite discrete-time case, it is a DTFT; and in the finite discrete-time case, it is a DFT (Sheng et al., 29 Aug 2025).

A closely related structural viewpoint comes from finite discrete time-varying graphs (TVGs), represented as

$$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$

Here, a dynamic edge is a quadruple $e=(u,t_a,v,t_b)$, and a temporal node is an ordered pair $(u,t_a)\in V\times T$. The model distinguishes spatial edges, temporal edges, mixed edges, and spatial-temporal self-loops; when time is totally ordered, edges may also be progressive or regressive. Regressive edges directly encode cyclic or periodic behavior, such as $(u,t_{p-1},v,t_0)$ when $t_{p-1}>t_0$ (Wehmuth et al., 2014).

This distinction is important because part of the TVGS literature assumes a fixed spatial graph together with an explicit temporal structure, while TVG models allow the topology itself to vary over time. In the TVG formalism, the set of temporal nodes $VT(H)=V\times T$ induces an isomorphism with a static directed graph over $VT(H)$, enabling adjacency-tensor, adjacency-matrix, and incidence-matrix representations. For practical cases, the asymptotic memory complexity is $X \in \mathbb{C}^{N \times T}$0 (Wehmuth et al., 2014).

In product-graph formulations used extensively in signal processing, time is often represented by a ring or cycle graph, and the joint Laplacian acts on a matrix signal as

$X \in \mathbb{C}^{N \times T}$1

with the temporal Laplacian frequently realized as a second-order difference operator. This gives a unified linear operator on the joint time-vertex domain and underlies many transform, filtering, and stationarity constructions (Grassi et al., 2016).

2. Joint spectral transforms and fractional generalizations

The central spectral operator for TVGS is the joint time-vertex Fourier transform (JFT). In matrix form, one common definition is

$X \in \mathbb{C}^{N \times T}$2

or equivalently, in product-graph notation,

$X \in \mathbb{C}^{N \times T}$3

In continuous notation, the JFT can be written as

$X \in \mathbb{C}^{N \times T}$4

These representations map the signal from $X \in \mathbb{C}^{N \times T}$5 into a joint spectral-frequency domain (Zhao et al., 2024).

The joint time-vertex fractional Fourier transform (JFRT) generalizes JFT by introducing two fractional orders: $X \in \mathbb{C}^{N \times T}$6 When $X \in \mathbb{C}^{N \times T}$7, it reduces to JFT; when $X \in \mathbb{C}^{N \times T}$8, it reduces to the identity. The JFRT inherits index additivity, reversibility, reduction to identity, reduction to ordinary JFT, and, for unitary graph Fourier transforms, unitarity. It is also separable and commutative, and if both the time domain and the graph are circular, it reduces to the classical 2D-Discrete Fractional Fourier Transform (Alikaşifoğlu et al., 2022).

Fractional analysis is used because ordinary time and graph spectral domains may be suboptimal for some data. In denoising experiments, minimal MSE was achieved at fractional orders not equal to $X \in \mathbb{C}^{N \times T}$9 in both domains: $X \in \mathbb{C}^{N \times T}$0 for Molene, $X \in \mathbb{C}^{N \times T}$1 for yearly NOAA, and $X \in \mathbb{C}^{N \times T}$2 for monthly NOAA, indicating that the best-performing domain did not coincide with ordinary JFT (Alikaşifoğlu et al., 2022).

A later extension is the dynamic multiple-parameter joint time-vertex fractional Fourier transform (DMPJFRFT), which assigns distinct fractional orders to each time step. In the Type-I-I case,

$X \in \mathbb{C}^{N \times T}$3

with $X \in \mathbb{C}^{N \times T}$4. This dynamic parameterization allows a separate set of fractional graph orders for each time instance, so the graph spectral decomposition adapts from slice to slice. The transform reduces to identity when all orders are zero, to JFRT when the orders are uniform, and to JFT when all orders are $X \in \mathbb{C}^{N \times T}$5. The reported computational cost per transformation is $X \in \mathbb{C}^{N \times T}$6 (Cui et al., 20 Nov 2025).

3. Joint stationarity and stochastic process models

A major statistical development is joint wide-sense stationarity (JWSS), also called joint stationarity. In this framework, a time-vertex process has constant mean and covariance given by a joint filter,

$X \in \mathbb{C}^{N \times T}$7

where $X \in \mathbb{C}^{N \times T}$8 is a non-negative real-valued joint power spectral density (JPSD). Equivalently, the covariance is jointly diagonalizable in the time-vertex Fourier basis. This replaces an $X \in \mathbb{C}^{N \times T}$9 covariance estimation problem with estimation of a structured two-dimensional spectrum (Loukas et al., 2016).

JWSS is strictly more general than vertex stationarity on a product graph and than separable products of time- and vertex-stationary models. Product-graph based stationarity restricts the JPSD to be a function only of $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$0, whereas JWSS allows an arbitrary function $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$1. This is a substantive distinction: the model admits arbitrary interactions in the joint frequency plane rather than only separable or diagonal structure (Loukas et al., 2016).

The practical consequences are explicit. For jointly stationary processes, one can reliably learn the covariance structure from as little as a single realization of the process, and solve MMSE recovery problems, such as interpolation and denoising, in computational time nearly linear on the number of edges and timesteps (Loukas et al., 2016). A sample JPSD estimator is

$$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$2

and a smoothed convolutional estimator introduces a bias-variance tradeoff through a bivariate window (Loukas et al., 2016).

The joint Wiener framework formulates inverse problems as

$$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$3

with $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$4. Under masking or additive noise and joint stationarity, this recovers the joint Wiener filter and yields scalable denoising, recovery, and semi-supervised learning (Perraudin et al., 2016).

Parametric stochastic modeling extends this line through graph ARMA processes. In one formulation,

$$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$5

with joint spectral transfer function

$$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$6

The JS-ARMA procedure first estimates a rough JPSD from incomplete realizations and then refines it by projecting onto the spectrum manifold of graph ARMA processes through convex relaxations. The reported estimation error in the JPSD and the missing-value recovery error decays as $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$7, where $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$8 is the number of realizations (Guneyi et al., 2023).

4. Localization, uncertainty principles, and concentrated dictionaries

Localization theory in TVGS formalizes how signal energy is distributed across the vertex-time and spectral-frequency domains. For a generalized graph signal $$H = (V, E, T), \qquad E \subseteq V \times T \times V \times T .$$9, the joint vertex-time spread over $e=(u,t_a,v,t_b)$0 is

$e=(u,t_a,v,t_b)$1

and the joint spectral-frequency spread over $e=(u,t_a,v,t_b)$2 is

$e=(u,t_a,v,t_b)$3

These are proportions of signal energy localized in chosen regions of the joint domain and joint spectrum (Zhao et al., 2024).

The corresponding uncertainty principle is expressed by

$e=(u,t_a,v,t_b)$4

together with three additional inequalities involving complements of $e=(u,t_a,v,t_b)$5 and $e=(u,t_a,v,t_b)$6. The feasible region for $e=(u,t_a,v,t_b)$7 contains allowed and forbidden zones, so a signal cannot be simultaneously highly concentrated in both domains unless the supports satisfy specific spectral-structural relationships (Zhao et al., 2024).

Perfect simultaneous localization occurs if and only if

$e=(u,t_a,v,t_b)$8

that is, $e=(u,t_a,v,t_b)$9 is an eigenvector of $(u,t_a)\in V\times T$0 with eigenvalue $(u,t_a)\in V\times T$1. In general this condition does not hold, so practical constructions focus on signals that are maximally concentrated in both domains. These are obtained as the leading eigenvectors of $(u,t_a)\in V\times T$2 under bandlimiting and orthogonality constraints, and they serve as jointly energy-concentrated atoms (Zhao et al., 2024).

The resulting joint vertex-time dictionary is built by combining atoms localized over multiple vertex and temporal subsets: $(u,t_a)\in V\times T$3 or, in compact notation,

$(u,t_a)\in V\times T$4

For real data, localization centers and widths are optimized by alternating minimization over sparse representation coefficients and localization parameters, with sparsity promoted by LASSO regularization. On traffic data, the joint energy concentrated dictionary (JECD) learning method yields the lowest relative square error (RSE) against the JFT basis, STVFT, and STVWT, particularly when the fraction of available samples is small and when additive noise is strong (Zhao et al., 2024).

A related extension formulates the same uncertainty geometry in a vertex-time setting and uses it for graph topology inference. In that work, the feasible region is again determined by operator eigenvalues, the dictionary atoms are eigenvectors of $(u,t_a)\in V\times T$5, and the graph-learning procedure incorporates an uncertainty-based concentration term. On a synthetic Erdős–Rényi graph, the proposed ECGL method attains correlation $(u,t_a)\in V\times T$6, adjacency error $(u,t_a)\in V\times T$7, precision $(u,t_a)\in V\times T$8, and recall $(u,t_a)\in V\times T$9, outperforming TV-GL and ESA-GL (Zhao et al., 3 Feb 2026).

5. Filtering, wavelets, and adaptive estimation

Beyond global spectral transforms, TVGS processing includes multiresolution and physics-informed constructions. Dynamic Graph Wavelets (DGW) are wavelet frames whose joint time-vertex evolution follows a dynamic process, including heat diffusion and wave propagation. For graph-based waves,

$(u,t_{p-1},v,t_0)$0

and in the joint spectral basis the wave kernel is

$(u,t_{p-1},v,t_0)$1

with stability condition $(u,t_{p-1},v,t_0)$2. A practically used causal damped wave kernel is

$(u,t_{p-1},v,t_0)$3

DGWs form a frame under explicit lower and upper bounds on $(u,t_{p-1},v,t_0)$4, and sparse recovery is posed as

$(u,t_{p-1},v,t_0)$5

solved by FISTA. On New Zealand seismic data, the method estimated earthquake epicenters with low localization error, including examples such as $(u,t_{p-1},v,t_0)$6 km at $(u,t_{p-1},v,t_0)$7 dB SNR (Grassi et al., 2016).

Optimal filtering in the joint domain can also be written as a bivariate polynomial in time and graph Laplacians: $(u,t_{p-1},v,t_0)$8 The MSE-optimal coefficients satisfy a time-vertex Wiener-Hopf equation,

$(u,t_{p-1},v,t_0)$9

and in fractional domains,

$t_{p-1}>t_0$0

Experiments report that maximum output SNR is not achieved at ordinary settings $t_{p-1}>t_0$1 or $t_{p-1}>t_0$2, but at intermediate fractional orders, and that the optimal fractional Fourier time-vertex filter outperforms the optimal time-vertex graph filter in ordinary domains and the optimal static graph filter in fractional domains, especially at low input SNR $t_{p-1}>t_0$3 dB (Ge et al., 2022).

A Tikhonov regularization formulation in the JFRT domain uses a joint fractional Laplacian and yields the closed-form filter

$t_{p-1}>t_0$4

On Molene, NOAA, and a dancer mesh dataset, fractional orders distinct from the ordinary transform improve denoising and clustering performance (Alikaşifoğlu et al., 2022).

Filter-bank generalizations introduce redundancy directly into the joint graph. The joint time-vertex oversampled graph Laplacian matrix supports bipartite extensions that preserve all edges of the original joint graph and enables two-channel oversampled graph filter banks with perfect reconstruction under spectral kernel conditions such as

$t_{p-1}>t_0$5

On the Minnesota Road Network, the oversampled method attains reconstruction MSE on the order of $t_{p-1}>t_0$6, and for the Yale Coat of Arms image improves SNR from $t_{p-1}>t_0$7 dB to $t_{p-1}>t_0$8 dB at $t_{p-1}>t_0$9 (Zhang et al., 14 Nov 2025).

Adaptive estimation broadens the TVGS scope beyond vertex-only signals. The AJVEE algorithm jointly estimates time-varying vertex and edge signals through a time-varying regression with a time-varying Laplacian,

$VT(H)=V\times T$0

and an ALMS-Hodge update

$VT(H)=V\times T$1

The method is reported as the first truly joint, online vertex-edge estimation algorithm and is validated on traffic and population mobility networks (Yan et al., 2022). In a related online graph-learning setting, AdaCGP estimates the graph shift operator from an adaptive time-vertex autoregressive model and achieves improvements in excess of $VT(H)=V\times T$2 for GSO estimation compared to adaptive vector autoregressive baselines, together with near-perfect precision after optimizing the forecast error (Jenkins et al., 2024).

6. Sampling and reconstruction theory

Sampling theory for TVGS starts from joint bandlimitedness. In finite settings, a signal can be general bandlimited (GBL), with $VT(H)=V\times T$3 nonzero JFT coefficients, or simultaneously bandlimited (SBL), with graph bandwidth $VT(H)=V\times T$4 and time bandwidth $VT(H)=V\times T$5, satisfying

$VT(H)=V\times T$6

A critical sampling set $VT(H)=V\times T$7 is a qualified sampling set with $VT(H)=V\times T$8, $VT(H)=V\times T$9, and $VT(H)$0. If $VT(H)$1 is qualified, then $VT(H)$2, $VT(H)$3, and $VT(H)$4. The theory therefore suggests assigning heterogeneous sampling pattern for each node in a network under the constraint of minimum resources (Yu et al., 2019).

For continuous correlated signals modeled as continuous time-vertex graph signals, graph correlation can explicitly reduce the minimum sampling rate. In the equal-bandwidth case with bandwidth $VT(H)$5, if $VT(H)$6 graph frequencies have zero bandwidth after graph Fourier projection, the minimum rate is

$VT(H)$7

which is lower than the naive independent rate $VT(H)$8. The feasible scheme uses dimension reduction and uniqueness sets to sample only the nonredundant joint components (Ni et al., 2022).

A more general critical sampling theory for jointly bandlimited TVGS derives lower bounds directly from the measure of joint spectral support. For continuous-time TVGS with joint bandwidth $VT(H)$9,

$X \in \mathbb{C}^{N \times T}$00

for infinite discrete-time TVGS,

$X \in \mathbb{C}^{N \times T}$01

and for finite TVGS,

$X \in \mathbb{C}^{N \times T}$02

Constructive multiband schemes show these bounds are achievable. In a synthetic $X \in \mathbb{C}^{N \times T}$03 FTVGS example with joint bandwidth $X \in \mathbb{C}^{N \times T}$04, the multiband sampling procedure uses $X \in \mathbb{C}^{N \times T}$05 samples out of $X \in \mathbb{C}^{N \times T}$06 and exactly recovers the original signal (Sheng et al., 29 Aug 2025).

Fractional spectral analysis yields a separate sampling theory. In the JFRFT-based framework, a TVGS is jointly supported and bandlimited if

$X \in \mathbb{C}^{N \times T}$07

Perfect localization and recovery hold if and only if

$X \in \mathbb{C}^{N \times T}$08

with recovery operator

$X \in \mathbb{C}^{N \times T}$09

Sampling-set selection strategies include MaxSigMin, MinTrac, MinPinv, MaxSig, and MaxVol, and localized sampling operators are introduced for large-scale processing. On seasonal US climate data with $X \in \mathbb{C}^{N \times T}$10 and $X \in \mathbb{C}^{N \times T}$11, the reported NMSE is as low as $X \in \mathbb{C}^{N \times T}$12 (Zhang et al., 22 May 2025).

When the spectral support is unknown and full access to all rows and columns is impractical, subset random sampling selects subsets of rows $X \in \mathbb{C}^{N \times T}$13 and columns $X \in \mathbb{C}^{N \times T}$14, forms a submatrix $X \in \mathbb{C}^{N \times T}$15, and then samples entries within that submatrix. A sufficient condition for high-probability recovery includes

$X \in \mathbb{C}^{N \times T}$16

and

$X \in \mathbb{C}^{N \times T}$17

On METR-LA traffic data, the proposed subset method achieves NRMSE $X \in \mathbb{C}^{N \times T}$18, $X \in \mathbb{C}^{N \times T}$19, and $X \in \mathbb{C}^{N \times T}$20 at total sampling ratios $X \in \mathbb{C}^{N \times T}$21, $X \in \mathbb{C}^{N \times T}$22, and $X \in \mathbb{C}^{N \times T}$23, outperforming SVT, TNNR, and LIMC under subset random sampling (Sheng et al., 2024).

A later reconstruction framework combines low-rank, sparsity, and smoothness priors (LSSP) under subset random sampling: $X \in \mathbb{C}^{N \times T}$24 subject to graph/time spectral representations and observation constraints. On METR-LA, with total sampling $X \in \mathbb{C}^{N \times T}$25, the reported NRMSE is approximately $X \in \mathbb{C}^{N \times T}$26 (Sheng et al., 29 Aug 2025).

7. Empirical domains, conceptual distinctions, and current directions

The empirical scope of TVGS is broad. Reported datasets include weather data from Molene and NOAA, traffic data from Sacramento and METR-LA, epidemic simulations, COVID-19 county case counts, seismological recordings from New Zealand, dancer mesh motion data, EEG signals, Minnesota Road Network epidemic simulations, rat-heart ventricular fibrillation recordings, and image/video benchmarks such as Yale Coat of Arms, REDS, and GoPro (Loukas et al., 2016). Across these domains, the recurring tasks are denoising, interpolation, forecasting, clustering, source localization, graph learning, and signal restoration.

Several conceptual distinctions recur in the literature. First, TVGS does not necessarily imply that the graph topology changes over time. Much of the signal-processing literature uses a fixed spatial graph and a temporal cycle or ring graph, whereas the TVG formalism represents genuinely time-varying topology through dynamic edges in $X \in \mathbb{C}^{N \times T}$27 (Wehmuth et al., 2014). Second, joint stationarity is not merely the conjunction of separate time and graph stationarity; the JWSS framework is strictly more general than product-graph stationarity because it allows arbitrary structure in the joint spectral plane $X \in \mathbb{C}^{N \times T}$28 (Loukas et al., 2016). Third, exact dual localization in vertex-time and spectral-frequency domains is exceptional: perfect simultaneous localization requires an eigenvalue-$X \in \mathbb{C}^{N \times T}$29 condition for $X \in \mathbb{C}^{N \times T}$30, so practical methods usually pursue maximal concentration rather than exact concentration in both domains (Zhao et al., 2024).

The current direction of the field is toward adaptive, uncertainty-aware, and structurally constrained joint processing. Recent developments include uncertainty-guided topology inference, oversampled edge-preserving filter banks, multiple-parameter fractional transforms with time-varying orders, adaptive GSO estimation, and subset sampling schemes that accommodate structural missingness and unknown spectral support (Zhao et al., 3 Feb 2026). This suggests an increasingly tight integration of harmonic analysis, stochastic modeling, sparse recovery, and graph learning within the common TVGS framework.