Guaranteed Multidimensional Time Series Prediction via Deterministic Tensor Completion Theory

Recovery theory for matrix/tensor completion in random cases has been extensively studied, yielding significant advancements in many fields [38] [39] [40] [41]. In contrast, deterministic observation patterns have garnered considerably less attention. However, the pattern of missing entries depends on the specific problem context in practice, which can lead to highly non-random occurrences. A case in point is the presence of missing entries in images, which is often attributed to occlusions or shadows present in the scene. To investigate the recovery theory for matrix/tensor completion under non-random missing conditions, previous researchers have conducted extensive exploration using various tools [28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. For instance, Singer et al. [29] apply the basic ideas and tools of rigidity theory to determine the uniqueness of low-rank matrix completion. Based on Plücker coordinates, Tsakiris [31] provides three families of patterns of arbitrary rank and arbitrary matrix size that allow unique or finite number of completions. Besides, Liu et al. [36] propose isomeric condition and relative well-conditionedness to ensure that any matrix can be recovered from sampling of matrix terms. Despite these advancements, most of the current theoretical work focuses on deterministic matrix completion, while the tensor case is not explored enough. Therefore, the theoretical framework for deterministic tensor completion still requires significant enhancement.

Matrix and tensor completion methods have garnered significant attention for time series forecasting [23] [26] [42]. These approaches generally fall into two categories: decomposition-based methods and nuclear norm-based methods. Matrix/tensor decomposition, widely used for collaborative filtering, provide a natural solution for handling missing data in time series forecasting [42] [43] [44]. These models approximate the incomplete time series using bilinear or multilinear decompositions with predefined rank parameters. To capture temporal dynamics, recent studies have introduced temporal constraint functions, such as autoregressive regularization terms, to regularize the temporal factor matrix [45] [23] [44]. Based on this, Takeuchi et al. [22] model spatiotemporal tensor data by introducing spatial autoregressive regularizer, providing additional predictive power in the spatial dimension.

As opposed to decomposition-based methods, nuclear norm-based approaches do not require predefined rank parameters [38] [40] [41]. However, directly applying nuclear norms to time series often yields suboptimal predictions [6]. A practical enhancement involves applying structural transformations to time series data to create structured matrices, enabling the application of nuclear norms for forecasting. For example, Gillard et al. [26] and Butcher et al. [27] used the Hankel transform to convert univariate time series into Hankel matrices and applied nuclear norm minimization for prediction. Similarly, Liu et al. [24] [25] utilized convolution matrices for encoding nuclear norms, providing insights into the predictive domain that can be exactly forecasted. Although these methods offer predictive insights, they involve converting time series into large matrices, resulting in the loss of inherent tensor structure and imposing high computational costs. To address these limitations, we propose the Temporal Convolution Tensor Nulclear Norm (TCTNN) model, which preserves the tensor structure, offers theoretical guarantees for predictions, and achieves superior performance in terms of both accuracy and computational efficiency.

We begin by providing a brief overview of the key concepts within the tensor Singular Value Decomposition (t-SVD) framework [46] [47] [48] [49] [50]. This includes the tensor-tensor product, tensor nuclear norm, and the conditions for tensor incoherence. Following this, we propose a novel deterministic sampling metric designed to establish an exact recovery theory for deterministic tensor completion.

For consistency, we use lowercase, boldface lowercase, capital, and Euler script letters to represent scalars, vectors, matrices, and tensors, respectively, e.g., $x\in\mathbb{R}$, $\bm{x}\in\mathbb{R}^{m}$, $X\in\mathbb{R}^{m_{1}\times m_{2}}$, and $\mathcal{X}\in\mathbb{R}^{m_{1}\times m_{2}\times\cdots\times m_{d}}$. For an order-$d$ tensor $\mathcal{X}$ of dimensions $m_{1}\times m_{2}\times\cdots\times m_{d}$, its elements, slices, and sub-tensors are represented as $\mathcal{X}_{(...)}$, such as $\mathcal{X}_{(i_{1},i_{2},\cdots,i_{d})}$ represents the $(i_{1},i_{2},\cdots,i_{d})$-th element and $\mathcal{X}_{(:,:,i_{3},\cdots,i_{d})}$ represents the $(i_{3},\cdots,i_{d})$-th face slice. To construct the t-SVD framework, some operators on $\mathcal{X}\in\mathbb{R}^{m_{1}\times m_{2}\times\cdots\times m_{d}}$ are summarized in Table I.

Since the multiplication of $\operatorname{bcirc}(\mathcal{A})$ and $\operatorname{bunfold}(\mathcal{B})$ is in the form of circular convolution, t-product can also be defined using circular convolution:

where $\star$ denotes the operator of circular convolution. Evidently, the t-product of tensors is an operation analogous to matrix multiplication. As a result, many properties of matrix multiplication can be extended to the t-product [48] [50].

In addition, the discrete Fourier transform (DFT) can convert the circular convolution operation into an element-wise product. Therefore, the t-product $\mathcal{A}*\mathcal{B}$ can be obtained through DFT. For a given tensor $\mathcal{X}$, its Fourier-transformed in last d-2 dimensions is expressed as

where $\times_{j}$ denotes the j-mode product [52] and $\mathrm{F}_{m_{j}}$ represents the DFT matrices of size $m_{j}\times m_{j}$ for $j=3,\cdots,d$. Then

where $\Delta$ denotes the face-wise product ($\mathcal{Z}=\mathcal{X}\Delta\mathcal{Y}\Leftrightarrow\mathcal{Z}_{(:,:,i_{3},\cdots,i_{d})}=\mathcal{X}_{(:,:,i_{3},\cdots,i_{d})}\mathcal{Y}_{(:,:,i_{3},\cdots,i_{d})}$ for all face slices).

Now we consider deterministic tensor completion problem, a broader research challenge. Let $\mathcal{M}\in\mathbb{R}^{m_{1}\times\cdots\times m_{d}}$ represent an unknow target tensor with tubal rank $r$ and skinny t-SVD $\mathcal{M}=\mathcal{U}*\mathcal{S}*\mathcal{V}^{\mathrm{T}}$. Suppose that we observe only the entries of $\mathcal{M}$ over a deterministic sampling set $\Omega\subseteq\left[m_{1}\right]\otimes\left[m_{2}\right]\otimes\cdots\otimes\left[m_{d}\right]$. The corresponding mask tensor is denoted by $\bar{\Omega}$, with the following entry definition:

The sampling operator $\mathcal{P}_{\Omega}$ is then defined as $\mathcal{P}_{\Omega}(\mathcal{M})=\bar{\Omega}\circ\mathcal{M},$ where $\circ$ denotes the Hadamard product.

To address the deterministic tensor completion problem, we employ a standard constrained Tensor Nuclear Norm (TNN) model:

where $\Omega$ is a deterministic sampling set. In the context of tensor completion under random sampling, it is often assumed that the data sampling follows a specific distribution, such as the Bernoulli distribution, then the sampling probability $p$ can naturally serve as a measure of sampling [48]. However, this approach is clearly ineffective for deterministic sampling. Inspired by the minimum row/column sampling ratio for matrix deterministic sampling [36], we introduce a tensor deterministic sampling metric $\rho$, termed the minimum horizontal/lateral sub-tensor sampling ratio. Here, the horizontal/lateral sub-tensors of a tensor are analogous to the rows/columns of a matrix, respectively. As mentioned above, the tensor t-product, spectral/nuclear norm within the t-SVD framework are defined in a matrix-like manner, providing a solid basis for defining tensor deterministic sampling metric by following matrix case.

For an intuitive understanding, we provide a schematic diagram illustrating horizontal/lateral mask sub-tensor sampling number in the three-dimensional case, as shown in Fig.2.

Tensor incoherence is a crucial theoretical tool in low-rank tensor recovery [38] [39] [50] [53]. It enforces or constrains the low-rank structure of the underlying tensor to avoid ill-posed problems, such as when most elements of $\mathcal{M}$ are zero, yet it is still low-rank.

We briefly show the exact recovery guarantee of deterministic low-rank tensor completion problem.

The above result demonstrates that minimizing the tensor nuclear norm can achieve exact deterministic tensor completion, provided that condition (13) is satisfied. The feasibility of exact recovery hinges upon the interplay between the sampling set and the tensor rank. As the low-rank structure of the tensor strengthens, the minimum slice sampling ratio required for exact deterministic tensor completion decreases correspondingly. Moreover, the proposed deterministic exact recovery theory can be extended to the random sampling case. The following corollary reveals that when the sampling set satisfies $\Omega\sim\operatorname{Ber}(p)$ and $p\geq 1-1/{2\mu r(r_{s}+1)}$, exact recovery holds with high probability.

Multidimensional time series prediction, as a specific instance of deterministic tensor completion, can be addressed using the aforementioned tensor nuclear norm minimization, which results in the following model:

where $\Omega=\left[t-h\right]\otimes\left[n_{1}\right]\otimes\cdots\otimes\left[n_{p}\right]$ represents the deterministic sampling set corresponding to the historical portion sample, $\mathcal{P}_{\Omega}(\cdot)$ is the projection operator. $\mathcal{M}\in\mathbb{R}^{t\times n_{1}\times n_{2}\times\cdots\times n_{p}}$ denotes the multidimensional time series to be completed.

Deterministic tensor completion theory can be applied to prediction tasks where the sampling set $\Omega$ is chosen as the historical data region in the prediction problem. However, in such cases, the minimum horizontal/lateral sub-tensor sampling ratio

which fails to meet the recovery conditions (13). This implies that tensor nuclear norm minimization cannot be directly applied to prediction tasks, as it would result in prediction outcomes that are entirely zero. The core issue lies in the concentrated nature of missing data in time series forecasting, which leads to insufficient sampling coverage. A natural solution to this problem is to implement a structural transformation that disperses the missing data, thereby increasing the minimum horizontal/lateral sub-tensor sampling ratio above zero. To address this, we introduce the temporal convolution tensor for multidimensional time series. This approach ensures that each horizontal/lateral sub-tensor contains a sufficient number of sampled elements, as illustrated in Fig.3.

In this study, circular convolution (III.2) is applied only along the time dimension of the multidimensional time series, a process referred to as temporal circular convolution. Following this, we obtain the temporal convolution tensor of the multidimensional time series.

To provide an intuitive understanding, we present a schematic diagram illustrating the conversion from a multidimensional time series (incomplete tensor) of size $4\times 3\times 2$ to the corresponding temporal convolution tensor of size $4\times 4\times 3\times 2$, as shown in Fig.3(a,b). Notably, the last horizontal sub-tensor of the incomplete tensor $\mathcal{P}_{\Omega}(\mathcal{M})$ consists entirely of zeros. However, after applying the temporal convolution transform $\mathcal{T}_{k}(\cdot)$, each horizontal and lateral slice of the incomplete temporal convolution tensor $\mathcal{T}_{k}(\mathcal{P}_{\Omega}(\mathcal{M}))$ contains a sufficient number of sampled entries. As depicted in Fig.3(b), the temporal convolution sampling set achieves a minimum horizontal/lateral sub-tensor sampling ratio

This significantly disperses the unobserved regions, enhancing the potential for axact predictions.

By applying transform $\mathcal{T}_{k}(\cdot)$ to the incomplete multidimensional time series $\mathcal{P}_{\Omega}(\mathcal{M})$, we obtain the incomplete temporal convolution tensor $\mathcal{P}_{{\Omega}_{\mathcal{T}}}(\mathcal{T}_{k}(\mathcal{M}))$ with a relatively dispersed distribution of missing values. In many cases, the completion of a tensor relies on its low-rankness, and this holds true for the temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})$ as well. Fortunately, The low-rankness of the temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})$ can often be satisfied in various scenarios, as it can be derived from the smoothness and periodicity of the time series $\mathcal{M}$. To reach a rigorous conclusion, we consider the rank-r approximation error of temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})\in\mathbb{R}^{t\times k\times n_{1}\times\cdots\times n_{p}}$, which is denoted by $\varepsilon_{r}(\mathcal{T}_{k}(\mathcal{M}))$ and defined as

where $\mathcal{Y}\in\mathbb{R}^{t\times k\times n_{1}\times\cdots\times n_{p}}$.

The notation $\left\lceil\cdot\right\rceil$ in aboves lemmas denotes the ceiling operation. Lemma IV.1 and Lemma IV.2 indicate that as the periodicity and smoothness of the original data matrix/tensor $\mathcal{M}$ along the time dimension increase, the low-rankness of its temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})$ becomes more pronounced. To illustrate this concept intuitively, we generated two multivariate time series: one exhibiting strong periodicity and the other characterized by temporal smoothness, as shown in Fig.4(a-1,a-2), respectively. Temporal convolution transform were applied to these time series matrices, and the resulting temporal convolution tensors are depicted in Fig.4(b-1,b-2). As shown in Fig.4(c-1,c-2), the tensor singular values of the temporal convolution tensors decrease rapidly. This trend highlights the enhanced low-rankness of the temporal convolution tensors, further substantiating the theoretical findings.

Given that the smoothness and periodicity of the original time series tensor $\mathcal{M}$ ensure the low-rank property of the temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})$, it is natural to leverage the convex relaxation of tensor rank—the tensor nuclear norm—to encode this low-rankness. This leads to the formulation of the following completion model:

to obtain $\hat{\mathcal{Y}}$ as the estimator of $\mathcal{T}_{k}(\mathcal{M})$. After that, we perform the inverse temporal convolution on $\hat{\mathcal{Y}}$, thereby obtaining $\mathcal{T}_{k}^{-1}(\hat{\mathcal{Y}})$ as the final preditior of $\mathcal{M}$.

This logical flow is depicted in Fig.3, outlining the steps to obtain the complete tensor from the incomplete tensor $\mathcal{P}_{\Omega}(\mathcal{M})$ in the prediction problem. The process consists of three main steps: 1.Transformation: Apply the temporal convolution transform $\mathcal{T}_{k}(\cdot)$ to the incomplete tensor $\mathcal{P}_{\Omega}(\mathcal{M})$, resulting in the incomplete temporal convolution tensor P$\mathcal{P}_{{\Omega}_{\mathcal{T}}}(\mathcal{T}_{k}(\mathcal{M}))$; 2. Completion: Solve model (20) to derive the complete temporal convolution tensor $\mathcal{T}_{k}(\mathcal{M})$ from the incomplete tensor $\mathcal{P}_{{\Omega}_{\mathcal{T}}}(\mathcal{T}_{k}(\mathcal{M}))$; 3. Inverse Transformation: Retrieve the complete tensor $\mathcal{M}$ by applying the inverse transform $\mathcal{T}_{k}^{-1}(\cdot)$. By consolidating these steps into an equivalent representation, we formulate the Temporal Convolution Tensor Nuclear Norm (TCTNN) model for multidimensional time series prediction, i.e.,

where $\Omega=\left[t-h\right]\otimes\left[n_{1}\right]\otimes\cdots\otimes\left[n_{p}\right]$ and $h$ is the forecast horizon.

We now develop the exact prediction gaurantee for time series data prediction by applying the deterministic tensor completion theory established in the preceding section. We first define the temporal convolution tensor incoherence conditions as follows:

Based on the temporal convolution tensor incoherence, we derive the exact prediction theory for the TCTNN model (21) from the deterministic tensor completion recovery theory.

Theorem IV.1 states that exact prediction is promising when the temporal convolution low-rankness is satisfied. As the low-rankness of the temporal convolution tensor increases and the prediction horizon $h$ decreases, achieving exact predictions becomes increasingly feasible. According to Lemmas IV.2 and IV.1, stronger periodicity and smoothness in the temporal dimension contribute to greater low-rankness in the temporal convolution tensor, thereby improving the likelihood of exact prediction.

In very recent, Liu et al. [25] propose the convolution matrix nulclear norm minimization (CNNM) model which can be formulated as

where $\mathbf{s}=[s_{t},s_{1},\cdots,s_{p}]$ is kernel size, $\mathcal{M}$ represents the time series to be completed, $\|\cdot\|_{*}$ denotes the nuclear norm of a matrix and $\mathcal{C}_{\mathbf{s}}(\cdot)$ denotes the transformation that converts the time series into the convolution matrix.

In contrast to the CNNM model, the proposed TCTNN model does not convert the original data tensor into convolution matrices along each dimension. Instead, we perform tensor convolution along the temporal dimension of the time series data. Our approach offers advantages in three aspects. 1. Tensor structure information: It preserves the tensor structure of multi-dimensional time series and avoids the information loss caused by converting tensor data into matrices. As demonstrated by the fact that $[\mathcal{T}_{k}(\mathcal{X})](:,1,:,:)=\mathcal{X}$, our temporal convolution tensors maintain the structure information of the multidimensional time series, unlike the convolution matrices, which satisfy $[\mathcal{C}_{\mathbf{s}}(\mathcal{X})](:,1)=vec(\mathcal{X})$. 2. Computational complexity: It avoids the computationally expensive task of processing the nuclear norm of very large matrices, significantly reducing the computational cost. For example, to predict a $(p+1)$-order multidimensional time series $\mathcal{M}\in\mathbb{R}^{n\times\cdots\times n}$, the computational complexity of each iteration of the CNNM model is $O(n^{3p+3})$, while the TCTNN model requires only $O(n^{p+3})$ per iteration. Clearly, the computational complexity of the former is nearly prohibitive. 3. Feature utilization: By performing convolution exclusively along the time dimension, our method is better suited for time series data. Many time series exhibit strong characteristics such as periodicity and smoothness only along the time dimension. Convolution along all dimensions could lead to improper feature capture, thereby reducing prediction accuracy.

We adopt the famous optimization framework alternating direction method of multipliers (ADMM) [54, 55, 56] to solve the model (21). First of all, we introduce the auxiliary variable $\mathcal{Y}=\mathcal{T}_{k}(\mathcal{X})$ to separate the temporal convolution transform $\mathcal{T}_{k}(\cdot)$, then the TCTNN model (21) can be rewritten as below:

For model (26), its partial augmented Lagrangian function is

where $\mathcal{N}$ are Lagrange multipliers, $\mu_{\ell}>0$ is a positive scalar, and $\mathcal{C}$ is only the multipliers dependent squared items. According to the ADMM optimization framework, the minimization problem of (27) can be transformed into the following subproblems for each variable in an iterative manner, i.e.,

and the multipliers is updated by

where ${\ell}$ denotes the count of iteration in the ADMM. In the following, we deduce the solutions for (28) and (29) respectively, each of which has the closed-form solution.

1) Updating $\mathcal{Y}^{{\ell}+1}$: Fixed other variables in (28), then it degenerates into the following tensor nuclear norm minimization with respect to $\mathcal{Y}$, i.e.,

the close-form solution of this sub-problem is given as

via the order-$(p+2)$ t-SVT as stated in Lemma III.2.

2) Updating $\mathcal{X}^{{\ell}+1}$: By deriving simply the KKT conditions, the closed-form solution of (29) is given by

where $\Omega^{\perp}$ is the complement of $\Omega$ and $\mathcal{T}_{k}^{-1}(\cdot)$ is the inverse operator of $\mathcal{T}_{k}(\cdot)$. The whole optimization procedure for solving proposed TCTNN model (21) is summarized in Algorithm 1.