Stationary time-vertex signal processing

This paper considers regression tasks involving high-dimensional multivariate processes whose structure is dependent on some known graph topology. We put forth a new definition of time-vertex wide-sense stationarity, or joint stationarity for short, that goes beyond product graphs. Joint stationarity helps by reducing the estimation variance and recovery complexity. In particular, for any jointly stationary process (a) one reliably learns the covariance structure from as little as a single realization of the process and (b) solves MMSE recovery problems, such as interpolation and denoising, in computational time nearly linear on the number of edges and timesteps. Experiments with three datasets suggest that joint stationarity can yield accuracy improvements in the recovery of high-dimensional processes evolving over a graph, even when the latter is only approximately known, or the process is not strictly stationary.


Introduction
One of the main challenges when modeling multivariate processes is to decouple the estimation variance from the problem size. Consider an N-variate process unfolding over T timesteps. If only mild assumptions are made, then the number of realizations needed to reliably estimate the first two moments is up to a logarithmic factor proportional to O(NT), i.e., the data size [1]. Assuming that the process is time wide-sense stationarity (TWSS) makes the length T of the process inconsequential. This is ideal for the univariate setting as it enables us to make relevant predictions even based on a single realization. If one additionally assumes that the signal autocorrelation is compactly supported, such that most data dependencies take place within a short time horizon, then the estimation variance can be reduced further by (roughly) splitting the observations into parts and considering each as an independent realization. This approach suffices when N is relatively small. For high-dimensional processes, however, one needs to incorporate additional assumptions to obtain meaningful predictions [2][3][4][5].
In this spirit, this paper focuses on high-dimensional processes that are supported on the vertex set and are *Correspondence: andreas.loukas@epfl.ch † Andreas Loukas and Nathanaël Perraudin contributed equally to this work. 1 Laboratoire de Traitement des Signaux 2, École Polytechnique Fédérale Lausanne, 1015 Lausanne, Switzerland Full list of author information is available at the end of the article statistically dependent on the edge set of some known graph topology. Whether examining epidemic spreading [6], how traffic evolves in the roads of a city [7], or neuronal activation patterns present in the brain [8], many of the high-dimensional processes one encounters are inherently constrained by some underlying network. This realization has been the driving force behind recent efforts to re-invent classical models by taking into account the graph structure, with advances in many problems, such as denoising [9] and semi-supervised learning [10,11], among others. Yet, standard models for processes (evolving) on graphs often fail to produce useful results when applied to real datasets. One of the main reasons for this shortcoming is that they model only a limited set of spatiotemporal behaviors. The well-used graph Tikhonov and total variation priors, for instance, assume that the signal varies slowly or in a piece-wise constant manner over edges, without specifying any precise relations [12][13][14]. Similarly, assuming that the graph Laplacian encodes the conditional correlations of variables, as is done with Gaussian Markov random fields [15], becomes a rigid model when the graph is known [16]. To capture the behavior of complex networked systems, such as transportation and biological networks, it is crucial to train expressive models, being able to reproduce a wide range of graph and temporal behaviors.

Contributions
This paper considers the statistical modeling of processes evolving on graphs. In particular, we investigate the relationship between two different hypotheses: TWSS and VWSS [17][18][19], which are individually helpful in reducing the variance of covariance estimation for time series and graph signals, respectively. We propose a combined multivariate hypothesis that we refer to as time-vertex wide-sense stationarity, or joint stationarity for short. The necessary first step of our analysis consists of reformulating the standard properties of stationarity (such as the relation of the covariance matrix and power spectral density to an appropriate Fourier transform) from the lens of time-vertex analysis [20,21]. This analysis is purposeful, yet not trivial as joint stationarity is more complicated than assuming stationarity on the product of two graphs 1 .
We use the hypothesis of joint stationarity to control variance and computational complexity in estimation and recovery tasks. Similar to [17], also here, one may reliably estimate the model parameters from few observations (e.g., see Fig. 4) and solve MMSE recovery problems in time linear on the number of edges and timesteps (e.g., see Fig. 3). Complimenting previous work, we also provide an analysis of the power spectral density (PSD) estimation, which brings insight into the inherent tradeoff between bias and variance. In addition, we experimentally demonstrate that assuming joint stationarity aids in recovery even when only an approximation of the graph is known, or the process is only approximately jointly stationary. These experiments corroborate that the joint stationarity hypothesis is a useful assumption, particularly in situations when the problem features a large number of variables but only a limited number of observations.
To test the utility of joint stationarity, we apply our methods on three diverse datasets: (a) a meteorological dataset containing the hourly temperature of 32 weather stations over 1 month in Molene, France [18], (b) a traffic dataset depicting high-resolution daily vehicle flow of 4 weekdays in the highways of Sacramento, and (c) simulated SIRS-type epidemics over Europe. Our experiments confirm that for high-dimensional processes evolving over graphs, assuming joint stationarity yields an improvement in recovery performance as compared to time-or vertex-based stationarity methods, even when the graph is only approximately known and the data violate the strict conditions of our definition.

Related work
There exists an extensive literature on multivariate stationary processes, developing the original work of Wiener et al. [22,23]. The reader may find interesting Bloomfield's book [24] focusing on spectral relations. We focus on two main approaches that relate to our work, graphical models and signal processing on graphs.

Graphical models
In the context of graphical models, multivariate stationarity has been used jointly with a graph in the work of [25,26]. Though relevant, we note that there is a key difference of these models with our approach: we assume that the graph is given, whereas in graphical models, the graph structure (or more precisely the precision matrix) is learned from the data. Knowing the graph allows us to search for more involved relations between the variables. As such, we are not restricted to the case that the conditional dependencies are given by the graph (and therefore that they are sparse), but allow non-adjacent variables to be conditionally dependent, modeling a broader set of behaviors. We also note that our approach is eventually more scalable. We refer to [16] for elements of connections between graphical models and graph signal processing.

Graph signal processing
The idea of studying the stationarity of a random vector w.r.t. a graph was first introduced in [18,27] and then in [17,19]. While these contributions have different starting points, they both roughly propose the same definition. Another more recent contribution relating to stationarity on graphs in the context of PSD estimation is [28]. Despite the relevance of these works, it is important to stress that the current paper is the first to consider a stationary hypothesis over graph signals varying in time. Moreover, the new results are non-trivial as they cannot be obtained by applying previous definitions on a product graph. In addition, some of the analysis presented here (particularly that of Section 4) is novel and can also be employed for the previously studied case of stationary graph signals. To make the connection with previous works transparent, in the following, every technical result (e.g., Lemma, Theorem, Proposition) that emerges as a generalization of [18,19,27] contains a reference in its heading pointing to the former claim.
The forecasting of time-evolving signals on graphs was also considered in [29][30][31][32]. Nonetheless, there are several differences with these works, with the most important being that we define joint stationarity and that we are not restricted to the causal case (where a process is reconstructed only from its past). Finally, it should be noted that some preliminary results of this work appeared in a conference paper [33]. This work extends the conference paper in many directions. We refine the definition of joint stationarity and explore how it relates to other well-known stationarity hypotheses. We propose a new PSD estimator and provide a theoretical analysis of the bias and variance of the old and new PSD estimators. Additionally, we study the complexity of the proposed solution and evaluate its merit w.r.t. two new datasets.

General notation
We use boldface symbols for matrices and vectors (e.g., A and a respectively) and calligraphic symbols for sets (e.g., V and E). Symbol j denotes the imaginary unit, I N is the N × N identity matrix, and 1 N is the all-ones vector of size N. We use brackets to index matrix elements and subscripts for matrix blocks: if A is of size n 1 × n 2 , then A [n 1 , n 2 ] is the element at the n 1 th row and n 2 th column and A n 1 ,n 2 is a (block) matrix. Vector a = vec(A) (without subscript) is the vectorized representation of A and a n is its nth column. Moreover, A is its transpose and A * is its transposed complex conjugate (meaning that A * T is the complex conjugate). If A is N × N Hermitian, its eigendecomposition is generically written as

Harmonic time-vertex analysis
We consider signals supported on the vertices V = {v 1 , v 2 ,. . . , v N } of a weighted undirected graph G= (V, E, W G ), with E the set of edges of cardinality E = |E| and W G the weighted adjacency matrix. Suppose that signal x t is sampled at T successive regular intervals of unit length. A real time-vertex signal X = [x 1 , x 2 , . . . , x T ] ∈R N×T is then the matrix having graph signal x t as its tth column. The frequency representation of a time-vertex signal X is given by the joint Fourier transform [14,21] (or JFT for short) with U G and U T being, respectively, the unitary graph Fourier transform (GFT) and discrete Fourier transform (DFT) matrices, whereas U * T is the complex conjugate of U T . In vector form, we have thatx = JFT{x} U * J x, where U J = U T ⊗ U G . As is often the case, we choose U G to be the eigenvector matrix of the combinatorial 2 graph is the all-ones vector of size N, and diag(W G 1 N ) is the diagonal degree matrix. Matrix U T is the eigenvector matrix of the Laplacian L T of a cyclic graph T : with With this in place,X[n, τ ] can be seen as the Fourier coefficient associated with the joint frequency [λ n , ω τ ], where λ n denotes the nth graph eigenvalue and ω τ the τ th angular frequency.
The JFT maintains a close connection with the product graph J [14,21]. The latter is the graph whose adjacency matrix is W J = W T ⊕ W G (this amounts to a Cartesian product between G and the ring graph T ). The connection is revealed if one realizes that the Laplacian L J = L T ⊕ L G of J carries the eigendecomposition L J = U J ( T ⊕ G )U J . It follows that computing the JFT (in vector form) is the same as computing the GFT of x w.r.t. graph J . The main issue with any 3 product graph interpretation is that it imposes a strict dependence between the eigenvalues of L G and L T (since the eigenvalues of L J are given by T ⊕ G ). As we will see in the next paragraph, to attain full generality, one needs to abandon the product graph. For an indepth discussion of JFT and its properties, we refer the reader to [34].

Joint time-vertex filtering
Filtering a time-vertex signal x with a joint filter h(L G , L T ) corresponds to element-wise multiplication in the joint frequency domain [λ, ω] by a function h : [34][35][36]. When a joint filter h(L G , L T ) is applied to x, the output is where G ∈ R N×N and ∈ R T×T are diagonal matrices with G [n, n] = λ n and [τ , and diag(vec(A)) creates a matrix with diagonal elements the vectorized form of A. The bi-variate notation h(·, ·) is meant to illustrate that joint filters operate independently on the two domains, something impossible 4 in the product graph framework [14,21]. For convenience, we will often overload notation and write h(θ n,τ ) to refer to the bivariate function h(λ n , ω τ ). Furthermore, we say that a joint filter is separable, if its joint frequency response h can be written as the product of a frequency response h 1 defined solely in the vertex domain and one h 2 in the time domain, i.e., h(θ) = h 1 (λ) · h 2 (ω).

Joint time-vertex stationarity
Let X ∈ R N×T be a real discrete periodic multivariate stochastic process with a finite number of timesteps T that is indexed by the vertex v i of graph G and time t. We refer to such processes as time-vertex processes, or joint processes for short. Our objective is to provide a definition of stationarity that captures statistical invariance of the first two moments of a joint process x = vec(X) ∼ D(x, ), i.e., the meanx = E [x] and the covariance = E [xx ] −xx . Crucially, the definition should do so in a manner that is faithful to the graph and temporal structure.

Definition
Typically, wide-sense stationarity is thought of as an invariance of the two first moments of a process w.r.t. translation. For the first moment, things are straightforward: stationarity implies a constant mean E [x] = c1, independently of the domain of interest. The second moment, however, is more complicated as it depends on the exact form translation takes in the particular domain. Unfortunately, for graphs, translation is a non-trivial operation and three alternative translation operators exist: the generalized translation [37], the graph shift [13], and the isometric graph translation [27]. Due to this challenge, there are currently three alternative (though akin) definitions of stationarity appropriate for graphs [17][18][19].
The ambiguity associated with translation on graphs urges us to seek an alternative starting point for our definition. Fortunately, there exists an interpretation which holds promise: up to its constant mean, a wide-sense stationary process corresponds to a white process filtered linearly on the underlying space. This "filtering interpretation" of stationarity is well known classically 5 as well as in the graph setting [19] and is equivalent to asserting that the second moment can be expressed as = h(L T ), where h(L T ) is a linear filter. Thankfully, not only filtering is elegantly and uniquely defined for graphs [37], but 4 Defining joint filters in terms of a product graph would imply that there is a fixed relation between angular and graph frequencies (determined by the product graph construction). As a result, in the product graph, framework filters are univariate functions. 5 As the correlation between two instants t 1 and t 2 depends only on the difference between these two instants , the covariance matrix has to be circulant, a property that is shared by linear filters. also stating that a process is graph wide-sense stationary if E[ x] = c1 N and = h(L G ), is a graph filter, is generally consistent 6 with current definitions [17][18][19].
This motivates us to also express the definition of stationarity for joint processes in terms of joint filtering:

Definition 1 (JWSS) A joint process x = vec(X) is called jointly wide-sense stationary (JWSS), if and only if
(a) The first moment of the process is constant is a non-negative real function referred to as joint power spectral density (JPSD).
Let us examine Definition 1 in detail. First moment condition. As in the classical case, the first moment of a JWSS process has to be constant over the time and the vertex sets, i.e.,X[i, t] = c for every i = 1, 2, . . . , N and t = 1, 2, . . . , T. For alternative choices of the graph Laplacian with a null-space not spanned by the constant vector, the first moment condition should be modified to requiring that the expected value of a JWSS process is in the null space of the matrix L T ⊕ L G (see Remark 2 [19] for a similar observation on stochastic graph signals).
Second moment condition. According to the definition, the covariance matrix of a JWSS process takes the form of a joint filter h(L G , L T ), and is therefore diagonalizable by the JFT matrix U J . It may also be interesting to notice that the matrix h(L G , L T ) can be expressed as follows where each block H t 1 ,t 2 of is an N ×N matrix defined as: and h ω τ (L G ) is the graph filter with frequency response h ω τ = h(λ, ω τ ). Being a covariance matrix, h(L G , L T ) must necessarily be positive-semidefinite; thus, h(·, ·) is real (the eigenvalues of every Hermitian matrix are real) and non-negative. Also, equivalently, every zero mean JWSS process x = vec(X) can be generated by joint filtering x = h(L G , L T ) 1/2 a white process with zero mean and identity covariance. The following proposition exploits these facts to provide an interpretation of JWSS processes in the joint frequency domain.
(For clarity, this and other proofs of the paper have been moved to the "Appendix".) We briefly present a few additional properties of JWSS processes that will be useful in the rest of the paper.
The proof follows easily by noting that the covariance of w is diagonalized by the joint Fourier basis of any graph This last equation tells us that the JPSD is constant, which implies that similar to the classical case, the energy of white noise is evenly spread across all joint frequencies.
A second interesting property of JWSS processes is that stationarity is preserved through a filtering operation.

Property 2 (Generalizes Theorem 2 [17], Property 1 [19]) When a joint filter f
Finally, we notice that for real processes X, which are the focus of this paper, the function h forming the joint filter should be symmetric w.r.t. ω, meaning that h(λ, ω) = h(λ, −ω). This property can be easily derived from the definition of the Fourier transform.

Relations to classical definitions
We next provide an in-depth examination of the relations between joint wide-sense stationarity, time and vertex stationarity, as well as their multivariate equivalents. For clarity, we order the rows/columns of the covariance matrix such that each t 1 ,t 2 block of size N × N measures the covariance between x t 1 and x t 2 (see (4)).

Standard definitions
As we discuss below, known definitions of stationarity in time/vertex domains are particular cases of joint stationarity.
TWSS ∩ VWSS ⊂ JWSS. The known versions of stationarity (TWSS, VWSS) are oblivious to any structure along one of the two dimensions of X. In this manner, assuming that X is TWSS amounts to interpreting each of the N time series as a separate realization of the same process with TPSD h T (ω). Similarly, if X is VWSS, then each graph signal x t is taken as a separate realization of a single stochastic graph signal with VPSD h G (λ) [17,19]. It is a simple consequence that, different from the JWSS hypothesis, assuming that X is both TWSS and VWSS is equivalent to limiting our scope to separable JPSD defined as the product of two univariate functions

Definitions based on the product graph
As explained in Section 2, the JFT can be interpreted as a graph Fourier transform taken over a product graph whose Laplacian is L J = L G ⊕ L T . This construction can give rise to two additional definitions for joint stationarity: VWSS on a product graph. The first is obtained by applying the VWSS definition of [17,19] on the graph associated with L J . The resulting model is not sufficiently general in order to generate the full spectrum of JWSS processes. The reason is that, whereas the JPSD h(λ, ω) can be any two-dimensional non-negative function, the JPSD of any VWSS process on L J is necessarily onedimensional (the eigenvalues of L J are the sums of all combinations of the eigenvalues of L G and L T )-see Fig. 1 for a pictorial demonstration and "Appendix: Univariate vs multivariate JPSD" for examples from real data. The same reasoning also holds for alternative products between graphs, such as the strong and Kronecker products [14].
Covariance diagonalized by the product graph Fourier transform. The second definition, which we refer to as JWSS-alternate, entails asserting that the covariance matrix can be diagonalized by the JFT, i.e., the eigenbasis of L J . This can be seen to differ from the JWSS definition only in case of graph Laplacian eigenvalue multiplicities: whenever the graph Laplacian features repeated eigenvalues, for Definition 1, the degrees of freedom of the JPSD h decrease, as necessarily h(λ 1 , ω) = h(λ 2 , ω) when λ 1 = λ 2 . This restriction is motivated by the following observation: for an eigenspace with multiplicity l l l l l Fig. 1 The joint stationarity hypothesis is more general than assuming either (standard) VWSS and TWSS or VWSS on a (Cartesian) product graph. The figure presents three examples of PSDs plotted as 2-dimensional function of λ 1 , λ 2 that, for simplicity, corresponds to the eigenvalues of two graphs. The second graph (time) is a ring. In the separable case (left), the PSD has to satisfy h(λ 1 , , making it unable to capture any dependencies between λ 1 and λ 2 . Using VWSS (middle) limits the PSD to h(λ 1 , λ 1 ) = h(λ 1 + λ 2 ) leading to constant values along the diagonal line λ 1 + λ 2 = c. Joint stationarity (right) can encode any PSD h(λ 1 , λ 2 ), as exemplified here greater than one, there exists an infinite number of possible eigenvectors corresponding to the different rotations in the space, and the JPSD is in general ill-defined. The condition h(λ 1 , ω) = h(λ 2 , ω) when λ 1 = λ 2 deals with this ambiguity, as it ensures that the JPSD is the same independently of the choice of eigenvectors. On the contrary, with JWSS-alternate, one should construct an arbitrary basis of each eigenspace with multiplicity and set 7 h(λ 1 , ω) = h(λ 2 , ω). This approach, which was followed in [38], features more degrees of freedom at the expense of the loss of filtering interpretation and higher computational complexity: one may not anymore use filters to estimate the JPSD (without reverting to Definition 1), whereas using the JFT to diagonalize the covariance scales like O N 3 + N 2 T + NT log(T) . On the contrary, in our setting, the PSD estimation complexity can be reduced to be close to linear in the number of edges E and timesteps T (see "Appendix: Implementation details of the JPSD estimator"). Nevertheless, we should mention that the differences mentioned above are mostly academic. Eigenvalue multiplicities occur mainly when graph automorphisms exist. In the absence of such symmetries (e.g., in the graphs used in our experiments), the two definitions yield the same outcome.

Multivariate definitions
On the other hand, joint stationarity can itself be derived as the combination of two multivariate versions of time/vertex stationarity, which we refer to respectively as MTWSS (see [25]) and MVWSS. Before formally defining them in Definitions 2 and 3, let us state our result formally:

Theorem 1 A joint process X is JWSS if and only if it is MTWSS and MVWSS.
To put this in context, we examine the two multivariate definitions independently.
(a) JWSS ⊂ MTWSS. The covariance matrix of a JWSS process has a block circulant structure, as t 1 ,t 2 = δ,1 = δ , where δ = t 1 − t 2 + 1. Hence, can be written as implying that correlations only depend on δ and not on any time localization. This property is shared by multivariate time wide-sense stationary processes: Similarly to the univariate case, the time power spectral density (TPSD) is defined to encode the statistics of the process in the spectral domain We then obtain the TPSD of a JWSS process by constructing a graph filter from h while fixing ω. Setting h ω τ (λ) = h(λ, ω τ ), the TPSD of a JWSS process isˆ τ = h ω τ (L G ). (b) JWSS ⊂ MVWSS. For a JWSS process, each block of has to be a linear graph filter, i.e., t 1 ,t 2 = γ t 1 ,t 2 (L G ), meaning that This is perhaps better understood when compared to the multivariate version of vertex stationarity defined below: It can be seen that every JWSS process must also be MVWSS, or equivalently JWSS ⊂ MVWSS.

Joint power spectral density estimation
The joint stationarity assumption can be useful in overcoming the challenges associated with dimensionality. The main reason is that for JWSS processes, the estimation variance is decoupled from the problem size. Concretely, suppose that we want to estimate the covariance matrix of a joint process x = vec(X) from K samples x (1) , x (2) , . . . , x (K) . As we show in the following, if the process is JWSS such that = h(L G , L T ), the JPSD estimation variance is O(1). This is a sharp decrease from the classical and MTWSS settings, for which K ≈ NT and K ≈ N realizations are necessary 8 , respectively.
This section presents two JPSD estimators. The first provides unbiased estimates at a complexity that is O N 3 T log(T) . The second estimator decreases further the estimation variance at the cost of a bounded bias and is approximated with (close to) linear complexity.

Sample JPSD estimator
We define the sample JPSD estimator for every graph frequency λ n and angular frequency ω τ as the estimatė In case the process does not have zero mean, it should be centered by subtracting the constant signal c 1 N 1 * T , where 8 The number of realizations needed for obtaining a good sample covariance matrix of an n-dimensional process is O(n log n) [1,39].
In that case, the unbiased estimator should involve division by K − 1, instead of K as we have in (7).

Analysis
For simplicity, in the following, we suppose that the process is correctly centered. As the Theorem 2 claims, the sample JPSD estimator is unbiased, and its variance decreases linearly with the number of samples K.

Theorem 2 For every distribution with bounded second and fourth order moments, the sample JPSD estimatorḣ(θ)
(a) is unbiased, i.e., E ḣ (θ) = h(θ), and where constant γ depends only on the distribution of x.
Proof For any θ =[λ, ω], the sample estimate iṡ withε (k) being independent realizations ofε, a zero mean complex random variable with unit variance. To see this, write x = h (L G , L T ) 1/2 , where the random vector has zero mean and identity covariance. Then, the complex random variableε is the JFT coefficient of corresponding to frequencies λ and ω. The bias follows by noting that E ˆ (k)ˆ * (k) = 1, for every k. The variance is computed similarly by exploiting the fact that different terms in the sum are independent as they correspond to distinct realizations and setting γ = E |ε| 4 .
For the standard case of a Gaussian joint process, we provide an exact characterization of the distribution.

Corollary 1 For every Gaussian JWSS process, the sample JPSD estimate follows a Gamma distribution with shape K/2 and scale 2h(θ)/K. The estimation error variance is equal to Var ḣ (θ) = 2 h 2 (θ)/K.
Proof We continue in the context of the proof of Theorem 2. For a Gaussian distribution,ε is centered and scaled Gaussian and thusε 2 is a chi-squared random variable with 1 degree of freedom. Our estimate is, therefore, a scaled sum of i.i.d. chi-squared variables and corresponds to a Gamma distribution. The corollary then follows directly.
Observe that the variance depends linearly on the fourth-order moment of |ε| (see proof of Theorem 2) and is inversely proportional to the number of samples, but it is independent of N and T. This implies that || −˙ 2 || can be made arbitrarily small using K = O(1) samples. In the following, we discuss how to achieve an even smaller variance by exploiting the properties of h(θ).

Convolutional JPSD estimator
When the number of available realizations K is small (even 1), one may make use of additional assumptions on to obtain reasonable estimates. To this end, we next present a parametric JPSD estimator that allows us to trade off bias for variance.
Before delving into JWSS processes, it is helpful to consider the purely temporal case. For a TWSS process, it is customary to assume that the autocorrelation function has support L that is a few times smaller than T. Then, cutting the signal into T L smaller parts and computing the average estimate reduces the variance (by a factor of T L ), without sacrificing frequency resolution. This basic idea stems from two established methods used to estimate the PSD of a temporal signal, namely Bartlett's and Welch's method [40,41]. Averaging across different windows is equivalent to smoothing the TPSD by convolving it with a window in the frequency domain: this results in attenuation of the correlation for long delays, enforcing localization in the time domain.

Estimator
Armed with this interpretation, we proceed by smoothing the JPSD with a user-specified bi-variate window g, such as a Gaussian or a disc window. The convolutional JPSD estimator computes the JPSD at joint frequency θ = (λ, ω) as: where c g (θ) n,τ g(θ − θ n,τ ) 2 is a normalization factor. For implementation specifics, including a discussion on the choice of the bivariate kernel g, we refer the reader to "Appendix: Implementation details of the JPSD estimator".
The convolutional JPSD estimator is related to known PSD estimators for TWSS and VWSS processes. The Dirac function is denoted by φ. We have that (a) for g(θ) = φ(λ) · g T (ω), we recover the classical TPSD estimator, applied independently for each λ. (b) For g(θ) = g G (λ) · φ(ω), we recover the VPSD estimator from [17] applied independently for each ω. Similar to the latter, the estimator can be closely approximated at a complexity that is linear w.r.t. the number of graph edges/nodes, and up to a logarithmic factor linear to the number of timesteps (see "Appendix: Implementation details of the JPSD estimator").

Analysis
To provide a meaningful bias analysis, we introduce a Lipschitz continuity assumption on the JPSD, matching the intuition that localized phenomena tend to have a smooth representation in the frequency domain.

Theorem 3 At θ, the convolutional JPSD estimatorḧ(θ)
(a) has bias where is the Lipschitz constant of h(θ), and (b) when the entries ofX are independent random variables, its variance is where Var ḣ (θ n,τ ) is the variance of the sample JPSD estimator at θ n,τ .
The derivations of the bias and variance are given in Lemmas 1 and 2, respectively.
We note two corner cases of interest. In the most convenient case, the JPSD is constant, and our estimator is unbiased (the Lipschitz constant is zero). On the other hand, if the JPSD fluctuates rapidly, the bias of the estimate will be significant unless g is close to a Dirac. Here, the sample estimator should be preferred.
We further consider as a theoretical example the case of a Gaussian JWSS process and a (spectral) disc window with bandwidth B, i.e., g B (θ) = 1 if ||θ|| 2 ≤ B 2 and 0 otherwise. Though perhaps not the most practical choice from a computational perspective, we consider here a disc window because it leads to simple and intuitive estimates.
Proof The results follow from Theorem 3 and Corollary 1 by noting that when a disc window is used, (a) c g (θ) = |S θ | and (b) g(θ − θ n,τ ) 2 = 1 for all n, τ in the window (there are |S θ | in total) and zero otherwise. The independence condition required by the variance clause of the theorem is satisfied sincex is Gaussian (as a rotationx = U * J x of a Gaussian vector) with diagonal covariance. The above result suggests that by selecting our window (bandwidth), we can trade off bias for variance. The trade-off is particularly beneficial as long as (a) the JPSD is smooth relatively to the disc size ( B 1) and (b) the graph eigenvalues are clustered (|S θ | 1 when h(θ) 0).

Recovery of JWSS Processes
This section considers the MMSE problem of recovering a JWSS process x = vec(X) from linear measurements y corrupted by a zero-mean JWSS process w: where the function f is linear on y, i.e., there exists a matrix W and a vector b such that f (y) = Wy + b. We remark that (a) for A binary diagonal and w = 0, (P0) is an interpolation problem, (b) for A = I and w white noise (P0) is a denoising problem, and (c) for A diagonal with A ii = 1 if i ≤ Nt and zero otherwise and w = 0 it corresponds to forecasting. We mainly consider the former two problems since, for forecasting, it is more computationally efficient to utilize autoregressive models [29].
The minimum mean-squared linear estimate is known to bė with the definitions y = A A * + w and xy = A * . Obtainingẋ therefore entails solving a linear system in matrix y that-naively approached-has O N 2 T 2 complexity. In addition, the condition number of y can be large, rendering direct inversion unstable. For instance, this may happen when one attempts to reverse any smoothing operation A that severely attenuates part of the signal's spectrum. We next discuss how to deal with these issues:

Decreasing the complexity
Thankfully, even if y is not always sparse, we can approximate its multiplication by a vector without actually computing it as (a) A is, for many applications (denoising, prediction, forecasting), sparse, and (b) per our assumption, and w are joint filters, and therefore, they can be implemented at complexity that is (up to logarithmic factors) linear to the number of edges E and timesteps T [20,34,36]. Therefore, if we employ an iterative method such as the (preconditioned) conjugate gradient to compute the solution, the complexity of each iteration will be linear on the problem size.

Singular or badly conditioned y
We choose the solution with the minimal residual by substituting the inverse −1 y in (11) with the pseudo-inverse + y . However, instead of solving the normal equationsẋ = xy 2 y −1 y (y −ȳ) +x, which has the effect of significantly increasing the condition number of our matrix, we suggest to employ the minimal residual conjugate gradient method for symmetric matrices [42]. For badly conditioned covariance matrices, an alternative solution is to rewrite the problem as a regularized least squares problem and solve it using the generalization of the fast iterative shrinkage-thresholding algorithm (FISTA) scheme [43][44][45]. This problem was shown to converge to the correct solution when w is white noise. More details about the optimization procedures can be found in [17]. Similarly, in the noiseless case one removes term Az − y 2 2 in (12) and introduces instead the constraint Az = y. The resulting optimization problem can be solved using a Douglas-Rachford scheme [46].

Joint power spectral density estimation
The first step in our evaluation is to analyze the efficiency of JPSD estimation. Our objective is dual. First, we aim to study the role of the different method parameters into the estimation accuracy and computational complexity, essentially providing practical guidelines for their usage. In addition, we wish to illustrate the usefulness of the joint stationarity assumption, even when the graph is only approximately known.

Variance-bias-complexity tradeoffs
To validate the analysis of Section 4 for the computational and accuracy trade-offs inherent to our JPSD estimation method, we performed numerical experiments with random geometric graphs of N = 256 vertices (we build a 10-nearest neighbor graph, weighted by a radial basis function kernel tuned so th at the average weighted degree is slightly above 7) and JWSS processes (T = 128 timesteps). Though our approach works with any JPSD, including high frequency ones, in this experiment, we consider a stochastic process generated by the discrete damped wave equation with a non-separable JPSD h(λ, ω) = exp(−|ω|/2) cos(ω acos(1 − λ))

Variance-bias
First, we examine the relation between the real JPSD h and the convolutional estimateḧ obtained using the "fast" method described in "Appendix: Implementation details of the JPSD estimator". We use the following metrics:

Loukas and Perraudin EURASIP Journal on Advances in Signal Processing
is the empirical expectation computed over 20 independent experiments. We remind the reader that there are two parameters influencing the performance of the convolutional JPSD estimator (see "Appendix: Implementation details of the JPSD estimator": the window size L corresponding to our assumption for the support length of the autocorrelation in time, and the number of graph filters F used to capture power density in the graph spectral dimension. As discussed in Theorem 3, the bias will be small as long as the JPSD is a smooth function (it has a small Lipschitz constant ), in which case one may opt for small L and F. Figure 2a-d report four key metrics for an exhaustive search of L, F combinations. We observe that large values of F and L generally reduce the estimation error (Fig. 2a) because they result in reduced bias (Fig. 2b). Nevertheless, setting the parameters to their maximum values is not suggested as the variance is increased (Fig. 2c).

Complexity
In Fig. 2d, we see that utilizing a large number of filters (i.e., large F) increases the average execution time. Figure 3 delves further into the issue of scalability. In particular, we vary the number of vertices from 1000 to 9000 and focus on a process with JPSD h(θ) = e −λ/λ max e −5 ω 2 .
We then examine the min/median/max execution time of the convolutional JPSD estimator for a for increasing problem sizes when ran in a desktop computer and repeated 10 times. We compare two implementations. The first, which naively performs the convolution in the spectral domain, uses the eigenvalue decomposition and therefore scales quadratically with the number of vertices. Due to its optimized code and simplicity, this should be the method of choice when N is small. For larger problems, we suggest using the fast implementation. As shown in the figure, this scales linearly with N (here E = O(N)) when the number of filters F and timesteps T are held constant. In this experiment, we set L to 64.

How to choose L and F?
Having no computational constrains, one should choose the parameter combination that minimized the Akaike information criterion (AIC) score AIC = 2FL − 2 ln ¨ , where¨ is the distribution dependent estimated likeli-hood¨ = P x|¨ and¨ is the estimated covariance based on the convolutional JPSD estimator with parameters L and F [47]. This procedure is often unfeasible as it is based on computing each model's log-likelihood and thus entails estimating one JPSD for each parameterization in consideration (as well as knowing the distribution type). We have found experimentally that setting F = min(N, 50) provides a good trade-off between computational complexity and error. On the other hand, we suggest setting L to an upper bound of the autocorrelation support. Figure 4 illustrates the benefit of a joint stationarity prior as compared to (a) an empirical covariance estimator which makes no assumptions about the data and (b) the MTWSS process estimator with optimal bandwidth [22]. As expected, accurate estimation is challenging when the number of realizations is much smaller than the number of problem variables (NT), returning errors above one for the empirical estimator. Introducing stationarity priors regularizes the estimation resulting in more stable estimates.

Learning from few realizations and a noisy graph
What is perhaps surprising is that even when the graph (and U G ) is known only approximately, estimating the second order moment of the distribution using the joint stationarity assumption is beneficial. To portray this phenomenon, we also plot the estimation error when using a noisy graph (we corrupted the weighted adjacency matrix by Gaussian noise, resulting in an SNR of 10 dB). Undoubtedly, introducing noise to the graph edges negatively affects estimation by introducing bias. Still, even with noise, the proposed method significantly outperforms purely time-based methods when less than NT realizations are available.

Recovery performance on three datasets
We apply our methods on three diverse datasets featuring multivariate processes evolving over graphs: (a) a weather dataset depicting the temperature of 32 weather stations over 1 month, (b) a traffic dataset depicting highresolution daily vehicle flow of 4 weekdays, and (c) SIRStype epidemics in Europe. Our experiments aim to show that joint stationarity is a useful model, even in datasets which may violate the strict conditions of our definition, and that it can yield a significant improvement in recovery performance, as compared to time-or vertex-based stationarity methods.

Experimental setup
We split the K realizations of each dataset into a training set of size p t K and a test set of size (1 − p t )K, respectively. The training set is used to estimate the JPSD. Then, in the first two experiments, we attempt to recover the values of p d NT variables randomly discarded from the test set. This corresponds to A being a binary diagonal matrix and In the third experiment, we instead consider a denoising problem with A = I and w being a random Gaussian vector. In each case, we report the RMSE for the recovered signal normalized by the 2 -norm of the original signal. We compare our joint method with the sample and convolutional JPSD estimators to univariate time/vertex stationarity [17]. These methods solve the statistical recovery problem under the assumption that signals are stationary in the time/vertex domains, but considering different vertices/timesteps as independent. These methods are known to outperform non-model based methods, such as Tikhonov regularization (ridge regression) and total-variation regularization (lasso) over the time or graph dimensions [12,13]. We also compare to the more involved MTWSS model [25] where the values at different vertices are correlated and the covariance is block circulant of size NT ×NT (see Definition 2). The latter is only shown for the weather dataset as the large number of variables present in the other datasets (e.g., ≈ 10 8 parameters for the traffic dataset) prohibited computation. We remark that the graph Laplacians we considered did not possess eigenvalue multiplicities, meaning that the results obtained using the JWSS-alternate definition are identical to that with JWSS using a sample JPSD estimator-thus, we do not include JWSS-alternate in our comparison.

Molene dataset
The French national meteorological service has published in open access a dataset 9 with hourly weather observations collected during the month of January 2014 in the region of Brest (France) [18]. The graph was built from the coordinates of the weather stations by connecting all the neighbors in a given radius with a weight function W G [i 1 , is the Euclidean distance between the stations i 1 and i 2 . Parameter k was adjusted to obtain an average degree around 5 (k however is not a sensitive parameter). We split the data in K = 15 consecutive periods of T = 48 h each. As sole pre-processing, we removed the mean (over time and stations) of the temperature 10 . We first test the influence of training set size p t , while discarding p d = 30% of the test variables. As seen in Fig. 5a, due to its large sample complexity, the MTWSS approach provides good recovery estimates when the number of realizations is large, approaching that of joint stationarity, but suffers for small training sets (though not shown in the figure, the relative mean error was 9.8 when only p t = 10% of the data was used for training). Due to their stricter modeling assumptions, univariate stationarity methods returned relevant estimates when trained from few realizations but exhibited larger bias. The convolutional JPSD estimator can be seen to improve upon the sample estimator when the amount of data used for JPSD estimation is small (less than 20%). For bigger training sets, the two estimators yield similar accuracy. Figure 5b reports the achieved errors for recovery problems with progressively larger percentage 5% ≤ p d ≤ 95% of discarded entries for a training percentage of p t = 20%. We can observe that the error trends are consistent across all cases.

Traffic dataset
The California department of transportation publishes high-resolution traffic flow measurements (number of vehicles per unit interval) from stations deployed in the highways of Sacramento 11 . We focused on 727 stations over four weekdays in the period 01-06 April 2016. Starting from the road connectivity network obtained by the OpenStreetMap.org, we constructed one time series for each highway segment by setting the flow over it to be a weighted average of all nearby stations, while abiding to traffic direction. This resulted in a graph of N = 710 vertices and a total of T = 24 × 12 measurements per day for K = 4 days. We used the convolutional JPSD estimator with parameters L = T/2 and F = 75, which were experimentally found to give good performance in the training set. Figure 6a and b depict the mean recovery errors when the training sets were 1 (p t = 25%) and 3 days (p t = 75%) respectively. The strong temporal correlations present in highway traffic were useful in recovering missing values. Considering both the temporal and spatial dimensions of the problem resulted in accurate estimates, with less than 0.04 error when p d =50% of the data were removed and the PSD was estimated from 1 day. As expected, the convolutional estimator is efficient in the case when the training set is small (1 out of 4 days used for training): assuming that the JPSD is smooth helps to reduce estimation variance and computational complexity but can lead to a slight decrease in accuracy when a large amount of training data is available.

SIRS epidemic
Our third experiment simulates the spread of an infectious disease over N = 200 major cities of Europe, as predicted by the susceptible-infected-recovered-susceptible (SIRS) model, one of the standard models used to study epidemics. We intend to examine the predictive power of the considered methods when dealing with different realizations of a non-linear and probabilistic process over a graph (the data are fictitious). We parameterized SIRS as follows: length of infection period, 2 days; length of immunity period, 10 days; probability of contagion across neighboring cities per day, 0.005; and total period, T = 180 days.
We generated a total of K = 10 infections, all having the same starting point. In contrast to the previous experiments, here, we attempt to recover the data after they have been corrupted with additive Gaussian noise. Figure 7a and b depict the mean recovery error as a function of the input signal-to-noise ratio (SNR), respectively, when p t = 50% and p t = 90% of the data were used for training. As in previous experiments, the joint stationarity attains better recovery. The difference becomes clearer for low SNR, in which case the error is decreased (roughly) by a factor of two w.r.t. the best alternative.

Conclusion
This paper proposed a new definition of wide-sense stationarity appropriate for multivariate processes supported on the vertices of a graph.
Our model presents two key benefits. First, the estimation and recovery of JWSS processes is efficient, both in terms of estimation variance and computational complexity. In particular, the JPSD of a JWSS process can be estimated from few observations at a complexity that is roughly linear to the number of graph edges and timesteps. After the PSD has been estimated, the linear MMSE recovery problems of interpolation and denoising This section describes how to approximate a convolutional estimate using a number of operations that is linear to ET. Before describing the exact algorithm, we note two helpful properties of the estimator. First, we can computë h(θ) by obtaining estimates for each X (k) independently and then averaging over k: h(θ) = 1 K c g (θ) k n,τ g(θ − θ n,τ ) 2 |JFT{X (k) }[n, τ ] | 2 As we will see in the following, the terms inside the outer sum can be approximated efficiently, avoiding the need for an expensive JFT. In addition, when the convolution window is separable, i.e., g(θ) = g G (λ)·g T (ω), as is assumed in this contribution, the joint convolution can be performed successively (and at any order) in the time and vertex domains where c g (θ) = c g T (ω) · c g G (λ). Exploiting this property, we treat the implementation of the two convolutions separately and the presented algorithms can be combined in any order.

Fast time convolution.
This is the textbook case of TPSD estimation that is solved by the Welch's method [41]. The method entails splitting each time series into equally sized overlapping segments and averaging over segments the squared amplitude of the Fourier coefficients. The procedure is equivalent to an averaging (over time) of the squared coefficients of a short-time Fourier transform (STFT), with half-overlapping windows w T defined such that DFTw T (t) = g T (ω) [51,52]. Let L be the support of the autocorrelation or equivalently the number of frequency bands. We suggest using the iterated sine window w T (t) sin 0.5π cos (πt/L) 2 if t ∈[−L/2, L/2] 0 o t h e r w i s e , as it turns the STFT into a tight operator. In order to get an estimate ofḧ at unknown frequencies, we interpolate between the L known points using splines [53].

Fast graph convolution.
Inspired by the technique of [17], we perform the graph convolution using an approximated graph filtering operation [54] that scales linearly to the number of graph edges E. In particular, . (13) We suggest using the Gaussian window with σ 2 = 2(F + 1)λ max /F 2 . As we did before, we only compute the above for F = O(1) different values of λ and approximate the rest using splines. As the eigenvalues are not known, we need a stable way to estimate c g G (λ). We obtain an unbiased estimate by filtering Q = O(1) random Gaussian signals on the graph ∈ R N ∼ N (0, I N ), such that with variance equal to 2 N n=1 g 4 (λ − λ n )/Q. We omit the analysis, as it is similar to that in Theorem 2. According to our numerical evaluation, the approximation error introduced by the latter estimator and spectral filtering is almost negligible for smooth JPSD.
O(TKF × E + QF × E) = O((TK + Q)EF) for the fast graph convolutions. Here, the TK and Q convolutions are performed in order to estimate the quantities at (13) and (15) for F different values of λ. (b) O(NK × T log(L)) for the fast time convolution, corresponding to NK STFT. Thus, in total the complexity of the fast convolutional JPSD estimator is O(TKFE + QEF + NKT log(L)). Furthermore, when Q, F, K, L are constants, the complexity simplifies to O(TE). We remark that, though asymptotically superior, the fast implementation can be significantly slower when the number of variables is small. Our experiments demonstrate that it should be preferred for N larger than a few thousands (see Fig. 3).

Univariate vs multivariate JPSD
As discussed in Section 3.2, one could potentially pose a VWSS hypothesis on a product graph to define joint stationarity, but the direct effect of such a choice is that the spectral domain becomes 1-dimensional instead of 2dimensional. To see why this is problematic, in Fig. 8, we plot the two different representations of the JPSD for the three datasets featured in our experiments. It can be seen that the 2D representation (corresponding to the JWSS hypothesis) is more structured than its 1D counterpart. More importantly, a JWSS hypothesis leads to a smoother JPSD: this is what our convolutional JPSD estimator employs to decrease the estimation variance.

Deferred proofs
Proof of Proposition 1 In order to simplify the notation in the next proof, we define the unravel function u r : Z 2 → Z that transforms the double indexes n, τ of the matrix indexing of X into its vector index of u(n, τ ) = (τ − 1)N + n, i.e., X[n, τ ] = vec(X)[u(n, τ )].
By construction of the JFT basis,X[ 0, 0] captures the DC-offset of a signal, and condition (a) is equivalent to stating that E[x] = c1 NT . Moreover, if the graph is connected and (a) holds, at least one of E X [n 1 , τ 1 ] and