Time-Varying Graph Signal Recovery Using High-Order Smoothness and Adaptive Low-rankness

A time-varying or spatial-temporal graph signal can be considered as a sequence of signals arranged chronologically, where each signal at a specific time instance is defined on a static or dynamically changing spatial graph.

Consider an undirected unweighted graph $G=(V,E),$ where $V$ is a set of $n$ vertices and $E$ is a set of edges. We assume a collection of time-varying graph signals $\{\mathbf{x}_{t}\}_{t=1,\ldots,m}$ with $\mathbf{x}_{t}\in\mathbb{R}^{n}$ are defined on $V$ with a time index $t$. Let $X=[\mathbf{x}_{1},\ldots,\mathbf{x}_{m}]\in\mathbb{R}^{n\times m}$ be the data set represented in matrix. The pairwise connections on the graph $G$ can be modeled by an adjacency matrix $A$, where the $(i,j)$-th entry of $A$ is one if there is an edge between vertices $i$ and $j$, and zero otherwise. This binary adjacency matrix can be extended to the non-binary case for a weighted graph, where each entry indicates the similarity between two vertices. Throughout the paper, we use a standard $k$ nearest neighbor (kNN) approach based on the Euclidean distance of data points to construct the adjacency matrix.

Given an adjacency matrix $A$, we further define the graph Laplacian matrix, $L=M-A\in\mathbb{R}^{n\times n},$ where $M$ is a diagonal matrix with its diagonal element ${M}_{ii}=\sum_{j}{A}_{ij}$. The graph Laplacian serves as a matrix representation of the graph structure and can be used to describe some important characteristics of a graph, such as node connectivity and similarity. For example, geographic locations in the form of coordinates, i.e., longitude and latitude, are typically used to calculate the pairwise distance and, thereby, the graph Laplacian for geospatial data. For some data sets without obvious graph domains, a preprocessing step of graph learning can be implemented; see xia2021graph for a comprehensive review of graph learning techniques.

Time-varying graph signal recovery aims to recover an underlying matrix from its partially observed entries that are possibly polluted by additive noise. Mathematically, a forward model is $Y=J\circ X+\mathcal{N}$, where $Y$ is the observed data, $J\in\{0,1\}^{n\times m}$ is a sampling matrix, and $\mathcal{N}$ is a random noise. In this work, we focus on recovering time-varying signals, represented by the matrix $X$, from incomplete noisy data $Y$ defined on static spatial graphs in the sense that the vertex set and the edges do not change over time. In addition, we adopt a symmetrically normalized graph Laplacian that is pre-computed based on geographic locations.

The recovery of graph signals from partial observations is an ill-posed problem due to missing data. Graph regularization plays a crucial role in developing a recovery model for time-varying signals by enforcing temporal correlation and/or describing the underlying graph topology. An intuitive approach for recovering time-varying graph signals is to apply interpolation methods to fill in the missing entries, such as natural neighborhood interpolation (NNI) sibson1981brief . Numerous recovery models with diverse smoothness terms have been proposed to further preserve the underlying geometry. For example, Graph Smoothing (GS) narang2013localized characterizes the smoothness of the signal using the graph Laplacian of $X$. Alternatively, temporal smoothness is incorporated in Time-Varying Graph Signal Recovery (TGSR) qiu2017time by formulating the graph Laplacian of $DX,$ where $D$ is a first-order temporal difference operator. The combination of the graph Laplacian of $X$ and the Tikhonov regularity of $DX$ was considered in perraudin2017towards . In contrast, the graph Laplacian of $DX$ with an additional low-rank regularity of $X$ was formulated as Low-Rank Differential Smoothness (LRDS) mao2018spatio . In the Tikhonov regularization, $\left\lVert XD\right\rVert_{F}^{2}=\operatorname*{tr}(XDD^{T}X^{T})$ implies that $DD^{T}$ is treated as the temporal graph Laplacian. In giraldo2022reconstruction , the graph Laplacian matrix $L$ is replaced by $(L+\epsilon I)^{r}$, where $I$ is the identity matrix and $r\geq 1$ for a high-order Sobolev spatial-temporal smoothness. Its main advantage lies in faster convergence, as this approach does not necessitate extensive eigenvalue decomposition or matrix inversion. Recently, another low-rank and graph-time smoothness (LRGTS) method has been proposed in liu2023time , where the sum of the nuclear norm and the Tikhonov regularizer on the second-order temporal smoothness is adopted to promote the low-rankness and the temporal smoothness, respectively.

All the models mentioned above can be condensed into one minimization framework:

where $D_{\theta}$ is a $\theta$-th order temporal difference operator, $L_{s}$ and $L_{r}$ are the spatial and temporal graph Laplacian matrices, respectively, $R(X)$ is the regularization term applied to $X$ describing its characteristics, and $\alpha\geq 0,\leavevmode\nobreak\ \beta\geq 0,\leavevmode\nobreak\ \gamma\geq 0$ are three parameters. Two common choices of $\theta$ are (1) $\theta=0$ that corresponds to $D_{\theta}=I$ and (2) $\theta=1$ used in TGSR. Additionally, $L_{s}$ can be a transformed version of the classical graph Laplacian $L$, e.g., $\tilde{L}=(L+\epsilon I)^{r}$ in the Sobolev method giraldo2022reconstruction , where $\epsilon>0$ and $r\geq 1$, which can be non-integer. The temporal graph Laplacian can be constructed by using the $\tau$-th order temporal difference operator, i.e., $L_{t}=D_{\tau}D_{\tau}^{T}$, for which case the temporal Laplacian can be expressed via the Frobenious norm $\operatorname*{tr}(XD_{\tau}D_{\tau}^{T}X^{T})=\left\lVert XD_{\tau}\right\rVert_{F}^{2}$ (see Tikhonov with $\tau=1$ and LRGTS with $\tau=2$). The regularization $R(X)$ can be chosen as the nuclear norm of $X$ if the underlying time-varying graph signal $X$ is low rank. Various models utilize different choices of $D_{\theta},L_{s}/\tilde{L},L_{t}$, and the regularization $R$. Leveraging the recent growth in deep learning, some time-varying graph signal recovery methods include unrolling technique kojima2023restoration , graph neural network (GNN) castro2023time , and joint sampling and reconstruction of time-varying graph signals xiao2023joint . In this work, we are dedicated to developing unsupervised time-varying graph signal recovery algorithms that do not involve or rely on data training.

Following the general framework (1), we propose a novel low-rank regularization $R(X)$ based on the error function (ERF) guo2021novel for sparse signal recovery (see Section 2.3). In addition, to handle non-Gaussian type of noise such as Laplace noise, we propose a variant model in which the Frobeinus norm based data fidelity term is replaced with the $\ell_{1}$-norm data fidelity (see Section 2.4). In Table 1, we provide a summary of the proposed models and relevant works pertaining to the general framework outlined in (1).

The major contributions of this work are described as follows.

1. 1.

   We develop a generalized time-varying graph signal recovery framework encompassing several state-of-the-art works as specific cases. We also develop two new models with a new regularization based on ERF.

2. 2.

   The proposed models combine high-order temporal smoothness and graph structures with the temporal correlation exploited by iteratively reweighted nuclear norm regularization.

3. 3.

   We propose an efficient algorithm for solving the proposed models. Convergence analysis has shown that the algorithm generates a sequence that converges to a stationary point of the problem.

4. 4.

   We conduct various numerical experiments, utilizing both synthetic and real-world datasets (specifically PM2.5 and sea surface temperature data), to validate the effectiveness of the proposed algorithm.

The subsequent sections of this paper are structured as follows. In Section 2, we introduce a pioneering framework for recovering time-varying graph signals, leveraging Sobolev smoothness and ERF regularization. Additionally, we put forth an efficient algorithm based on the alternating direction method of multipliers (ADMM) and iterative reweighting scheme. A comprehensive convergence analysis of the proposed algorithm is also provided. In Section 3, we present numerical experiments conducted on synthetic and real-world datasets sourced from environmental and epidemic contexts. Finally, Section 4 encapsulates our conclusions and outlines potential avenues for future research.

To enhance the low-rankness of a matrix, weighted nuclear norm minimization (WNNM) has been developed with promising performance in image denoising gu2014weighted . Specifically, the weighted nuclear norm (WNN) is defined as

where $\sigma_{i}(L)$ is the $i$-th singular value of $L$ in the decreasing order and the weight vector $\mathbf{w}=(w_{i})$ is in the non-decreasing order with $w_{i}\geq 0$ being the $i$-th weight. Choosing the weights is challenging in sparse and low-rank signal recovery problems. Iteratively reweighted L1 (IRL1) candes2008enhancing was proposed for the sparse recovery problem, where the weight is updated based on the previous estimate. It can solve many problems with complicated sparse regularizations, exhibiting improved sparsity and convergence speed.

In this work, we introduce a novel ERF-weighted nuclear norm based on the ERF regularizer guo2021novel and use linearization to obtain WNN. For any real matrix $X$ with $n$ singular values $\sigma_{1}(X)\geq\ldots\geq\sigma_{n}(X)$, the ERF-weighted nuclear norm is

where $\sigma$ serves as a filtering parameter. To solve the ERF-nuclear norm regularized minimization problem, we use iterative reweighting (linearization) to get WNN with adaptive weights.

Inspired by the Grünwald-Letnikov fractional derivative podlubny1999fractional , we introduce the total $\theta$-th order temporal forward difference matrix with a zero boundary condition, as shown below

Here the coefficients $\{C(i)\}_{i=0}^{k}$ are defined as

where $\Gamma(x)$ is the Gamma function. Notice that if $\theta$ is a positive integer, $k$ can be deterministic. For example, if $\theta=1$, then $k=1$ and we have $C(0)=1$ and $C(1)=-1$, which is reduced to the first-order finite difference case. If $\theta=2$, then it reduces to the temporal Laplacian operator. Generally if $\theta=n$, then only the first $n+1$ coefficients $\{C(i)\}_{i=0}^{n}$ are nonzero and thereby $k=n+1$. For any fractional value $\theta$, we have to choose the parameter $k.$ The difference matrix (4) is built upon the zero boundary condition, while other types of boundary conditions, e.g., Newmann and periodic boundary conditions, can also be used. Alternatively, we can use low-order difference schemes for boundary conditions, e.g., the first-order forward difference based on the first $m-1$ time points and the zeroth order for the last time point.

We propose the following ERF regularized time-varying graph signal recovery model

Here we use the least squares as a data fidelity term, the Sobolev smoothness of time-varying graph signals giraldo2022reconstruction as the graph regularization, and an ERF-based regularization defined in (3) for temporal low-rank correlation.

We apply ADMM with linearization to solve the problem (5). First, we introduce an auxiliary variable $Z$ to rewrite the problem (5) into an equivalent constrained problem:

Since the proximal operator of $\|\cdot\|_{\text{erf}}$ is difficult to compute, we apply linearization on the ERF term to obtain a WNN when solving the subproblem for $Z$. The ADMM iterates as follows,

where $\rho>0$ is a stepsize that affects the convergence; please refer to Theorem 2.1 for more details. We derive closed-form solutions for both $Z$- and $X$-subproblems in (6). Specifically for the $Z$-subproblem, it can be updated via the singular value thresholding operator, i.e.,

where $U\Sigma V^{T}$ is the singular value decomposition of $X+\widehat{Z}$, and $\operatorname{diag}(\cdot)$ is a diagonalization operator turning a vector into a diagonal matrix with the same entries as the vector. Here the shrink operator $\operatorname{shrink}(x,\xi)=\operatorname*{sign}(x)*\max(|x|-\xi)$ is implemented entrywise, where $\operatorname*{sign}(x)$ returns the sign of $x$ if $x\neq 0$ and zero otherwise.

In the $X$-subproblem, we can rewrite the second term of the objective function as

where $\otimes$ is the Kronecker product. Thus, the $X$-subproblem has the closed-form solution as

where $\widehat{J}=\operatorname{diag}(\operatorname*{vec}(J))$. Note that $\widehat{J}^{T}Y=Y$ since $\widehat{J}$ is a diagonal matrix with binary entries in the diagonal, whose nonzero entries correspond to the sampled spatial points. Furthermore, considering that the matrix $\widehat{J}+\alpha A^{T}A+\rho I$ is symmetric and positive definite, we perform its Cholesky factorization as $\widehat{J}+\alpha A^{T}A+\rho I=\tilde{L}\tilde{L}^{T}$. Subsequently, we leverage forward/backward substitution as a substitute for matrix inversion, thereby reducing computational time. The pseudo-code of the proposed approach for minimizing the model (5) is given in Algorithm 1.

In real-world applications, the type of noise could be unknown, and it is possible to encounter a mixture of different types of noise. To enhance the robustness against noise, we propose the second model,

Compared with (5), this new model utilizes the $\ell_{1}$-norm data fidelity to accommodate various types of noise. Because of the $\ell_{1}$ term, we need to introduce an additional variable $V$ to make the subproblems easy to solve. The equivalent constrained problem is

Therefore, the ADMM with linearization on the ERF term has the following subproblems

For the $V$-subproblem, we get the closed-form solution expressed via the shrinkage operator

Similar to Algorithm 1, the solution of the $Z$-subproblem is given by (7) with $\rho$ replaced by $\rho_{2}$. For the $X$-subproblem, we get the closed-form solution

The entire algorithm is described in Algorithm 2.

For simplicity, we define

and hence the augmented Lagrangian function is given by

The function $f$ is convex and continuously differentiable. In addition, $\nabla f$ is Lipschitz continuous with a constant $L$.

Following the work of qiu2017time , we generate $N=100$ nodes randomly from the uniform distribution in a $100\times 100$ square area. The graph weight is determined using $k$-nearest neighbors. Specifically, the weight between any two nodes is inversely proportional to the square of their Euclidean distance. We consider $k=5$ and visualize the corresponding graph in Fig. 1.

Denote the weight matrix by $W,$ its degree matrix $D,$ and the graph Laplacian $L$ has eigen-decomposition $L=U\Lambda U^{T},$ where $\Lambda=\text{diag}(0,\lambda_{2},\cdots,\lambda_{N})$. We further define $L^{-1/2}=U\Lambda^{-1/2}U^{T}$ where $\Lambda^{-1/2}=\text{diag}(0,\lambda_{2}^{-1/2},\cdots,\lambda_{N}^{-1/2})$. Starting from $x_{1},$ we generate the time-varying graph signal

where $f_{t}$ is an i.i.d. Gaussian signal rescaled to $\|f_{t}\|_{2}=\kappa$ and $\kappa$ corresponds to a temporal smoothness of the signal. Stacking $\{x_{t}\}$ as a column vector, we obtain a data matrix $X=[x_{1},x_{2},\cdots,x_{T}]$. We generate a low-rank data matrix obtained by starting with an empty matrix $X$ and repeating $X=[X,x_{1},\cdots,x_{10},x_{10},x_{9},\cdots,x_{1}]$ 10 times, thus also getting a $100\times 200$ data matrix. The measurement noise at each node is i.i.d. Gaussian noise $\mathcal{N}(0,\eta^{2}),$ where $\eta$ is the standard deviation.

Parameter tuning. For the proposed Algorithm 1, we fix the following parameters: $k=3$ and $\theta=1.8$ in the definition of fractional-order derivative (4); $\sigma=10^{3}$ in the definition of the ERF regularization (3); $\epsilon=0.1$ and $r=3$ in the Sobolev graph Laplacian; and the step size $\rho=10^{-6}$ in the ADMM iterations (6). In each set of experiments, we carefully tune two parameters $(\alpha,\beta)$ that determine the weights for the spatial-temporal smoothness and the low-rankness, respectively, in the proposed model (5). We choose the best combination of $(\alpha,\beta)$ among $\alpha\in\{0,10^{-5},10^{-4},10^{-3},10^{-2},10^{-1},1,10\}$ and $\beta\in\{0,10^{-8},10^{-7},10^{-6},10^{-5},10^{-4},10^{-3},10^{-2},10^{-1},1,10\}$. As demonstrated in Table 1, some competing methods are special cases of the proposed models, and hence, we only tune the parameters $\alpha,\beta$ for these methods while keeping other parameters fixed.

Reconstruction errors with respect to sampling rates. We begin by evaluating the performance of competing methods under different sampling rates. The smoothness level is set as $\kappa=1,$ while the standard deviation of the Gaussian noise is $\eta=0.1.$ The reconstruction performance is evaluated via RMSE, defined in (18), showing that the recovery errors of all the methods decrease with the increase of the sampling rates. The comparison results are visualized in Fig. 2. The proposed method achieves significant improvements over the competing methods. Surprisingly, LRDS, equipped with the nuclear norm, does not yield stable reconstruction performance in the low-rank case.

Reconstruction errors with respect to noise levels. We then investigate the recovery performance under different noise levels by setting the noise variance $\eta^{2}=\{0.01,0.1,0.2,0.4,0.6,0.8,1\}.$ In this set of experiments, we fix the sampling rate as 40% and smoothing level $\kappa=1.$ The noise level affects the magnitude of the least-squares fit, and as a result, we adjust the search window of $\alpha\in\{0,10^{-3},10^{-2},10^{-1},1,10,10^{2},10^{3},10^{4}$. The parameter $\beta$ remains the same: $\beta\in\{0,10^{-8},10^{-7},10^{-6},10^{-5},10^{-4},10^{-3},10^{-2},10^{-1},1,10\}$. The results are presented in Fig. 3, demonstrating the superior performance of the proposed Algorithm 1 under various noise levels.

In the real data experiments, we first test the daily mean Particulate Matter (PM) 2.5 concentration dataset from California provided by the US Environmental Protection Agency https://www.epa.gov/outdoor-air-quality-data. We used the data captured daily from 93 sensors in California for the first 200 days in 2015. The constructed graph is depicted in Fig. 4. In Fig. 5, we compare the average recovery accuracy of all the comparing methods over 50 trials when the sampling rates are $0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45$. In Table 2, we also compare the performance of Algorithm 1 and Algorithm 2, which shows Algorithm 2 can improve the accuracy of Algorithm 1 under some sampling rates with longer time in general.

Next, we test the sea surface temperature dataset, which was captured monthly by the NOAA Physical Sciences Laboratory (PSL). The data set can be downloaded from the PSL website https://psl.noaa.gov/. We use a subset of 200 time points on the Pacific Ocean within 400 months. The constructed graph is illustrated in Fig. 6. We see from Fig. 7 that the proposed algorithm outperforms other methods significantly and consistently across all sampling rates. In Table 3, we also compare the performance of Algorithm 1 and Algorithm 2, which indicates Algorithm 2 can improve the accuracy of Algorithm 1 under certain sampling rates but with more computational time in general.

Using the sea surface temperature data, we conduct an ablation study of the proposed model (5) without the smoothing regularization by setting $\alpha=0$ or without the low-rank ERF term by setting $\beta=0$. We plot the RMSE curves with respect to the sampling rates and the noise levels in Fig. 8, showing that the ERF regularization has a larger influence on the performance compared to the Sobolev-base graph Laplacian regularization.

Using the same sea surface temperature data, we investigate whether the proposed model (5) is sensitive to the parameters $(r,\epsilon)$ in defining the Sobolev-graph Laplacian and $\sigma^{2}$ in defining the ERF regularization. Fig. 9 shows that the proposed approach is not sensitive to various degrees of smoothness controlled by $r$ and $\epsilon$. Although the ERF regularization plays an important role in the recovery performance, as illustrated in the ablation study, the proposed model is not sensitive to the choice $\sigma^{2}$ as long as it is larger than 10,000.

In addition, we compare the proposed Algorithm 1 and Algorithm 2 using the sea surface temperature data and show the results in Tables 2 and 3. One can see that the two algorithms lead to similar RMSE, but Algorithm 2 is slower overall. We therefore prefer to use Algorithm 1 unless the data is heavily polluted by the non-Gaussian type of noise, such as Laplace noise.