Subset Random Sampling and Reconstruction of Finite Time-Vertex Graph Signals

TL;DR

This paper introduces a subset random sampling method for finite time-vertex graph signals with unknown spectral support, along with a reconstruction framework that leverages low-rank, sparsity, and smoothness priors, validated through experiments.

Contribution

It proposes a practical subset sampling scheme and a novel reconstruction framework for FTVGS with unknown spectral support, addressing real-world sampling constraints.

Findings

High-probability reconstruction guarantees under certain conditions

Effective reconstruction demonstrated through experiments

Framework leveraging low-rank, sparsity, and smoothness priors

Abstract

Finite time-vertex graph signals (FTVGS) provide an efficient representation for capturing spatio-temporal correlations across multiple data sources on irregular structures. Although sampling and reconstruction of FTVGS with known spectral support have been extensively studied, the case of unknown spectral support requires further investigation. Existing random sampling methods may extract samples from any vertex at any time, but such strategies are not friendly in practice, where sampling is typically limited to a subset of vertices and moments. To address this requirement, we propose a subset random sampling scheme for FTVGS. Specifically, we first randomly select a subset of rows and columns to form a submatrix, followed by random sampling within that submatrix. In theory, we provide sufficient conditions for reconstructing the original FTVGS with high probability. Additionally, we introduce a reconstruction framework incorporating low-rank, sparsity, and smoothness priors (LSSP), and verify the feasibility of the reconstruction and the effectiveness of the framework through experiments.

Tables

Table 1. TABLE I: Differences between three schemes for random sampling of FTVGS.

	Random sampling of MC	CCS	Subset random sampling
Assumption	rank-$r$, $\{{\mu }_1,{\mu }_2\}$-incoherence		rank-$r$, smooth, $\{{\mu }_1,{\mu }_2\}$-incoherence
Sufficient sample complexity	$O\left(r^2{n}_2\log \left(n_2\right)\right)$ [42] ($n_2=\max \{N,T\}$)	$O\left(r^2{n}_2{\log }^2\left(n_2\right)\right)$ [32]	$O\left(r^{2}N{\log }^2\left(n\right)\right)$

Table 2. TABLE II: Existing FTVGS reconstruction algorithms.

	$f_r\left(\widehat{𝐗}\right)$	$f_c({𝐅}_{𝒢},{𝐅}_{𝒯})$	$f_s\left(\widehat{𝐗}\right)$
SVT [34]	nuclear norm ${\Vert \widehat{𝐗}\Vert }_{\ast }={\sum }_i{\sigma }_i\left(\widehat{𝐗}\right)$	-	-
NC-MC [36]	non-convex regularizer $g\left({\lambda }_i\left({\widehat{𝐗}}^T\widehat{𝐗}\right)\right)$, $g\left(x\right)=\log \left(x^{1/2}+1\right)$	-	-
ReLaSP [46]	nuclear norm ${\Vert \widehat{𝐗}\Vert }_{\ast }$	${\Vert {𝐅}_{𝒢}\Vert }_1={\max }_j{\sum }_i^N\left|{𝐅}_{𝒢}(i,j)\right|$ or ${\Vert {𝐅}_{𝒯}\Vert }_1$	-
LIMC [47]	nuclear norm ${\Vert \widehat{𝐗}\Vert }_{\ast }$	-	${\Vert \widehat{𝐗}{𝐃}_2\Vert }_F^2+{\Vert {\widehat{𝐗}}^T{𝐃}_2^{\prime }\Vert }_F^2$ ${𝐃}_2^{\prime }$ is the second-order graph differential operator
LRDS [29]	nuclear norm ${\Vert \widehat{𝐗}\Vert }_{\ast }$	-	$\text{tr}\left({\left(\widehat{𝐗}{𝐃}_1\right)}^T𝐋\widehat{𝐗}{𝐃}_1\right)$
LRGTS [30]	nuclear norm ${\Vert \widehat{𝐗}\Vert }_{\ast }$	-	${\Vert \widehat{𝐗}{𝐃}_2\Vert }_F^2+\text{tr}\left({\widehat{𝐗}}^T𝐋\widehat{𝐗}\right)$

Table 3. TABLE III: Comparison of NRMSE for reconstruction algorithms on synthetic data.

${\alpha }_{RC}$	$90\%$	$80\%$	$70\%$	$60\%$
${\alpha }_{{𝐗}_{RC}}$	$90\%$	$80\%$	$70\%$	$60\%$
${\alpha }_{total}$	$72.65\%$	$50.80\%$	$34.61\%$	$21.71\%$
ICURC [32]	0.4378	0.6009	0.7149	0.8022
MC-NC [36]	0.4545	0.5922	0.7048	0.7976
LRDS [29]	0.3651	0.4889	0.5750	0.6555
LIMC [47]	0.1529	0.2146	0.2852	0.3049
ReLaSP [46]	0.1306	0.1805	0.2269	0.2506
LSSP	0.1118	0.1615	0.1939	0.2330

Table 4. TABLE IV: Comparison of NRMSE for reconstruction algorithms on METR-LA.

${\alpha }_{RC}$	$90\%$	$80\%$	$70\%$	$60\%$
${\alpha }_{{𝐗}_{RC}}$	$90\%$	$80\%$	$70\%$	$60\%$
${\alpha }_{total}$	$72.81\%$	$51.37\%$	$34.29\%$	$21.55\%$
ICURC [32]	0.4337	0.5966	0.7110	0.8075
MC-NC [36]	0.4375	0.5982	0.7207	0.8063
LRDS [29]	0.3931	0.5025	0.5917	0.6669
LIMC [47]	0.1409	0.2218	0.2739	0.3491
ReLaSP [46]	0.1311	0.2105	0.2296	0.2731
LSSP	0.1248	0.1934	0.2121	0.2457

Equations

$$\mathbf{Q}(i,m)=\left\{\begin{array}[]{ll}1,&{\text{if $i\in\mathcal{V}_{\mathcal{G}}$ is the initial vertex of edge $e_{m}$};}\\ -1,&{\text{if $i$ is the terminal vertex of edge $e_{m}$};}\\ 0,&{\text{otherwise}.}\end{array}\right.$$

$$\mathbf{X}(:,j{)}^T\mathbf{LX}(:,j)=\frac{1}{2}\sum _{i\in {\mathcal{V}}_{\mathcal{G}}}\sum _{l\sim i}(\mathbf{X}(l,j)-\mathbf{X}(i,j){)}^2,$$

$${\mathbf{D}}_1=\left[\begin{array}{cccc}-1& & & \\ 1& -1& & \\ & 1& \ddots & \\ & & \ddots & -1\\ & & & 1\end{array}\right],{\mathbf{D}}_2=\left[\begin{array}{cccc}1& & & \\ -2& 1& & \\ 1& -2& & \\ & 1& \ddots & \\ & & \ddots & -2\\ & & \ddots & 1\end{array}\right],$$

$$(\mathbf{X}{\mathbf{D}}_1)(i,j)=-\mathbf{X}(i,j)+\mathbf{X}(i,j+1),$$

$$(\mathbf{X}{\mathbf{D}}_2)(i,j)=\mathbf{X}(i,j)-2\mathbf{X}(i,j+1)+\mathbf{X}(i,j+2).$$

$$▽\mathbf{X}(i,j)=\left[\begin{array}{c}({\mathbf{Q}}^T\mathbf{X})(i,j)\\ (\mathbf{X}{\mathbf{D}}_1)(i,j)\end{array}\right].$$

$$\max \left\{\left|\mathbf{X}(:,j{)}^T\mathbf{LX}(:,j)\right|,||▽\mathbf{X}(i,j)||,\left|(\mathbf{X}{\mathbf{D}}_2)(i,j)\right|\right\}\le C$$

$$\mathbf{X}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T,$$

$$\frac{N}{r}||\mathbf{U}|{|}_{2,\infty }^2\le {\mu }_1(\mathbf{X}),\frac{T}{r}||\mathbf{V}|{|}_{2,\infty }^2\le {\mu }_2(\mathbf{X}),$$

$$\mathcal{I}\subseteq {\mathcal{V}}_{\mathcal{G}},\mathcal{J}\subseteq {\mathcal{V}}_{\mathcal{T}}.$$

$$\mathcal{S}=\{(i,j):i\in \mathcal{I},j\in \mathcal{J}\}.$$

$$|\mathcal{I}|\ge 3r{\mu }_1(\mathbf{X})\frac{\ln (2r/\delta )}{{ϵ}^2},|\mathcal{J}|\ge 3r{\mu }_2(\mathbf{X})\frac{\ln (2r/\delta )}{{ϵ}^2}$$

$${\mathbf{X}}_R=\mathbf{U}(\mathcal{I},:)\Sigma {\mathbf{V}}^T={\mathbf{U}}_R{\Sigma }_R{\tilde{\mathbf{V}}}_R^T{\mathbf{V}}^T={\mathbf{U}}_R{\Sigma }_R{\mathbf{V}}_R^T,$$

$${\mathbf{X}}_{RC}$$

$$={\mathbf{U}}_{RC}{\Sigma }_{RC}{\mathbf{V}}_{RC}^T,$$

$$||{\mathbf{U}}_{RC}|{|}_{2,\infty }\le \kappa (\mathbf{X})\sqrt{\frac{{\mu }_1(\mathbf{X})r}{(1-\eta )|\mathcal{I}|}},$$

$$||{\mathbf{V}}_{RC}|{|}_{2,\infty }\le \frac{\kappa (\mathbf{X})}{1-\eta }\sqrt{\frac{{\mu }_2(\mathbf{X})Nr}{|\mathcal{I}||\mathcal{J}|}},$$

$$|\mathcal{S}|\ge \frac{32\beta (\kappa (\mathbf{X}){)}^4{r}^{2}N}{(1-\eta {)}^3}{\mu }_1(\mathbf{X}){\mu }_2(\mathbf{X})\frac{|\mathcal{I}|+|\mathcal{J}|}{|\mathcal{I}|}{\log }^2(2n),$$

$$(1-\delta {)}^2(1-p{)}^2-\frac{6\log n}{(|\mathcal{I}|+|\mathcal{J}|{)}^{2\beta -2}}-n^{2-2{\beta }^{1/2}}.$$

$$\underset{\hat{\mathbf{X}}}{\min }{f}_{\text{obj}}=f_r(\hat{\mathbf{X}})+{\gamma }_c{f}_c({\mathbf{F}}_{\mathcal{G}},{\mathbf{F}}_{\mathcal{T}})+{\gamma }_s{f}_s(\hat{\mathbf{X}})$$

$$\displaystyle\begin{array}[]{c@{\quad\quad\quad}l@{\quad}l}s.t.&\hat{\mathbf{X}}=\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{G}},\\ &\hat{\mathbf{X}}^{T}=\Psi_{\mathcal{T}}\mathbf{F}_{\mathcal{T}},\\ &\hat{\mathbf{X}}(i,j)=\mathbf{X}(i,j),(i,j)\in\mathcal{S},\end{array}$$

$$f_{\text{obj}}$$

$$+{\gamma }_{\mathcal{T}}||{\mathbf{W}}_{\mathcal{T}}\odot {\mathbf{F}}_{\mathcal{T}}|{|}_1+{\gamma }_d||\hat{\mathbf{X}}{\mathbf{D}}_2|{|}_F^2,$$

$$|\mathcal{I}|=\frac{\mathrm{Round}({\alpha }_{RC}N)}{N},|\mathcal{J}|=\frac{\mathrm{Round}({\alpha }_{RC}T)}{T},$$

$$\mathrm{NRMSE}({\mathbf{X}}_a,{\mathbf{X}}_b)=\frac{||{\mathbf{X}}_a-{\mathbf{X}}_b|{|}_F}{||{\mathbf{X}}_a|{|}_F}.$$

$$|\mathcal{I}|\ge 3N{\xi }_1\frac{\ln (2r/\delta )}{{ϵ}^2}$$

$$|\mathcal{I}|\ge 3r{\mu }_1(\mathbf{X})\frac{\ln (2r/\delta )}{{ϵ}^2}.$$

$${\xi }_2=||{\mathbf{V}}_R|{|}_{2,\infty }^2=||\mathbf{V}{\tilde{\mathbf{V}}}_R|{|}_{2,\infty }^2\stackrel{(a)}{=}||\mathbf{V}|{|}_{2,\infty }^2\le \frac{{\mu }_2(\mathbf{X})r}{T},$$

$$|\mathcal{J}|\ge 3r{\mu }_2(\mathbf{X})\frac{\ln (2r/\delta )}{{ϵ}^2}.$$

$$||{\mathbf{U}}_R|{|}_{2,\infty }$$

$$=\kappa (\mathbf{X})\frac{1}{{\sigma }_{\min }(\mathbf{U}(\mathcal{I},:))}\sqrt{\frac{{\mu }_1(\mathbf{X})r}{N}}.$$

$$\sqrt{\frac{(1-\eta )|\mathcal{I}|}{N}}\le {\sigma }_{\min }(\mathbf{U}(\mathcal{I},:))$$

Full text

I Introduction

Real-world data often exhibits a non-regular structure and is time-varying, as seen in sensor networks [1, 2], brain networks [3, 4], and social networks [5, 6]. Graph Signal Processing (GSP) has become widespread, extending classical signal processing techniques rooted in the Fourier Transform to include the Graph Fourier Transform (GFT) [7], graph fractional Fourier transform [8], graph signal sampling [9], and graph aggregated sampling [10]. Time-vertex graph signals (TVGS) naturally extend graph signals by associating each graph vertex with a time-series signal [11]. This structure enables data processing from both vertex and time perspectives. Consequently, TVGS is well-suited for representing and analyzing irregular time-varying data, leading to the development of TVGS processing [12, 13, 14], including graph filtering [15, 16], data inference [17, 18], and spatiotemporal graph learning [19, 20]. Applications in domains such as computer vision [21], anomalous sound detection [22], and oceanography [23] stand to benefit significantly from advances in TVGS processing.

Efficient sampling and reconstruction of TVGS are crucial for accurately capturing and analyzing real-world phenomena. Since the topology of TVGS handles data from both the vertex and temporal perspectives simultaneously, the volume of data grows exponentially [13, 14]. This presents significant challenges in data acquisition, storage, and processing. Therefore, sampling theory is a critical aspect of time-point graph signal processing, aiming to recover the entire original signal from partial observations. This paper focuses on the sampling and reconstruction problem of finite time-vertex graph signals (FTVGS). In FTVGS, each vertex is associated with a finite-length time sequence, allowing the signal to be modeled as a matrix [24].

Discrete signal processing (DSP) and GSP have been combined to introduce the joint time-vertex Fourier transform (JFT), which provides the spectral expansion of FTVGS [12]. Based on the known JFT spectrum, some progress has been made in sampling and reconstructing FTVGS [25, 24, 26, 27]. We call these methods deterministic sampling. However, in practice, we often cannot obtain the spectral information of the signal before sampling. For example, in sensor networks, each row of a matrix represents the time series collected by a sensor, while each column represents a traditional graph signal. Before deploying the sensors, the JFT spectrum of the FTVGS is unknown, making sampling methods that rely on known JFT spectra impractical.

In this case, we consider using random sampling to address the challenge of unknown JFT spectra. One approach is to sample the entire FTVGS through an unbiased stochastic process, followed by formulating an optimization problem for reconstruction based on assumptions such as low-rank or smoothness [28, 29, 30]. However, this method assumes that all elements of the FTVGS are sampled with the same probability, which is often unrealistic due to factors such as participant willingness or environmental conditions. Lorenzo et al. designed a probabilistic sampling strategy where each vertex at every time instant is sampled with a given probability [31]. Cai et al. proposed the cross-concentrated sampling (CCS) strategy to ensure certain elements remain unobserved [32]. However, these sampling schemes do not reduce the number of nodes and time instants generating samples.

In practice, we prefer to sample only a subset of rows and columns, while others produce no samples at all due to hardware, environmental, or cost limitations. We propose a subset random sampling scheme for FTVGS, significantly reducing the number of rows and columns producing samples. That is, we randomly select a subset of rows and columns from the matrix to form a submatrix and then randomly sample within the submatrix.

The subset random sampling scheme is straightforward but results in samples with entire rows and columns missing, which we call structural missing. In [33], we first proposed the subset random sampling scheme and established theoretical guarantees for reconstructing the original FTVGS under structural missing. This paper provides more detailed theoretical analyses and proofs regarding the sufficient conditions for guaranteed reconstruction with high probability.

However, most existing algorithms for FTVGS reconstruction from randomly observed samples rely on the uniform sampling assumption. Some studies reconstruct the original signal by relaxing the rank minimization problem [34, 35, 36]. Qiu et al. extended the definition of smooth signals from static graph to FTVGS, utilizing the smoothness of temporal differences for FTVGS reconstruction [28]. To accelerate algorithm convergence, Giraldo et al. proposed a reconstruction method based on the Sobolev smoothness of FTVGS [37]. The analysis and comparison of reconstruction algorithms are detailed in Table I.

Existing algorithms are rarely designed for structural missing data. To address the reconstruction of signals from samples restricted to the submatrix, we introduce a reconstruction framework with low-rank, sparse, and smooth priors (LSSP), which was not discussed in our previous work [33]. To our knowledge, no existing reconstruction algorithm has leveraged all three signal priors simultaneously. The algorithm within the LSSP framework efficiently recovers the original FTVGS even when rows and columns are missing.

In brief, the contributions of this paper are summarized as follows:

• •

We propose a more practical subset random sampling scheme for FTVGS. In this scheme, samples are only observed from the submatrix of FTVGS, effectively reducing the number of sampled vertices (rows) and instants (columns).

• •

We prove sufficient conditions for reconstructing the original FTVGS from samples obtained through subset random sampling schemes, including lower bounds on the number of selected rows and columns (Lemma 1), as well as a lower bound on the number of observed samples (Theorem 1).

• •

We introduce an LSSP reconstruction framework that deeply exploits the prior information of the FTVGS to address the challenges posed by structural missing. The feasibility of the reconstruction and the effectiveness of the framework are verified through experiments.

The rest of the paper is organized as follows. In Section II, we introduce the theoretical and algorithmic studies on random sampling and reconstruction of FTVGS. In Section III we introduce the FTVGS model. In Section IV, we present the subset random sampling scheme for FTVGS and discuss the conditions for reconstructing the original FTVGS. In Section V, we design an LSSP reconstruction framework tailored to the samples with structural missing. Numerical results are provided in Section VI, and conclusions are drawn in Section VII.

The notation in this paper is as follows. Flowery letters are used to denote sets. For a matrix $\mathbf{X}\in\mathbb{C}^{N\times T}$, let $\mathbf{X}(i,j)$ denotes the $(i,j)$-th entry of $\mathbf{X}$, and $\mathbf{X}(\mathcal{I},:)$, $\mathbf{X}(:,\mathcal{J})$, and $\mathbf{X}(\mathcal{I},\mathcal{J})$ denote the row submatrix with row indices $\mathcal{I}$, the column submatrix with column indices $\mathcal{J}$, and the submatrix of $\mathbf{X}$ composed of rows indexed by $\mathcal{I}$ and columns indexed by $\mathcal{J}$, respectively. $\sigma_{i}(\mathbf{X})$ represents the $i$-th largest singular value of $\mathbf{X}$, while $\lambda_{i}(\mathbf{X})$ is the $i$-th smallest eigenvalue of $\mathbf{X}$. $||\mathbf{X}||_{F}=\left(\sum_{i}\sum_{j}|\mathbf{X}(i,j)|^{2}\right)^{1/2}$ denotes the Frobenius norm.

II Related Works

Sampling methods that utilize prior information about the signal to determine specific sample indices are known as deterministic sampling. In contrast, methods that do not predefine the sampling locations are referred to as random sampling. This section provides a brief review of research related to both sampling approaches, with the main focus of this paper being the exploration of random sampling theory for FTVGS.

II-A Deterministic sampling of FTVGS

Based on the known JFT spectral support, the separate sampling scheme applies graph sampling in the vertex and time domains of FTVGS, respectively, resulting in a submatrix containing a part of rows and columns as samples [25]. However, there remains the potential to reduce the number of samples further. Sheng et al. jointly processed the FTVGS spectra in both vertex and time domains, established the lower bound of the total sampling density based on the JFT spectrum, and proposed a feasible multi-band sampling scheme [24].

In addition, Wei et al. designed a reconstruction scheme for FTVGS using subspace prior knowledge and proposed a generalized sampling framework for FTVGS by applying the sampling framework of shift-invariant spaces [27]. Xiao et al. proposed a sampling method for FTVGS based on known graph spectrum support and introduced a robust reconstruction algorithm leveraging the smoothness of the signal [26].

The above-mentioned studies are sampling and reconstruction of FTVGS based on the spectrum of signals. However, in practice, the spectral information is unknown prior to sampling. Therefore, we consider random sampling of FTVGS when the JFT spectrum is unavailable.

II-B Random Sampling of FTVGS

Random sampling of FTVGS is similar to the problem of random sampling of matrices. The pioneering work [38] explored the matrix completion (MC) problem based on the Bernoulli sampling model, where the sample index set is obtained through Bernoulli sampling. Moreover, [38] and [39] proved a lower bound on the number of samples required for accurate recovery with high probability.

Subsequent research on MC expanded beyond Bernoulli sampling to uniform sampling with and without replacement [39, 40]. In addition, [41] demonstrated the equivalence of Bernoulli sampling and uniform sampling in MC. Moreover, the lower bound on the number of samples required for signal reconstruction is further reduced in [42].

Numerous studies on sampling and reconstruction of FTVGS have adopted an approach based on sampling via an unbiased stochastic process, leading to reconstruction algorithms that leverage priors such as smoothness [28, 29, 30]. Although these studies provide insight into the sampling and reconstruction of FTVGS, they assume that all elements of the matrix are sampled with the same probability (Fig. 1 (a)), which is not always feasible in practice. For example, in sensor networks, practitioners typically aim to deploy sensors at a subset of candidate locations while minimizing operational durations.

A subsequent study combined uniform sampling with CUR decomposition, proposing CCS [32]. Specifically, CCS randomly samples within row and column submatrices of FTVGS to obtain the samples (Fig. 1 (b)). While CCS ensures that some elements of the FTVGS remain unobserved, it still requires that every vertex and instant produce samples.

In contrast to existing methods, our proposed subset random sampling scheme restricts the sample indices to a subset of rows and columns of FTVGS, ensuring that no samples are produced from the unselected rows and columns (Fig. 2 ). This is motivated by practical scenarios where sampling is often desired only at certain vertices and instances. For example, in sensor networks, we aim to deploy as few sensors as possible. Besides, hardware limitations may prevent sampling at every instance.

To address structurally missing samples in subset random sampling, we introduce an LSSP reconstruction framework and validate its feasibility and effectiveness. The LSSP framework integrates low-rank property, smoothness, and the correlations between rows and columns of FTVGS.

III Model

An undirected graph can be represented as $\mathcal{G}=(\mathcal{V}_{\mathcal{G}},\mathcal{E}_{\mathcal{G}})$, where $\mathcal{V}_{\mathcal{G}}$ is the set of vertices with $|\mathcal{V}_{\mathcal{G}}|=N$, $\mathcal{E}_{\mathcal{G}}=\{e_{1},\dots,e_{M}\}$ is the set of edges. We consider every vertex in $\mathcal{V}_{\mathcal{G}}$ relates to a time sequence of length $T$, which can be represented by a directed cyclic graph $\mathcal{T}=(\mathcal{V}_{\mathcal{T}},\mathcal{E}_{\mathcal{T}})$. Such signals are FTVGS [13], whose graphical model is $\mathcal{G}\times\mathcal{T}$ (as shown in Fig. 2), where $\times$ is the Cartesian product. An FTVGS is given to the matrix $\mathbf{X}\in\mathbb{C}^{N\times T}$. Each row of $\mathbf{X}$ is a finite time sequence on a vertex, and each column is a static graph signal.

Suppose each edge of $\mathcal{G}$ is assigned an orientation, which is arbitrary but fixed. We can define the incidence matrix $\mathbf{Q}\in\mathbb{R}^{N\times M}$ as

[TABLE]

Then the gradient on $\mathcal{G}$ can be defined as $\mathbf{Q}^{T}\mathbf{X}$ and the graph Laplacian matrix is $\mathbf{L}_{\mathcal{G}}=\mathbf{Q}\mathbf{Q}^{T}$ [43]. Let $l\sim i$ indicate that there is an edge between vertices $l,i\in\mathcal{V}_{\mathcal{G}}$. A typical measure of how much a signal varies on $\mathcal{G}$ for every $j\in\mathcal{V}_{\mathcal{T}}$ is the quadratic

[TABLE]

which is the sum of the neighborhood variation of all vertices, i.e., graph total variation (TV).

Then we introduce first and second-order temporal differential operators

[TABLE]

where $\mathbf{D}_{1}\in\mathbb{R}^{T\times(T-1)}$ and $\mathbf{D}_{2}\in\mathbb{R}^{T\times(T-2)}$. For each $i\in\mathcal{V}_{\mathcal{G}}$, the first and second order differences of the signal on $\mathcal{T}$ can be represented as [44]

[TABLE]

and

[TABLE]

Subsequently, the gradient of $\mathbf{X}$ in the graph and time domains is

[TABLE]

Our objective is to propose a subset random sampling scheme for FTVGS and to prove the lower bound on the number of samples required for random sampling on the submatrix of FTVGS. So we begin with a formal expression of the widely used assumptions.

Assumption 1**.**

[38, 45]** (low-rank, smoothness) An FTVGS $\mathbf{X}\in\mathbb{C}^{N\times T}$ is a rank-$r$ matrix, where $r<\min\{N,T\}$. For all $i\in\mathcal{V}_{\mathcal{G}},j\in\mathcal{V}_{\mathcal{T}}$,

[TABLE]

for some constant $C$.

The low rank implies that $\mathbf{X}$ can be expressed as a matrix with only $r$ non-zero rows (or columns) through a row (or column) linear transformation. This is consistent with the assumption of bandlimited JFT spectra in the sampling and reconstruction theory of FTVGS, as the JFT is a type of linear transformation [24]. In particular, the JFT spectrum can be expressed as $\Psi_{\mathcal{G}}^{T}\mathbf{X}\Psi_{\mathcal{T}}$, where $\Psi_{\mathcal{G}}$ and $\Psi_{\mathcal{T}}$ are the GFT basis matrix and discrete Fourier transform (DFT) basis matrix, respectively. The JFT spectral bandlimitness then implies that each row and column of $\Psi_{\mathcal{G}}^{T}\mathbf{X}\Psi_{\mathcal{T}}$ is tending to be sparse.

The quadratic $\mathbf{X}(:,j)^{T}\mathbf{L}\mathbf{X}(:,j)$ reflects how the $\mathbf{X}(:,j)$ varies with respect to the underlying graph structure, while the time-domain difference $(\mathbf{X}\mathbf{D}_{2})(i,j)$ captures the changes in the $\mathbf{X}(i,:)$ over time. The gradient $\triangledown\mathbf{X}(i,j)$ indicates both the difference between an element and its neighboring vertices and the variation across time. The fact that quadratic on the graph $\mathcal{G}$, second-order difference on $\mathcal{T}$, and gradient are all bounded implies that $\mathbf{X}$ is smooth. That is, there are no abrupt changes in $\mathbf{X}$ in either graph or temporal domains. Correspondingly, in the context of TVGS processing theory, the smoothness of $\mathbf{X}$ somehow indicates that its JFT spectrum exhibits low-frequency characteristics.

Since $\mathbf{X}$ is rank-$r$, its thin singular value decomposition (SVD) is

[TABLE]

where $\mathbf{U}\in\mathbb{R}^{N\times r}$, $\bm{\Sigma}\in\mathbb{R}^{r\times r}$ and $\mathbf{V}\in\mathbb{R}^{T\times r}$. Then we have the incoherence assumption.

Assumption 2**.**

[38, 39]** (incoherence) An FTVGS $\mathbf{X}\in\mathbb{C}^{N\times T}$ satisfying 1 is $\{\mu_{1},\mu_{2}\}$-incoherence for some constants $1\leq\mu_{1}(\mathbf{X})\leq\frac{N}{r},1\leq\mu_{2}(\mathbf{X})\leq\frac{T}{r}$ provided

[TABLE]

where $||\mathbf{U}||_{2,\infty}:=\max_{i}\left(\sum_{j}(\mathbf{U}(i,j))^{2}\right)^{1/2}$ denotes the largest row-wise $\ell_{2}$-norm..

Coherence indicates whether the energy of $\mathbf{X}$ is uniformly distributed across columns and rows. Low coherence means that the structure of $\mathbf{X}$ is relatively uniform, making it easier to obtain uniform samples during random sampling.

IV Subset Random Sampling

IV-A Sampling Scheme

As mentioned in the sensor network example in Section I, not all sensors can reliably collect data at every moment. In practice, it is preferable to activate only a subset of sensors and sample at a portion of the moments. To address this need, we propose a subset random sampling scheme for signals, as illustrated in Fig. 2.

Formally, we first randomly select a subset of rows and columns from the FTVGS without replacement, i.e.,

[TABLE]

Let $\mathbf{X}(\mathcal{I},\mathcal{J})$ be the submatrix of FTVGS $\mathbf{X}$ composed of rows indexed by $\mathcal{I}$ and columns indexed by $\mathcal{J}$, the samples observed from $\mathbf{X}$ be indexed by a set $\mathcal{S}$, and $|\mathcal{S}|$ denote the number of samples. Subsequently, we uniformly observe samples within the submatrix $\mathbf{X}(\mathcal{I},\mathcal{J})$ with replacement (Fig. 2), i.e.,

[TABLE]

Denote $\mathbf{X}(i,j),(i,j)\in\mathcal{S}$ as $\mathbf{X}_{\mathcal{S}}$.

Unlike existing methods, we would like to sample randomly on the submatrix $\mathbf{X}(\mathcal{I},\mathcal{J})$ instead of $\mathbf{X}$, ensuring that no samples are observed on unselected rows and columns. We summarize the sampling process as Algorithm 1.

The subset random sampling scheme is straightforward but violates the common assumption in MC problems that each row and column must produce at least one sample. To the best of our knowledge, theoretical guarantees for reconstructing the original FTVGS when the samples are structurally missing in entire rows and columns have not been explored. Therefore, we then study subset random sampling from a theoretical perspective.

IV-B Sufficient Conditions

The conditions for FTVGS reconstruction consist of two parts: the lower bounds on the number of retained rows and columns, and the lower bound on the total number of samples observed. First, consider the lower bounds for $|\mathcal{I}|$ and $|\mathcal{J}|$. If too few rows or columns are retained, even when $|\mathcal{I}|<r$ or $|\mathcal{J}|<r$, it becomes challenging to ensure accurate FTVGS reconstruction, even if we observe the entire $\mathbf{X}(\mathcal{I},\mathcal{J})$.

The following lemma ensures that the rank of $\mathbf{X}(\mathcal{I},\mathcal{J})$ remains $r$ with high probability, which is crucial for the subsequent theoretical derivation. Lemma 1 also provides the lower bounds on the number of selected rows and columns. The proof of Lemma 1 is presented in APPENDIX A.

Lemma 1**.**

For an FTVGS $\mathbf{X}\in\mathbb{C}^{N\times T}$ satisfying 1 and 2, let $\mathbf{X}_{R}$ and $\mathbf{X}_{RC}$ be obtained by Algorithm 1. Then if

[TABLE]

for some $0<\delta<1$ and $0<\epsilon<1$, $\mathrm{rank}(\mathbf{X}_{R})=r$ holds with probability at least $1-\delta$ and $\mathrm{rank}(\mathbf{X}_{RC})=r$ holds with probability at least $(1-\delta)^{2}$.

With Lemma 1, we obtain the submatrix $\mathbf{X}_{RC}$. Following Algorithm 1, we will subsequently sample randomly on $\mathbf{X}_{RC}$. We wish to prove a lower bound on the number of samples required, which depends not only on the rank of $\mathbf{X}_{RC}$ but also on its incoherence properties. To this end, we prove Lemma 2, which demonstrates how the properties of the original signal $\mathbf{X}$ are converted to those of its submatrix $\mathbf{X}_{RC}$.

Let $\mathbf{X}_{R}=\mathbf{U}(\mathcal{I},:)\Sigma\mathbf{V}^{T}$. Then $\mathbf{U}(\mathcal{I},:)\Sigma$ has the thin SVD $\mathbf{U}_{R}\Sigma_{R}\tilde{\mathbf{V}}^{T}_{R}$, where $\tilde{\mathbf{V}}_{R}\in\mathbb{R}^{r\times r}$ is an orthonormal matrix and $\mathbf{U}_{R}\in\mathbb{R}^{r\times|\mathcal{I}|}$. The thin SVD of $\mathbf{X}_{R}$ is

[TABLE]

where $\mathbf{V}_{R}=\mathbf{V}\tilde{\mathbf{V}}_{R}$ is a matrix with orthogonal columns. Likewise, we have $\mathbf{X}_{RC}=\mathbf{U}_{R}\Sigma_{R}\mathbf{V}_{R}(\mathcal{J},:)^{T}$ and the thin SVD $\Sigma_{R}\mathbf{V}_{R}(\mathcal{J},:)^{T}=\tilde{\mathbf{U}}_{RC}\Sigma_{RC}\mathbf{V}_{RC}^{T}$, where $\tilde{\mathbf{U}}_{RC}\in\mathbb{R}^{r\times r}$ is an orthonormal matrix and $\mathbf{V}_{RC}\in\mathbb{R}^{r\times|\mathcal{J}|}$. The thin SVD of $\mathbf{X}_{RC}$ is

[TABLE]

where $\mathbf{U}_{RC}=\mathbf{U}_{R}\tilde{\mathbf{U}}_{RC}$ is a matrix with orthogonal columns.

Lemma 2**.**

Under the assumption of Lemma 1, if $\mathcal{I}$ and $\mathcal{J}$ satisfy the inequality of Lemma 1, then the following conditions hold with probability at least $(1-p)^{2},p=r\left[\frac{e^{-\eta}}{(1-\eta)^{1-\eta}}\right]^{\log r}$:

[TABLE]

[TABLE]

where $\kappa(\mathbf{X})$ is the condition number of $\mathbf{X}$ and $\eta\in[0,1)$.

The proof of Lemma 2 is presented in APPENDIX  B. Lemma 2 shows how the properties of the original signal $\mathbf{X}$ are converted to those of its submatrix $\mathbf{X}(\mathcal{I},\mathcal{J})$.

Having Lemma 1 and Lemma 2, we then introduce the second part of the conditions for FTVGS reconstruction: the lower bound on the number of samples $|\mathcal{S}|$. This forms the main theorem of this paper, following the theory of MC problems. But here we observe $\mathcal{S}$ from $\mathbf{X}_{RC}$, rather than from $\mathbf{X}$.

Theorem 1**.**

For an FTVGS $\mathbf{X}\in\mathbb{C}^{N\times T}$ satisfying 1 and 2, suppose that $\mathbf{X}_{RC}$ and $\mathcal{S}$ is derived from Algorithm 1. Then if

[TABLE]

where $\beta>1$, $\kappa(\mathbf{X})$ is the condition number of $\mathbf{X}$, $\eta\in[0,1)$ and $n=\max{\{|\mathcal{I}|,|\mathcal{J}|\}}$, the $\mathbf{X}$ can be recovered from $\mathcal{S}$ with probability at least

[TABLE]

The proof of Theorem 1 is presented in APPENDIX C.

Lemma 1 and Theorem 1 provide the conditions for reconstructing the original FTVGS from the subset random sampled samples. Among the probabilities in Theorem 1, $(1-\delta)^{2}$ is the probability that guarantees $\mathrm{rank}(\mathbf{X}_{RC})=\mathrm{rank}(\mathbf{X})=r$ and $(1-p)^{2}$ is the probability that guarantees the incoherence of $\mathbf{X}_{RC}$.

We summarize the differences between random sampling of MC, CCS, and subset random sampling in Table I. It can be observed that when $N$ and $T$ are sufficiently large, uniform random sampling requires a relatively smaller number of samples, while the sample complexities of CCS and subset random sampling are comparable. It is because uniform sampling imposes no constraints, assuming each vertex at each instant is sampled with equal probability, whereas CCS and subset random sampling restrict observable positions. *It is worth noting that subset random sampling does not have a higher sample complexity than CCS, although it constrains the vertices and instants that produce samples to a subset. * For smaller $N$ and $T$, the parameters of subset random sampling may increase the required number of samples, resulting in a higher lower bound. This is attributed to the fact that increased sampling constraints necessitate more samples to uniquely determine the original signal.

Our proposed subset random sampling scheme restricts the sample indices within $(\mathcal{I},\mathcal{J})$, meaning that no samples will be drawn from unselected rows and columns. However, this comes at the cost of some reconstruction error, as detailed in Section VI. Because there are rows and columns without samples, MC methods and CCS cannot complete the reconstruction. Methods based on TV inpainting and our proposed reconstruction framework can recover the FTVGS, but cannot completely eliminate the error.

V Reconstruction Framework

In Section IV, we established the conditions for reconstructing the original FTVGS from samples indexed by $\mathcal{S}$ through theoretical analysis. In this section, together with Section VI, we aim to verify the feasibility of FTVGS reconstruction from the perspective of the reconstruction algorithm.

MC has been a topic of extensive research, and well-established frameworks now exist. To address the challenge of reconstructing a complete signal from a limited number of randomly observed samples, an effective approach is to incorporate priors or regularization. A common prior assumes that the original matrix is low-rank or approximately low-rank. By modeling the MC process as a low-rank matrix approximation, it can then be solved using optimization methods.

Our Algorithm 1 meets the practical need to reduce the number of vertices and moments to be sampled, but it introduces the issue of structural missing, where entire rows and columns are absent (see Fig. 2). Reconstructing the FTVGS from these structurally incomplete samples goes beyond the scope of traditional MC conditions and algorithms. This is because the low-rank prior and its variants are insufficient in regularizing structural incompleteness, necessitating additional regularization to compensate for the limitations of low-rankness in the presence of structural missing.

To tackle the challenges posed by structural missing, we aim to more fully exploit the prior information of the FTVGS. Based on our assumptions about the FTVGS, we derive three types of regularization:

$f_{r}(\mathbf{X})$ stands for rank regularization of the FTVGS. Based on the general 1, we assume that the original signal exhibits high redundancy and can be represented in a more compact form, implying that FTVGS is either low-rank or approximately low-rank. To some degree, low-rank characterizes the global prior of the FTVGS, so the underlying matrix is regularized by the low-rank prior. 2. 2.

$f_{c}(\mathbf{F}_{\mathcal{G}},\mathbf{F}_{\mathcal{T}})$ denotes a priori regularization based on the spatial and temporal correlations of FTVGS. The low-rank property of FTVGS indicates strong correlations between its rows and columns, allowing $\mathbf{X}$ to be represented in a low-dimensional form through linear transformations. Based on TVGS processing theory, we apply the GFT and DFT bases to obtain a low-dimensional representation of $\mathbf{X}$. In the graph domain, we have $\mathbf{X}=\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{G}}$, where $\mathbf{F}_{\mathcal{G}}$ is the GFT spectrum, and each $\mathbf{F}_{\mathcal{G}}(:,j)$ is the bandlimited spectrum of the graph signal $\mathbf{X}(:,j)$. In the time domain, $\mathbf{X}^{T}=\Psi_{\mathcal{T}}\mathbf{F}_{\mathcal{T}}$, where each $\mathbf{F}_{\mathcal{T}}(:,i)$ represents the bandlimited spectrum of the time sequence $\mathbf{X}(i,:)$. To incorporate the sparse prior of FTVGS, we introduce sparse regularization of the graph frequency and the time frequency. 3. 3.

$f_{s}(\mathbf{X})$ represents smoothing regularization. In 1, we also introduced the smooth feature of the FTVGS. This implies that, in both the graph and time domains, the difference between an entry and its neighboring entries is typically small. Commonly studied regularization includes TV and its variants, as well as the smoothness of temporal difference signals.

Now, the method based on three kinds of regularization can be modeled as the following optimization problem

[TABLE]

where $\gamma_{c}$ and $\gamma_{s}$ are non-negative regularization parameters.

Some existing studies have not investigated the theoretical guarantees of Algorithm 1, but many reconstruction algorithms incorporating multiple priors have been proposed. In the LSSP reconstruction framework, we summarize some of these existing methods in Table II.

In Table II, the low-rank variant serves as a global regularization by constraining the low-rank structure of the underlying matrix [34, 36]. Building on the low-rank prior, ReLaSP introduced methods that impose $\ell_{1}$ regularization on the graph domain spectrum, encouraging each $\mathbf{F}_{\mathcal{G}}(:,j)$ to be sparse [46]. Although they only constrained a single domain, this approach can easily be extended to the graph and time domains framework, where $f_{c}(\mathbf{F}_{\mathcal{G}},\mathbf{F}_{\mathcal{T}})=||\mathbf{F}_{\mathcal{G}}||_{1}+||\mathbf{F}_{\mathcal{T}}||_{1}$. In this case, each $\mathbf{F}_{\mathcal{G}}(:,j)$ and each $\mathbf{F}_{\mathcal{T}}(:,i)$ tend to be sparse.

On the other hand, some studies introduced a smoothness prior in addition to the low-rank prior. LIMC constrains the smoothness constraints of FTVGS in both the time and graph domains using a modified second-order TV [47]. Similarly, LRGTS aims to leverage second-order smoothness information of FTVGS in both domains, denoting the smoothing in the graph domain by $\text{tr}(\hat{\mathbf{X}}^{T}\mathbf{L}\hat{\mathbf{X}})$ [30]. LRDS defines the smoothness of the temporal differences of FTVGS as $\text{tr}((\hat{\mathbf{X}}\mathbf{D}_{1})^{T}\mathbf{L}\hat{\mathbf{X}}\mathbf{D}_{1})$ and uses it as the smoothing regularization [29].

Since no existing research has designed optimization problems that combine $f_{r}(\hat{\mathbf{X}})$ with both $f_{c}(\mathbf{F}_{\mathcal{G}},\mathbf{F}_{\mathcal{T}})$ and $f_{s}(\hat{\mathbf{X}})$, we present an optimization problem that simultaneously incorporates all three regularizations for a more comprehensive comparison and analysis. Combined with the reweighting method [35, 48], we express the regularization function of the optimization problem as

[TABLE]

where $g(x)=\log(x^{1/2}+1)$, $\odot$ denotes the element-wise multiplication of two matrices, $\mathbf{W}_{\mathcal{G}}$ and $\mathbf{W}_{\mathcal{T}}$ are adjustable weight matrices for sparsity terms.

In considering smoothing regularization, it is found that in practice, FTVGS tends to exhibit stronger smoothness in the time domain and less smoothness in the graph domain [28]. Therefore, we adopt a temporal smoothing prior. It is noteworthy that Eq. 4 does not entirely disregard the low graph-domain smoothness. The term $||\mathbf{F}_{\mathcal{G}}||_{1}$ makes the graph signal band-limited at each instance, implying inherent correlations among signals at vertices. This implicitly constrains graph-domain smoothness in some sense. The solution to the optimization problem based on Eq. 4 is detailed in APPENDIX D, and the resulting algorithm is summarized in Algorithm 2.

Overall, since traditional MC algorithms cannot handle the problem of structurally missing samples, we delve deeper into the prior information of FTVGS. Starting from basic assumptions, we classify the priors into three types and organize them into an LSSP reconstruction framework. Within this framework, we first introduce some existing reconstruction algorithms and then present Eq. 4, which incorporates all three kinds of regularizers.

VI Experiment

In this section, we verify the feasibility of reconstructing the FTVGS from samples indexed by $\mathcal{S}$ by applying different reconstruction methods to real-world datasets. Furthermore, we analyze different algorithms within the LSSP framework based on the experimental results.

In the sampling phase, we follow the subset random sampling scheme (Algorithm 1), which ensures that no sample is observed from $\mathbf{X}(\mathcal{V}_{\mathcal{G}}\backslash\mathcal{I},\mathcal{V}_{\mathcal{T}}\backslash\mathcal{J})$. For each $\mathbf{X}$, we sample at different ratios, where $\alpha_{RC}$ represents the sampling ratio of randomly selected rows or columns, i.e.,

[TABLE]

where ${\rm Round}(\cdot)$ represents rounding numbers to the nearest integer. $\alpha_{\mathbf{X}_{RC}}$ represents the ratio of uniform sampling on $\mathbf{X}_{RC}$, i.e., $|\mathcal{S}|={\rm Round}(\alpha_{\mathbf{X}_{RC}}|\mathcal{I}||\mathcal{J}|)/(|\mathcal{I}||\mathcal{J}|)$. And $\alpha_{total}=|\mathcal{S}|/(NT)$ is the total sampling ratio.

For reconstruction in all experiments, the parameters of Algorithm 2 were set as follows: $\zeta=0.001,\rho_{1}=\rho_{2}=\rho_{3}=0.01$, $\mu_{1}^{(1)}=\mu_{2}^{(1)}=\mu_{3}^{(1)}=\frac{1}{2\max(\sigma(\mathbf{X}_{\mathcal{S}}))}$, where $\sigma(\mathbf{X}_{\mathcal{S}})$ is the singular values of $\mathbf{X}_{\mathcal{S}}$.

VI-A Competing Algorithms and Performance Metrics

We compare the Algorithm 2 with five algorithms: ICURC [32], MC-NC [36], ReLaSP [46], LIMC [47], and LRDS [29].

The ICURC algorithm integrates error minimization, low-rank constraints, and CUR decomposition. In a recent study on random sampling theory, CCS [32], ICURC demonstrated strong performance in reconstructing the original signal. Thus, we include it as a comparison algorithm.

The four algorithms, MC-NC, ReLaSP, LIMC, and LRDS, listed in Table II, all leverage partial prior information of the signals to formulate the optimization problem. Optimization problem based on Algorithm 2 is constructed within the LSSP framework, incorporating all three types of regularities.

To compare the reconstruction effects of different algorithms, we define the normalized root mean square error (NRMSE) to quantify the error between signals $\mathbf{X}_{a}$ and $\mathbf{X}_{b}$,

[TABLE]

VI-B Test on Synthetic Data

We generate a $128\times 128$ FTVGS $\mathbf{X}$ by $\mathbf{X}=\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{J}}\Psi_{\mathcal{T}}^{T}$, where $\Psi_{\mathcal{G}}$ and $\Psi_{\mathcal{T}}$ are given GFT basis and FT basis, respectively. As shown in Fig. 3 (a), $\mathbf{F}_{\mathcal{J}}$ represents the generated JFT spectrum with 100 nonzero rows, i.e., 100 nonzero low-frequency spectra. Each non-zero spectral component (each row) in $\mathbf{F}_{\mathcal{J}}$ is randomly generated, with the low-frequency bandwidth randomly selected between 28 and 88. Subsequently, Fig. 3 (b) displays $\mathbf{F}_{\mathcal{T}}=\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{J}}$, demonstrating that the spectrum of each temporal signal $\mathbf{X}(i,:)$ is bandlimited. In Fig. 3 (c), $\mathbf{F}_{\mathcal{G}}=\mathbf{F}_{\mathcal{J}}\Psi_{\mathcal{T}}^{T}$ is also constrained to 100 non-zero rows, indicating that the graph spectrum of each graph signal $\mathbf{X}(:,j)$ is sparse. All visualizations represent the squared absolute values of the elements.

We generate a total of 100 FTVGS. Under different sampling ratio settings, we calculate the mean NRMSE for various reconstruction algorithms, as compared in Table III.

VI-C Test on Real-world Data

The traffic data METR-LA is collected from loop detectors in the highways of Los Angeles County [49]. We chose data collected by 207 sensors from March 1st, 2012 to June 27th, 2012 for this experiment. Each sensor is regarded as a vertex $v_{i}\in\mathcal{V}_{\mathcal{G}}$, and we denote the data of every 512 consecutive timesteps as $\mathbf{X}$ for simulation. The METR-LA is divided into 100 FTVGS in total. That is, each $\mathbf{X}$ is an FTVGS with $N=207$ vertices and $T=512$ timesteps.

For the METR-LA dataset, the energy of $\mathbf{F}_{\mathcal{T}}$ and $\mathbf{F}_{\mathcal{G}}$ are shown in Fig. 4 (a) and (b), respectively. From Fig. 4 (a), each sensor acquired signal $\mathbf{X}(i,:)$ has a spectrum $\mathbf{F}_{\mathcal{T}}(i,:)$ that tends to be low-frequency, which implies both temporal smoothness and spectral sparsity. Furthermore, Fig. 4 (b) reveals that the graph spectrum of each graph signal $\mathbf{X}(:,j)$ also exhibits sparsity but does not demonstrate strong smoothness.

Under different sampling ratio settings, we calculate the average NRMSE for various reconstruction algorithms, as compared in Table IV

To investigate the trade-off between sampling rate and reconstruction accuracy, we conducted additional experiments under a wider range of total sampling ratios. The results for different methods are summarized in Fig 5.

It can be seen that in the task of reconstructing FTVGS from samples with structural missing, LSSP can achieve better results. Focusing on the results of the LSSP, the reconstruction error increases steadily as the sampling rate decreases from $85.82\%$ to $51.37\%$. However, when the total sampling rate further drops from $42.12\%$ to $34.29\%$, the NRMSE remains almost unchanged (0.2119 to 0.2121). This indicates that reducing the sampling rate to around $34.29\%$ does not significantly degrade reconstruction performance, making it a practical choice for costefficient sampling. On the other hand, further reductions below $34.29\%$ result in a sharp increase in error, suggesting that this point serves as a boundary beyond which reconstruction quality deteriorates rapidly. Therefore, we recommend sampling rates of approximately $\alpha_{RC}=70\%,\alpha_{\mathbf{X}_{RC}}=70\%,\alpha_{total}=34.29\%$ as the optimal balance between sampling cost and reconstruction accuracy. We visualize the comparison between the original FTVGS and the signal reconstructed by different algorithms in Fig. 6

VI-D Discussion

Now, we analyze the experimental results in relation to the FTVGS properties shown in Figs. 3 and 4.

$f_{r}(\mathbf{X})$ only: The results in Tables III and IV demonstrate that low-rank matrix approximation models are ineffective in handling data with structural missing. Algorithms like ICURC and MCNC, which rely solely on the low-rank prior of the FTVGS, cannot reconstruct the elements in $\mathcal{V}_{\mathcal{G}}\backslash\mathcal{I}$ and $\mathcal{V}_{\mathcal{T}}\backslash\mathcal{J}$. 2. 2.

$f_{r}(\mathbf{X})+f_{s}(\mathbf{X})$: LRDS and LIMC use two different smoothing priors, respectively. The smoothing regularizer of LRDS does not take into account the time domain, i.e., it fails to reconstruct the elements in $\mathcal{V}_{\mathcal{T}}\backslash\mathcal{J}$. LIMC uses the modified second-order TV regularizer, which takes into account the information of the neighboring entries, and thus gives a better result in reconstructing both the missing rows and columns. 3. 3.

$f_{r}(\mathbf{X})+f_{c}(\mathbf{F}_{\mathcal{G}},\mathbf{F}_{\mathcal{T}})$: From Tables III and IV, we can see that each column of $\mathbf{F}_{\mathcal{G}}$ is sparse and each row of $\mathbf{F}_{\mathcal{T}}$ is either low-frequency or tends to be low-frequency, indicating that $\mathbf{F}_{\mathcal{T}}(i,:)$ is also sparse. In the ReLaSP experiments, $\mathbf{F}_{\mathcal{G}}$ and $\mathbf{F}_{\mathcal{T}}$ are used to constrain the correlations between rows and columns, respectively. This allows missing rows and columns to be reconstructed through these correlations. 4. 4.

$f_{s}(\mathbf{X})$ v.s. $f_{c}(\mathbf{F}_{\mathcal{G}},\mathbf{F}_{\mathcal{T}})$: Algorithms like LIMC reconstruct the FTVGS based on a smoothing prior. However, as shown in Fig. 4 (b), the FTVGS does not exhibit strong smoothness in the graph domain. Each $\mathbf{F}_{\mathcal{G}}(:,j)$ is only sparse rather than low-rank, meaning its energy is not concentrated solely in the low-frequency (smaller $\lambda$) region. Consequently, LIMC may over-smooth the reconstructed signal in the vertex domain. In contrast, ReLaSP achieves better reconstruction results by imposing a sparse regularity constraint on the graph spectrum $\mathbf{F}_{\mathcal{G}}$. 5. 5.

LSSP: Based on 4), Eq. 4 selects $||\mathbf{F}_{\mathcal{G}}||_{1}$ as the regularizer. Additionally, as shown in Fig. 4 (a), the FTVGS is low-frequency in the time domain, but real-world signals tend to exhibit sparsity even within low-frequency bands. Therefore, for the regularity of $\mathbf{F}_{\mathcal{T}}$, we incorporate both a smoothing regularizer $||\hat{\mathbf{X}}\mathbf{D}_{2}||_{F}^{2}$ and sparsity regularizer $||\mathbf{F}_{\mathcal{T}}||_{1}$. Ultimately, Algorithm 2 achieves better reconstruction results when addressing structurally missing samples. This demonstrates that, within the LSSP framework, leveraging signal priors to refine the optimization problem enhances reconstruction accuracy.

VII Conclusion

In this work, motivated by practical applications, we propose a subset random sampling scheme for FTVGS, which involves uniform sampling from the submatrix $\mathbf{X}_{RC}$ instead of the entire matrix $\mathbf{X}$. Theoretically, we prove sufficient conditions that guarantee the reconstruction of the original FTVGS with high probability under the subset random sampling scheme. Additionally, for data with randomly missing elements and entire rows and columns, we introduce the targeted reconstruction algorithm LSSP. We validate the effectiveness of the LSSP algorithm on a real-world dataset.

Appendix A Proof of the Lemma 1

We first prove that $\mathrm{rank}(\mathbf{X}_{R})=r$ holds with probability at least $1-\delta$ when we sample $|\mathcal{I}|$ rows of the signal $\mathbf{X}$ uniformly without replacement. The lower bound of $|\mathcal{I}|$ is related to $r$ and $\mu_{1}(\mathbf{X})$ as shown in Eq. 5. From the thin SVD of $\mathbf{X}$ we know that $\mathbf{U}\in\mathbb{R}^{N\times r}$ is a matrix with orthogonal columns. Selecting the rows of $\mathbf{X}$ is essentially row sampling on $\mathbf{U}$.

Theorem 4.1 and Corollary 4.2 in [50] proved that $\mathrm{rank}(\mathbf{U}(\mathcal{I},:))=r$ holds with probability at least $1-\delta$ if

[TABLE]

for some $0<\delta<1$ and $0<\epsilon<1$, where $\xi_{1}=||\mathbf{U}||_{2,\infty}^{2}$.

According to 2, $\mu_{1}(\mathbf{X})r/N\geq||\mathbf{U}||_{2,\infty}^{2}=\xi_{1}$, so we have

[TABLE]

Further, we have $\mathbf{X}_{R}=\mathbf{U}(\mathcal{I},:)\Sigma\mathbf{V}^{T}$, where $\Sigma\mathbf{V}^{T}$ is a matrix with orthogonal rows. Then $\mathrm{rank}(\mathbf{X}_{R})=r$ holds with probability at least $1-\delta$.

We secondly prove that $\mathrm{rank}(\mathbf{X}_{RC})=r$ holds with probability at least $(1-\delta)^{2}$ when we sample $|\mathcal{J}|$ columns of the signal $\mathbf{X}_{R}$ uniformly without replacement. The lower bound of $|\mathcal{J}|$ is related to $r$ and $\mu_{2}(\mathbf{X})$ as shown in Eq. 6.

The thin SVD of $\mathbf{X}_{R}$ is $\mathbf{X}_{R}=\mathbf{U}_{R}\Sigma_{R}\mathbf{V}^{T}_{R}$, as derived in Eq. 1. The same argument yields $|\mathcal{J}|\geq 3T\xi_{2}\ln{(2r/\delta)}/\epsilon^{2}$ for some $0<\delta<1$ and $0<\epsilon<1$, where $\xi_{2}=||\mathbf{V}_{R}||_{2,\infty}^{2}$ [50].

We have

[TABLE]

where $\tilde{\mathbf{V}}_{R}$ is an orthogonal matrix, for any $i$, $||\mathbf{V}(i,:)\tilde{\mathbf{V}}_{R}||_{2}=||\mathbf{V}(i,:)||_{2}$ and (a) holds. So we get

[TABLE]

Since $\mathbf{U}(\mathcal{I},:)\Sigma$ is a matrix with orthogonal columns with probability at least $1-\delta$. Then $\mathrm{rank}(\mathbf{X}_{RC})=r$ holds with probability at least $(1-\delta)^{2}$.

Appendix B Proof of the Lemma 2

We demonstrate that the submatrix $\mathbf{X}_{RC}$ exhibits similar incoherence. Using $\mathbf{U}_{R}$ and $\mathbf{V}_{R}$ as intermediaries, we establish upper bounds on $\mathbf{U}_{RC}$ and $\mathbf{V}_{RC}$, illustrating how the properties of the original FTVGS $\mathbf{X}$ are transferred to its submatrix $\mathbf{X}_{RC}$.

According to the thin SVD of $\mathbf{X}_{R}$, Lemma 4.2 and Lemma 4.3 in [51] give

[TABLE]

According to Lemma 3.4 in [52],

[TABLE]

with failure probability at most $p=r\left[\frac{e^{-\eta}}{(1-\eta)^{1-\eta}}\right]^{\log r},\eta\in[0,1)$.

Then we substitute Eq. 8 into Eq. 7 to get

[TABLE]

holds with probability at least $1-p$. Eq. 2 yields

[TABLE]

where $\tilde{\mathbf{U}}_{RC}$ is an orthogonal matrix, for any $i$ we have $||\mathbf{U}_{R}(i,:)\tilde{\mathbf{U}}_{RC}||_{2}=||\mathbf{U}_{R}(i,:)||_{2}$ and (b) holds.

Next, we consider $||\mathbf{V}_{RC}||_{2,\infty}$. According to Eq. 2, we have $\mathbf{V}_{RC}^{T}=\Sigma_{RC}^{-1}\mathbf{U}_{RC}^{T}\mathbf{X}_{RC}$. Thus, for any $i$,

[TABLE]

where $\mathbf{e}_{i}$ is a unit vector whose $i$-th element is 1, $\mathbf{U}_{RC}$ is a column orthogonal matrix and $\sigma_{1}(\mathbf{U}_{RC})=1$. Besides $\mathbf{X}_{RC}=\mathbf{X}(:,\mathcal{J})=\mathbf{U}_{R}\Sigma_{R}\mathbf{V}_{R}(\mathcal{J},:)^{T}$, so

[TABLE]

where $\mathbf{U}_{R}$ is a column orthogonal matrix and $\sigma_{1}(\mathbf{U}_{R})=1$. Then combining Eq. 9 and Eq. 10, we have

[TABLE]

where Lemma 4.2 in [51] gives (c).

Then, according to $\mathbf{X}_{RC}=\mathbf{U}(\mathcal{I},:)\Sigma\mathbf{V}(\mathcal{J},:)^{T}$ and $||\Sigma^{\dagger}||_{2}=\frac{1}{\sigma_{\mathrm{min}}(\mathbf{X})}=||\mathbf{X}^{\dagger}||_{2}$, we have

[TABLE]

Thus, we substitute Eq. 12 into Eq. 11 to get

[TABLE]

with probability at least $(1-p)^{2}$, where Lemma 3.4 in [52] gives (d).

Appendix C Proof of the Theorem 1

At last, we provide a lower bound on the number of samples (Theorem 1), which depends only on the submatrix size and properties of $\mathbf{X}$. We demonstrate that the sampling process, including the selection of rows/columns to obtain the submatrix and random sampling within $\mathbf{X}_{RC}$, allows for successful reconstruction.

According to the Cauchy-Schwartz inequality, we have

[TABLE]

and for all $(i,j)$

[TABLE]

Since $\mathcal{I}$ and $\mathcal{J}$ are chosen uniformly from $\mathcal{V}_{\mathcal{G}}$ and $\mathcal{V}_{\mathcal{T}}$ based on Algorithm 1, combining Eq. 13 and Lemma 2 we have

[TABLE]

for $\eta\in[0,1)$ with probability at least $(1-p)^{2}$.

Let $\mu_{0}=\frac{(\kappa(\mathbf{X}))^{2}}{(1-\eta)^{3/2}}\sqrt{\frac{\mu_{1}(\mathbf{X})\mu_{2}(\mathbf{X})Nr^{2}}{|\mathcal{I}|}}$. The matrix $\mathbf{U}_{RC}\mathbf{V}_{RC}^{T}$ has a maximum entry bounded by $\mu_{0}\sqrt{r/(|\mathcal{I}||\mathcal{J}|)}$. Let $n=\max{\{|\mathcal{I}|,|\mathcal{J}|\}}$. Applying $\mu_{0}$ and $n$ to Theorem 2 in [39], we obtain that if

[TABLE]

for some $\beta>1$, $\mathbf{X}_{RC}$ can be reconstructed with probability

[TABLE]

Based on 1, we know that $||\triangledown\mathbf{X}(i,j)||<C$, so $|\mathbf{Q}^{T}\mathbf{X}(i,j)|<\infty$ and $|\mathbf{X}\mathbf{D}_{1}(i,j)|<\infty$. We define that

[TABLE]

for some positive parameter $a$, then

[TABLE]

Similarly, we have

[TABLE]

and

[TABLE]

According to Eq. (5.7) in [45], we have

[TABLE]

Then, combining 1 and Eqs. 15, 16 and 17 we have

[TABLE]

Thus, we can reconstruct $\mathbf{X}$ from $\mathbf{X}_{RC}$ by solving the TV inpainting model according to Theorem 4.4 in [53].

As a result, if Eq. 14 holds, we are able to reconstruct $\mathbf{X}$ from $\mathcal{S}$ with probability at least

[TABLE]

Appendix D Algorithm for opt. prob. based on Eq. 4

In Eq. 4, the signal spectral coefficients are constrained using the $\ell_{1}$-norm. Unlike the non-convex $\ell_{0}$-norm, the $\ell_{1}$-norm is convex, making the optimization problem easier to solve. However, since the $\ell_{1}$-norm is defined as the sum of the absolute values of non-zero elements, it imposes sparsity constraints but also affects the magnitude of these elements. To overcome the influence of the signal magnitudes, a reweighting scheme is introduced [35, 48].

To enhance robustness and accelerate the solution of the optimization problem, we incorporate the error matrix $\mathbf{E}$ into the constraints [54]. As a result, the optimization problem based on Eq. 4 is

[TABLE]

where $\mathbf{E}$ is an auxiliary error matrix and $\mathscr{P}_{\mathcal{S}}(\cdot)$ be the projection operator onto $\mathcal{S}$,

[TABLE]

The LSSP algorithm consists of three layers of loops. The outermost loop is the reweighting framework, with $p$ denoting the iteration index. For example, the variable $\mathbf{F}_{\mathcal{G}}$ at the $p$-th iteration is represented as $\mathbf{F}_{\mathcal{G}}^{\{p\}}$. The middle loop updates each variable, with $k$ denoting the iteration index, and the result of the $k$-th iteration is $\mathbf{F}_{\mathcal{G}}^{(k)}$. During each iteration of the middle loop, the updates of $\mathbf{F}_{\mathcal{G}}$ and $\mathbf{F}_{\mathcal{T}}$ are required to solve the optimization problems Eq. 21 and Eq. 25, leading to the innermost loop. In this loop, indexed by $q$, $\mathbf{F}_{\mathcal{G}}^{(k),q}$ represents the $q$-th iteration of $\mathbf{F}_{\mathcal{G}}$ during the $k$-th variable update.

The reweighting scheme is an effective iterative method that significantly enhances signal reconstruction by reinforcing the sparsity of the coefficients. The key idea behind reweighted $\ell_{1}$ norm minimization is to apply large weights to suppress zero or near-zero elements and small weights to promote non-zero entries. As suggested in [48], the weights in reweighted $\ell_{1}$ norm minimization are inversely proportional to the magnitude of the coefficients. The weight matrix is updated as follows

[TABLE]

for parameter $\zeta>0$.

In the $p$-th iteration, given the weight matrices $\mathbf{W}_{\mathcal{G}}^{\{p\}}$ and $\mathbf{W}_{\mathcal{T}}^{\{p\}}$, we solve Eq. 4 subject to the equation constraints.

The non-convex regularizer is optimized using an iteratively reweighted least squares method [36]. In each iteration $k$, A weight matrix $\mathbf{R}^{(k)}$ is updated based on $\hat{\mathbf{X}}^{(k-1)}$, and then $\hat{\mathbf{X}}^{(k)}$ is updated according to $\mathbf{R}^{(k)}$ .

By applying the augmented Lagrangian method, we transform Eq. 4 into an unconstrained optimization problem, where the Lagrangian function is expressed as:

[TABLE]

where $\mathbf{Y}_{1}$, $\mathbf{Y}_{2}$ and $\mathbf{Y}_{3}$ are the Lagrange multiplier matrices, $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$ are the penalty factors, and $\mathbf{X}_{\mathcal{S}}$ is the incomplete observed matrix.

To efficiently solve Eq. 20 with multiple unknown variables, we adopt the alternating direction method of multipliers (ADMM), optimizing each variable while keeping the others fixed.

Optimizing $\mathbf{F}_{\mathcal{G}}$ means solving the following optimization problem:

[TABLE]

Instead of solving Eq. 21 directly, we estimate $\mathbf{F}_{\mathcal{G}}$ using the following approximation model. Let $f_{1}(\mathbf{F}_{\mathcal{G}})=\langle\mathbf{Y}_{1}^{(k)},\hat{\mathbf{X}}^{(k)}-\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{G}}\rangle+\frac{\mu_{1}^{(k)}}{2}||\hat{\mathbf{X}}^{(k)}-\Psi_{\mathcal{G}}\mathbf{F}_{\mathcal{G}}||_{F}^{2}$. By introducing a proximal variable $\mathbf{A}_{1}$ [55], we have

[TABLE]

where $\tau_{1}>0$ is a constant, and

[TABLE]

Now we can solve the $\text{prox}_{\tau_{1}}(\mathbf{A}_{1})$ by nonuniform soft thresholding on $\mathbf{A}_{1}-\tau_{1}\triangledown f_{1}(\mathbf{A}_{1})$ [35], i.e.,

[TABLE]

where

[TABLE]

is the shrinkage operator. Then, the proximal variable $\mathbf{A}$ can be update by

[TABLE]

where $b_{1}^{q}$ is a positive parameter with $b_{1}^{1}=1$. The solution to this proximal problem is $\mathbf{F}_{\mathcal{G}}^{(k+1)}$, i.e., $\mathbf{F}_{\mathcal{G}}^{(k+1)}=\mathbf{F}_{\mathcal{G}}^{(k),Q}$

Consider $\mathbf{F}_{\mathcal{T}}$, the optimization problem is

[TABLE]

where $f_{2}(\mathbf{F}_{\mathcal{T}})=\langle\mathbf{Y}_{2}^{(k)},(\hat{\mathbf{X}}^{(k)})^{T}-\Psi_{\mathcal{T}}\mathbf{F}_{\mathcal{T}}\rangle+\frac{\mu_{2}^{(k)}}{2}||(\hat{\mathbf{X}}^{(k)})^{T}-\Psi_{\mathcal{T}}\mathbf{F}_{\mathcal{T}}||_{F}^{2}$. In like manner, $\mathbf{F}_{\mathcal{T}}^{(k+1)}$ can be approximated as follows:

[TABLE]

where $\mathbf{A}_{2}$ is a proximal variable, $\tau_{2}>0$ is a constant, $b_{2}^{q}$ is a positive parameter with $b_{2}^{1}=1$, and

[TABLE]

Eventually, we get $\mathbf{F}_{\mathcal{T}}^{(k+1)}=\mathbf{F}_{\mathcal{T}}^{(k),Q}$.

To solve $\hat{\mathbf{X}}$, we first eigenvalue decompose $h(\hat{\mathbf{X}})$ as

[TABLE]

Then $\mathbf{R}$ is updated by

[TABLE]

where $g^{\prime}(x)=\frac{1}{2\sqrt{x}(\sqrt{x}+1)}$. Let

[TABLE]

By taking the derivative with respect to $\hat{\mathbf{X}}$ and setting it to zero, we have

[TABLE]

The optimal solution of $\hat{\mathbf{X}}$ can be obtained by

[TABLE]

To solve $\mathbf{E}$, the Lagrangian function becomes

[TABLE]

In the same way, we take the derivative with respect to $\mathbf{E}$ and set it to zero

[TABLE]

When updating the error matrix $\mathbf{E}$, we only need to update the entries in $\mathcal{S}^{c}$ (the complement of $\mathcal{S}$). Thus, we have

[TABLE]

where 0 is an all-zero matrix.

Furthermore, the Lagrange multiplier matrices and penalty factors can be updated by

[TABLE]

where $\rho_{1}>1$, $\rho_{2}>1$, and $\rho_{3}>1$ are constants that ensure the increase of the penalty factors.