# Spatio-Temporal Traffic Data Tensor Restoration Method Based on Direction Weighting and P-Shrinkage Norm

Authors: Jiangnan Wu; Yongmei Zhao; Hongmei Zhang; Gang Hu; Hang Zeng; Song Li
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Due to the influence of data collection methods and external complex factors, missing traffic data is inevitable. However, complete traffic information is a necessary input for route planning and forecasting tasks. To reduce the impact of missing data problems, this paper uses the low-rank tensor completion framework based on T-SVD to complete the missing spatio-temporal traffic data, the aim is to recover a low-rank tensor from a tensor with partial observation terms, and the WLRTC-P model is proposed. We use the idea of direction weighting to solve the dependence of the original model on the data input direction, extract each direction correlation information of the tensor spatio-temporal traffic data, and use the p-shrinkage norm to replace the tensor average rank minimization problem, and the study shows that the p-shrinkage norm is tighter than the tensor nuclear norm and, finally, uses the alternating direction method of multipliers to solve this model. Experiments on two publicly available spatio-temporal traffic datasets verified the conjecture of data input direction’s influence on the completion accuracy, and compared with the existing classical model methods, WLRTC-P has high precision and generalization ability.


---

## Body

## 1. Introduction

With the rapid development of Internet technology and traffic informatization, the scale of traffic data [1] is getting larger and larger, and intelligent transportation system (ITS) has been widely used in the management and control of roads as an important means to improve the efficiency of the transportation system and achieve sustainable development [2, 3] of transportation. In recent decades, the development and research on ITS, which integrates road, vehicle, traffic participants, and high-tech factors, have been flourishing and gaining great progress, which is the future direction of transportation system development.As the core of ITS, the traffic information acquisition system collects road traffic flow data from various sensors to form a complete urban traffic data set. At present, there are four main ways of collecting traffic flow data on urban expresswaysinduction coil, remote microwave, floating vehicle detection, and video detection. However, due to the acquisition equipment is high-precision equipment, vulnerable to the influence of external environments such as weather, temperature, sensors, communication equipment failure, and other factors, resulting in the final acquisition of traffic flow data abnormal situations, therefore, the problem of data missing [4, 5] in sensors and traffic information system is very common. The missing problem seriously affects a series of activities in intelligent transportation, such as monitoring traffic situations, predicting traffic flow [6–8], and deploying traffic planning.Therefore, the restoration of missing data on traffic flow is necessary, and how to recover the missing data quickly has become a hot issue in the field of intelligent transportation. For the existing common missing data phenomenon in the traffic data collection process, domestic and foreign scholars have proposed a series of solutions to recover the missing data and achieved certain results. However, in the field of traffic data recovery, there are still many open problems worth studying and exploring, for example, with the rapid development of urbanization, the urban traffic situation becomes more complex, and there are more and more factors affecting the traffic condition. Finding and using the most relevant features for data modeling among many factors and further recovering the missing data are a hot research topic in the field of traffic data recovery at present. For this reason, effective completion of missing traffic data is of great research significance at theory and practice levels.


## 2. Related Work

Spatio-temporal traffic data recovery methods can be divided into three categoriesprediction-based, statistical learning, and other machine learning methods. Among them, prediction-based methods will generate large data processing errors when repairing the missing data set, while most other machine learning-based methods have low interpretability; therefore, the existing research mainly focuses on using statistical learning methods for traffic data restoration. From the perspective of time development, statistical learning-based methods have mainly gone through three stagesvector-based, matrix-based, and tensor-based, and the research shows that the repair accuracy based on matrix and tensor is higher.Among the matrix-based methods, Bayesian principal component analysis (BPCA) and probabilistic principal component analysis (PPCA) proposed by Qu et al. [9, 10] in the early stage are the most representative. Subsequently, Tan et al. [11, 12] first introduced the tensor mode to the traffic data completion and proposed the TDI model based on Tucker decomposition, showing that the tensor model could cover more spatio-temporal traffic information, who studied the mode correlation information of tensor traffic data and tried to recommend the most suitable tensor model for the traffic flow data completion. After that, the tensor-based completion methods began to be applied in the traffic data completion field. Then, Liu et al. [13, 14] proposed three classical low-rank tensor completion models, SiLRTC (simple low-rank tensor completion model), FaLRTC (fast low-rank tensor completion model), and HaLRTC (high accuracy low-rank tensor completion model), respectively, which gave a generalized computational method for the completion problem and highlighted the accuracy advantage of the model and the characteristics of high-dimensional data in comparison with matrix completion. Tan et al. [15] further combined Tucker with Wopt to propose the Tucker–Wopt model, and it was shown to successfully reconstruct tensors with noise and up to 95% missing data; Chen et al. [16, 17] then extended the Bayesian decomposition of the matrix to the tensor model and combined it with the CP decomposition of the tensor proposed the BGCP model, and after adding the time regularization parameter proposed the BTTF mode. Chen et al. [18–20] have since resolved and optimized the HaLRTC model and proposed LRTC-TNN (low-rank tensor completion model based on truncated nuclear norm), LSTC-Tubal (low-tubal-rank smoothing tensor completion), and LATC (low-rank autoregressive tensor completion) model, and experiments on real spatio-temporal traffic datasets have obtained good results. Zemin et al. [21] proposed the second t-SVD-based LRTC framework using the concept of tensor average rank, which transfers the tensor pattern to the Fourier domain for solving and did not need to expand the tensor to break the spatio-temporal correlation information of multidimensional channels. Song et al. [22] proposed a weighted residual model of the tensor (T-WTNNR) on this basis, which was solved by the gradient descent method and obtained good results in the image domain. Cl et al. [23] proposed the LRTC-p model based p-shrinkage norm, p-shrinkage norm replaced the original tensor nuclear norm, and proved the way of p-shrinkage norm is more compact and can achieve better results. Kong et al. [24] proposed the Schatten-p norm and used the Schatten-p norm to replace the multiple TNN or nuclear norm of tensors, showing strong applicability in the two LRTC frameworks.However, the LRTC framework based on multiple TNNs requires the tensor to be expanded into multiple weighting matrices for the completion computation when recovering the tensor, which may break the modal correlation of the tensors and lose key information; the t-SVD-based LRTC framework operates only on the frontal slicing of the input data when processing the data and has no prior in spatio-temporal traffic data restoration. To address these problems, we propose a weighted optimization model based on the p-shrinkage norm on the second LRTC framework, which solved the above problems and applied the completion of traffic data.The main contributions of this paper are twofold1) we proposed to apply the t-SVD-based LRTC framework for modeling, integrating the weighting idea of the first framework, extracting the correlation information in each direction of the tensor, reducing the influence of dependence due to the data input direction, applying the p-shrinkage norm to preserve the strong correlation information inside, defining the model as WLRTC-P, and finally, using the alternating direction multiplier method to solve it.2) Two publicly available spatio-temporal traffic datasets are selected for experimental analysis, and the missing rate is set to 20%–80% under different missing scenarios. Experiment1 verified the influence problem of data input direction on the completion accuracy; Experiment 2 compared with the existing high-precision completion models and highlighted the accuracy advantage and generalizability capability of the WLRTC-P model.The rest of the paper is organized as followsin section3, we briefly introduced the tensor basis and the classical tensor completion model. In section 4, the tensor p-shrinkage norm and the external weighted way are introduced, and the WLRTC-P model is proposed and solved. In section 5, the problem of the effect of data input direction is verified, experimented on several publicly available datasets, and compared with several state-of-the-art baseline models. In section 6, this study is summarized.


## 3. Tensor Notation and Tensor Completion Model

### 3.1. Tensor Basis

In this section, some basic notations are first introduced, and then, some necessary definitions [25] are briefly provided for use.Use bold Euler letters to denote the tensors, e.g.,A; use bold capital letters to indicate the matrices, e.g., A; use bold lowercase letters to indicate vectors, e.g., a; scalars are represented by lowercase letters, e.g., a; a unit matrix of size n×n is denoted by In, and the real and complex domains are denoted by ℝ and ℂ, respectively.For a 3-way tensorA∈ℂn1×n2×n3, the elements are represented in the same way as matrices, such as Aijk denotes the i,j,k-th term of the tensor A. For the slice fibers of the tensor, Ak,:,:, A:,k,:, and A:,:,k denote the i-th horizontal, lateral, and frontal slices, respectively, and the slices of the tensor are in the form of a matrix. Besides, the frontal slices A:,:,k is denoted by Ai, and conjA is noted as the conjugate complex of A∈ℂn1×n2×n3, which is meant to take the complex conjugate of all elements in A.The inner product of tensorsA∈ℂn1×n2×n3 and ℬ∈ℂn1×n2×n3 is defined as A,ℬ=∑i=1n3Ai,Bi=A,ℬ=∑i1=1n1∑i2=1n2∑i3=1n3ai1,i2,i3bi1,i2,i3; the trace of a tensor A∈ℂn1×n2×n3 is defined as trA=∑i=1n3Ai; the Frobenius norm is defined as AF=∑i1=1n1∑i2=1n2∑i3=1n3ai1,i2,i32; and the expansion and collapse mapping operators are defined, respectively, as(1)unfoldA=A1A2⋮An3foldunfoldA=A.Theunfold is an expansion operator that projects the tensor A to the matrix, and the size after the expansion is n1n3×n2.ForA∈ℂn1×n2×n3, A¯=fftA,,2 represents the result of the discrete Fourier transformation (DFT) of A along the third dimension, and the tensor A¯ of the Fourier domain can be recovered to the real domain A by the inverse discrete Fourier transform (IDFT), A=fftA¯,2. For A¯∈ℂn1×n2×n3, we define A¯ as a block diagonal matrix in which each frontal slice A¯i of A¯ lies on the diagonal, denoted as(2)A¯=bdiagA¯=A¯1A¯2OA¯n3.The block circulant matrix of tensorA is defined bcircA and is noted as(3)bcircA=A1An3…A2A2A1…A3⋮⋮⋱⋮An3An3−1…A1.An important property here is that the block circulant matrix can be block diagonalized in the Fourier domain as(4)A¯=Fn3⊗In1⋅bcircA⋅Fn3−1⊗In2,where Fn3 is denoted as a discrete Fourier transform matrix of size n3×n3, In1 denotes a matrix of size ni×ni, and ⊗ denotes the Kronecker product.Using the properties of the discrete Fourier transform, forA¯∈ℂn1×n2×n3, the frontal slice is taken by default, so the discrete Fourier transform is carried out along the third dimension, and then, the conjugate symmetry of the real-valued signal transformed to the Fourier domain is(5)A¯1∈ℝn1×n2,conjA¯i=A¯n3−i+2,i=2,…,n3+12.The nature of conjugate symmetry of real-valued signals in the frequency domain is elaborated using the properties of the discrete Fourier transform, which will help to avoid redundant calculations.Definition 1.
(t-product [26])The t-product between A∈ℂn1×d×n3 and ℬ∈ℂd×n2×n3 is defined as a 3-way tensor of size n1×n2×n3.(6)C=A∗ℬ=foldbcircA⋅unfoldℬ.
T-product is identical to matrix multiplication in form except for cyclic convolution and expansion operations. It is worth noting that t-product reduces to matrix multiplication whenn3=1.Definition 2.
(Conjugate transpose and orthogonality of a tensor)Given a tensorA∈ℂn1×n2×n3, its conjugate transpose is A∗∈ℂn2×n1×n3, which is obtained by conjugate transposing each slice of A and then reversing the order of the transposed frontal slices 2 through n3. Given a tensor A∈ℂn2×n1×n3, if it satisfies the orthogonality condition, then there is A∗∗A=A∗A∗=ℐ and ∗ is the t-product.Definition 3.
(Unit tensor and F-diagonal tensor)Given a unit tensor of sizeA∈ℝI×I×n, the first frontal slice of the unit tensor is a unit matrix of size I×I, and all other forward slices are zero. Given a tensor A∈ℂn1×n2×n3, if each of its frontal slices is a diagonal matrix, then this tensor is called an F-diagonal tensor.Theorem 1.
(t–SVD[26])For a 3-way tensor A∈ℂn1×n2×n3, there exists tensor singular value decomposition:(7)A=U∗S∗V⊤,where U∈ℂn1×n2×n3 and V∈ℂn1×n2×n3 are two orthogonal tensors, S∈ℝn1×n2×n3 is an F-diagonal tensor. According to (4), t-SVD can be performed efficiently based on matrix singular value decomposition in the Fourier domain. In the Fourier domain, each frontal slice can be independently singular value decomposed, A¯=U¯∗S¯∗V¯. Applying this property, we can obtain Us¯,S¯,V¯ in the Fourier domain and, finally, use the inverse discrete Fourier transform to obtain U¯S¯V¯. The tensor singular value decomposition is shown in Figure 1.
The process of decomposing a third-order tensor through a tensor t-product is called t-SVD, and the algorithm is summarized as follows:
Algorithm 1: t-SVD.
InputA∈ℝn1×n2×n3(1)
A¯=fftA,,3(2)
fori =  1,…,n3+1/2 do(3)
U¯:,:,i,S¯:,:,i,V¯:,:,i=SVDA¯:,:,i(4)
end for(5)
fori =  1,…,n3+1/2 do(6)
U¯i=conjU¯n3−i+2,S¯i=conjS¯n3−i+2,V¯i=conjV¯n3−i+2(7)
end for(8)
U=ifftU,,3S=ifftS,,3V=ifftV,,3OutputU∈ℝn1×n2×n3,S∈ℝn1×n2×n3,V∈ℝn1×n2×n3Figure 1
Tensor singular value decomposition.Definition 4.
(Tensor tubal rank and average rank [21])Given a tensor A∈ℝn1×n2×n3, assuming a t-SVD of U∗S∗V⊤ for tensor A, the tensor tubal rank of A is noted as ranktA, define the number of nonzero terms on S; i.e.,(8)ranktA=#i,Si,i,≠0=#i,Si,i,1≠0.
The concept of the average rank of a tensorgiven a tensorA∈ℝn1×n2×n3, the average rank is noted as rankaA, defined as follows :(9)rankaA=1n3bcircA.Definition 5.
(Spectral norm and nuclear norm)for a given tensorA∈ℝn1×n2×n3, its tensor spectral norm is defined as A=bcircA.
For the tensor nuclear norm, assuming that the t-SVD of the tensor isU∗S∗V⊤ for tensor A, the tensor nuclear norm is defined as the sum of the tensor singular values, as(10)A∗=∑i=1rSi,i,1=∑i=1rσiA.
But the tensor nuclear norm is only determined by the first frontal sliceSi,i,1, which is different from the nuclear norm of the matrix. Note that ∥A∥∗ is induced by the t-product, it is the dual norm of the tensor spectral norm. A∗ is the convex envelope of the tensor average rank, and using (3), the following relation (11) can be obtained:(11)A∗=1n3bcircA∗=1n3bdiagA¯∗.Theorem 2.
(Tensor singular value thresholding [27])given a tensor A∈ℂn1×n2×n3, its t-SVD is A=U∗S∗V⊤.
Given parametersτ and A, the tensor singular value thresholding is defined as (12)(12)DτA=U∗Sτ∗V⊤,Sτ=ifftmaxS^−τ,0,3.
For anyτ > 0 and A∈ℂn1×n2×n3, tensor singular value thresholding is related to the nuclear norm as follows:(13)DτA=nargminxɛRn1×n2×n3τpxp∗+pA−xpF2.


### 3.2. Tensor Complementary Model

The tensor completion framework based on t-SVD was first proposed by Zemin et al. and introduced the problem of minimizing the average rank of tensor, which can be regarded as a high-dimensional extension of matrix completion. Similar to the classical HaLRTC model problem, it aims to accurately fill the positions of unobserved elements with the correlation information inside the data. Experiments on image and video processing show that the algorithm based on low-rank tensor completion is numerically superior to matrix completion and other tensor completion methods. However, the LRTC framework based on t-SVD does not have a priori in the problem of spatio-temporal traffic data, therefore, for partially observed spatio-temporal traffic data, constructed as tensor structureℳ∈ℝn1×n2×n3, and applied the concept of tensor average rank to model and analysis, as follows:(14)minX:rankaX,s.t:pΩM=pΩX,where X is a tensor to be recovered; ℳ is an observation tensor with missing elements; minX:rankaX is the problem of finding the average rank minimization of the completion tensor. The constraints require that while solving the X rank minimization, the position elements of the original tensor ℳ with observations are equal to the elements of the relative positions of the completion tensor X, and Ω is the set representation of the observation elements of the tensor.(15)pΩxmntxmnt,ifm,n,t∈Ω,0,otherwise.For an arbitrary tensorX, where PΩX denotes the set of element positions with observations, and PΩ¯ denotes the set of element positions without observations, and the two operators PΩ and PΩ¯ have a complementary relationship with each other, PΩ¯X+PΩ¯X=X.The rank minimization problem in (14) for tensor completion is NP-hard and computationally intractable. According to Definition 5, the convex envelope of X average rank is its nuclear norm X∗, therefore, using the tensor nuclear norm to replace the tensor average rank minimization problem, and as the convex envelope of tensor average rank minimization:(16)minX:X∗,s.t:PΩM=PΩX,where min:X∗X means minimizing the nuclear norm of the tensor X, thus, transforming the original problem into a solvable form.


## 3.1. Tensor Basis

In this section, some basic notations are first introduced, and then, some necessary definitions [25] are briefly provided for use.Use bold Euler letters to denote the tensors, e.g.,A; use bold capital letters to indicate the matrices, e.g., A; use bold lowercase letters to indicate vectors, e.g., a; scalars are represented by lowercase letters, e.g., a; a unit matrix of size n×n is denoted by In, and the real and complex domains are denoted by ℝ and ℂ, respectively.For a 3-way tensorA∈ℂn1×n2×n3, the elements are represented in the same way as matrices, such as Aijk denotes the i,j,k-th term of the tensor A. For the slice fibers of the tensor, Ak,:,:, A:,k,:, and A:,:,k denote the i-th horizontal, lateral, and frontal slices, respectively, and the slices of the tensor are in the form of a matrix. Besides, the frontal slices A:,:,k is denoted by Ai, and conjA is noted as the conjugate complex of A∈ℂn1×n2×n3, which is meant to take the complex conjugate of all elements in A.The inner product of tensorsA∈ℂn1×n2×n3 and ℬ∈ℂn1×n2×n3 is defined as A,ℬ=∑i=1n3Ai,Bi=A,ℬ=∑i1=1n1∑i2=1n2∑i3=1n3ai1,i2,i3bi1,i2,i3; the trace of a tensor A∈ℂn1×n2×n3 is defined as trA=∑i=1n3Ai; the Frobenius norm is defined as AF=∑i1=1n1∑i2=1n2∑i3=1n3ai1,i2,i32; and the expansion and collapse mapping operators are defined, respectively, as(1)unfoldA=A1A2⋮An3foldunfoldA=A.Theunfold is an expansion operator that projects the tensor A to the matrix, and the size after the expansion is n1n3×n2.ForA∈ℂn1×n2×n3, A¯=fftA,,2 represents the result of the discrete Fourier transformation (DFT) of A along the third dimension, and the tensor A¯ of the Fourier domain can be recovered to the real domain A by the inverse discrete Fourier transform (IDFT), A=fftA¯,2. For A¯∈ℂn1×n2×n3, we define A¯ as a block diagonal matrix in which each frontal slice A¯i of A¯ lies on the diagonal, denoted as(2)A¯=bdiagA¯=A¯1A¯2OA¯n3.The block circulant matrix of tensorA is defined bcircA and is noted as(3)bcircA=A1An3…A2A2A1…A3⋮⋮⋱⋮An3An3−1…A1.An important property here is that the block circulant matrix can be block diagonalized in the Fourier domain as(4)A¯=Fn3⊗In1⋅bcircA⋅Fn3−1⊗In2,where Fn3 is denoted as a discrete Fourier transform matrix of size n3×n3, In1 denotes a matrix of size ni×ni, and ⊗ denotes the Kronecker product.Using the properties of the discrete Fourier transform, forA¯∈ℂn1×n2×n3, the frontal slice is taken by default, so the discrete Fourier transform is carried out along the third dimension, and then, the conjugate symmetry of the real-valued signal transformed to the Fourier domain is(5)A¯1∈ℝn1×n2,conjA¯i=A¯n3−i+2,i=2,…,n3+12.The nature of conjugate symmetry of real-valued signals in the frequency domain is elaborated using the properties of the discrete Fourier transform, which will help to avoid redundant calculations.Definition 1.
(t-product [26])The t-product between A∈ℂn1×d×n3 and ℬ∈ℂd×n2×n3 is defined as a 3-way tensor of size n1×n2×n3.(6)C=A∗ℬ=foldbcircA⋅unfoldℬ.
T-product is identical to matrix multiplication in form except for cyclic convolution and expansion operations. It is worth noting that t-product reduces to matrix multiplication whenn3=1.Definition 2.
(Conjugate transpose and orthogonality of a tensor)Given a tensorA∈ℂn1×n2×n3, its conjugate transpose is A∗∈ℂn2×n1×n3, which is obtained by conjugate transposing each slice of A and then reversing the order of the transposed frontal slices 2 through n3. Given a tensor A∈ℂn2×n1×n3, if it satisfies the orthogonality condition, then there is A∗∗A=A∗A∗=ℐ and ∗ is the t-product.Definition 3.
(Unit tensor and F-diagonal tensor)Given a unit tensor of sizeA∈ℝI×I×n, the first frontal slice of the unit tensor is a unit matrix of size I×I, and all other forward slices are zero. Given a tensor A∈ℂn1×n2×n3, if each of its frontal slices is a diagonal matrix, then this tensor is called an F-diagonal tensor.Theorem 1.
(t–SVD[26])For a 3-way tensor A∈ℂn1×n2×n3, there exists tensor singular value decomposition:(7)A=U∗S∗V⊤,where U∈ℂn1×n2×n3 and V∈ℂn1×n2×n3 are two orthogonal tensors, S∈ℝn1×n2×n3 is an F-diagonal tensor. According to (4), t-SVD can be performed efficiently based on matrix singular value decomposition in the Fourier domain. In the Fourier domain, each frontal slice can be independently singular value decomposed, A¯=U¯∗S¯∗V¯. Applying this property, we can obtain Us¯,S¯,V¯ in the Fourier domain and, finally, use the inverse discrete Fourier transform to obtain U¯S¯V¯. The tensor singular value decomposition is shown in Figure 1.
The process of decomposing a third-order tensor through a tensor t-product is called t-SVD, and the algorithm is summarized as follows:
Algorithm 1: t-SVD.
InputA∈ℝn1×n2×n3(1)
A¯=fftA,,3(2)
fori =  1,…,n3+1/2 do(3)
U¯:,:,i,S¯:,:,i,V¯:,:,i=SVDA¯:,:,i(4)
end for(5)
fori =  1,…,n3+1/2 do(6)
U¯i=conjU¯n3−i+2,S¯i=conjS¯n3−i+2,V¯i=conjV¯n3−i+2(7)
end for(8)
U=ifftU,,3S=ifftS,,3V=ifftV,,3OutputU∈ℝn1×n2×n3,S∈ℝn1×n2×n3,V∈ℝn1×n2×n3Figure 1
Tensor singular value decomposition.Definition 4.
(Tensor tubal rank and average rank [21])Given a tensor A∈ℝn1×n2×n3, assuming a t-SVD of U∗S∗V⊤ for tensor A, the tensor tubal rank of A is noted as ranktA, define the number of nonzero terms on S; i.e.,(8)ranktA=#i,Si,i,≠0=#i,Si,i,1≠0.
The concept of the average rank of a tensorgiven a tensorA∈ℝn1×n2×n3, the average rank is noted as rankaA, defined as follows :(9)rankaA=1n3bcircA.Definition 5.
(Spectral norm and nuclear norm)for a given tensorA∈ℝn1×n2×n3, its tensor spectral norm is defined as A=bcircA.
For the tensor nuclear norm, assuming that the t-SVD of the tensor isU∗S∗V⊤ for tensor A, the tensor nuclear norm is defined as the sum of the tensor singular values, as(10)A∗=∑i=1rSi,i,1=∑i=1rσiA.
But the tensor nuclear norm is only determined by the first frontal sliceSi,i,1, which is different from the nuclear norm of the matrix. Note that ∥A∥∗ is induced by the t-product, it is the dual norm of the tensor spectral norm. A∗ is the convex envelope of the tensor average rank, and using (3), the following relation (11) can be obtained:(11)A∗=1n3bcircA∗=1n3bdiagA¯∗.Theorem 2.
(Tensor singular value thresholding [27])given a tensor A∈ℂn1×n2×n3, its t-SVD is A=U∗S∗V⊤.
Given parametersτ and A, the tensor singular value thresholding is defined as (12)(12)DτA=U∗Sτ∗V⊤,Sτ=ifftmaxS^−τ,0,3.
For anyτ > 0 and A∈ℂn1×n2×n3, tensor singular value thresholding is related to the nuclear norm as follows:(13)DτA=nargminxɛRn1×n2×n3τpxp∗+pA−xpF2.


## 3.2. Tensor Complementary Model

The tensor completion framework based on t-SVD was first proposed by Zemin et al. and introduced the problem of minimizing the average rank of tensor, which can be regarded as a high-dimensional extension of matrix completion. Similar to the classical HaLRTC model problem, it aims to accurately fill the positions of unobserved elements with the correlation information inside the data. Experiments on image and video processing show that the algorithm based on low-rank tensor completion is numerically superior to matrix completion and other tensor completion methods. However, the LRTC framework based on t-SVD does not have a priori in the problem of spatio-temporal traffic data, therefore, for partially observed spatio-temporal traffic data, constructed as tensor structureℳ∈ℝn1×n2×n3, and applied the concept of tensor average rank to model and analysis, as follows:(14)minX:rankaX,s.t:pΩM=pΩX,where X is a tensor to be recovered; ℳ is an observation tensor with missing elements; minX:rankaX is the problem of finding the average rank minimization of the completion tensor. The constraints require that while solving the X rank minimization, the position elements of the original tensor ℳ with observations are equal to the elements of the relative positions of the completion tensor X, and Ω is the set representation of the observation elements of the tensor.(15)pΩxmntxmnt,ifm,n,t∈Ω,0,otherwise.For an arbitrary tensorX, where PΩX denotes the set of element positions with observations, and PΩ¯ denotes the set of element positions without observations, and the two operators PΩ and PΩ¯ have a complementary relationship with each other, PΩ¯X+PΩ¯X=X.The rank minimization problem in (14) for tensor completion is NP-hard and computationally intractable. According to Definition 5, the convex envelope of X average rank is its nuclear norm X∗, therefore, using the tensor nuclear norm to replace the tensor average rank minimization problem, and as the convex envelope of tensor average rank minimization:(16)minX:X∗,s.t:PΩM=PΩX,where min:X∗X means minimizing the nuclear norm of the tensor X, thus, transforming the original problem into a solvable form.


## 4. WLRTC-P Model

### 4.1. P-Shrinkage Norm

First, describe the definition of the matrix p-shrinkage norm [28]given a matrix X∈ℝm×n, Gp as a map of ℝm×n⟶ℝ, the proximal operator of the p-shrinkage norm is Gpμ, and the p-shrinkage norm is defined as(17)GpμX=GpμσX=∑i=1ngpμσiX,where GpμσX means taking the p-shrinkage norm of the matrix, σiX means taking the i-th singular value of X, and the singular value size ordering of matrix X is σ1⩾σ2⩾⋯⩾σminm,n⩾0.Extending the p-shrinkage norm to the tensor model, for tensorX∈ℝn1×n2×n3, the p-shrinkage norm under the tensor singular value decomposition can be defined as follows:(18)Gpμx=1n3∑k=1n3GpμσX¯k=1n3∑k=1n3∑i=1minn1,n2gpμσX¯k.Corollary 1.
Optimal solution for any−∞<p<1, X∈ℝn1×n2×n3, and α, ρ>0 problem:(19)Spα/ρx=αGpμx+ρ2X−ℳF2.
Given the generalized singular value thresholding [29], defined the threshold operator for the p-shrinkage norm as Spμ, we have(20)Spα/ρx=USpα/ρΣijkV⊤,Spα/ρX is a threshold decomposition for each slice of the tensor, and the whole process can be done in the Fourier domain, so there have(21)Spα/ρΣijk=DiagσiX¯−αρσiX¯p−1,0+,Σijk denotes the tensor singular value of X, and for any a∈ℝ, the additive function .+ is defined as a+=maxa,0. P-shrinkage norm is meant to preserve the characteristics of large singular values, which shrinks less for larger singular values and penalizes excessively for smaller ones, and can be seen as another way of threshold shrinkage.
Based on the inference of [23], when −∞<p<1, GpμX is a nonconvex envelope of the average rank of the tensor in the unit sphere of the spectral norm, and it is tighter than the tensor nuclear norm, and we have pxp∗≤GpμX≤rankaX.


### 4.2. WLRTC-P Model

The original low-rank tensor completion framework is based on t-SVD as in (14), when there is data input, the nuclear norm of the tensor is only related to the first frontal slice of the tensor, and the lateral and horizontal slices are not processed; however, the frontal slice of the tensor contains data information in only two dimensions, which is contrary to our idea of using tensor modeling. Thus, to avoid the dependence of the model on the data input direction, and to solve the problem that one-way slicing does not fully utilize the correlation information within the data, inspired by the use of multiple TNN instead of the tensor rank minimization problem by LIU et al., take the weighted nuclear norm of the three-way slice to instead of the tensor average rank minimization problem, and as much as possible to reduce the influence of data input direction. Define the three tensors [30] size as shown in Figure 2, X1∈ℝn1×n2×n3, X2∈ℝn2×n3×n1, X3∈ℝn3×n1×n2, using three tensors weighting to replace the original tensor, weighting is desired to retain each direction data channel’s correlation information to improve the completion accuracy. So, we weighted the original tensor completion model and performed parameter optimization to recreate the model. Therefore, for the original problem (14), the model can be rewritten as follows:(22)minXi:∑i=13αixi∗,s.t:pΩxi=PΩxi.Figure 2
Three ways to slice a tensor (frontal slicing, lateral slicing, and horizontal slicing).The sum of∑i=13αi weights is 1. For the external three-way weighting, the tensor nuclear norm is taken in three directions, and the nuclear norm in each direction is given by(23)∑i=13αiXi∗=∑i=13αimibcircXi∗=∑i=13αimiX¯i∗,where m1=n3,m2=n1,m3=n2, X¯i are the corresponding block-diagonal matrices obtained from each direction tensor X by (2).Adding the idea of the p-shrinkage norm, the model is rebuilt and defined as WLRTC-P. To reduce the dependence effect, introduced three third-order auxiliary tensor variablesℳii=123 and an additional set of constraints ℳi=Xi, and transformed the above problem into a tractable form as follows:(24)minX,D,Z,W:αiGpμXi,s.tPΩMi=PΩYi,Mi=Xii=1,23.To solve this optimization problem, a straightforward and widely used approach is the Alternating Direction Method of Multipliers (ADMM [31]) framework. The ADMM framework is mainly applied to solve large-scale problems and to solve optimization problems with multiple nonsmooth terms in the objective, and the general idea is to decompose the large global problem into multiple smaller, easier-to-solve local subproblems by decomposing the coordination process, and to obtain the optimal solution of the global problem by coordinating the solutions of the subproblems. Defined the enhanced Lagrangian function, as follows:(25)LρXi,ℳi,Ti,ρ=αiGpXi+ρ2Xi−ℳiF2+Ti,Xi−ℳi,where Z∈ℝn1×n2×n3 is the Lagrange multiplier, ρ is the penalty parameter after adding the constraint, and the auxiliary variable Tili=1,2,3 is used for the dual update under ADMM. According to the solution of the ADMM framework, decompose xiMilTilp into local subproblems, respectively, and the optimal solution of the problem is approximated by alternating updates, fixing the other position elements in the process, and the order of iterative updates is as follows:(1)
First updateXil+1, fixing the parameters of ℳil and Til.(27)Xil+1=argminXℒXil,Mil,Til,ρ=argminXαiGpμXil+ρ2Xil−MilF2+Xil,Til=argminXαiGpμXil+ρ2Xil−Mil−1ρTilF2=argminXSpαi/ρkMil−1ρTil.Spαi/ρkℳil−1/ρTil is the tensor singular value thresholding decomposition, where μ=αi/ρk.(2)
Updateℳil+1, fixing the parameters of Xil+1 and Ti1.(28)Mil+1=argminMLMil,Xil+1,Til,ρ=argminM∑i=13ρ2Xil+1−MilF2+Mil,Til=argminM∑i=13ρ2MilF2−PXil+1,Mil+Mil,Til=∑i=13Xil+1+1PTil.Adding constraintPΩMi=PΩYi, we have(29)Mil+1=∑i=13Xil+1+1ρTilΩ¯+PΩ¯Mil.(3)
Finally solve forTil+1:(30)Til+1=Til+ρXil+1−ℳil+1.Set the number of iterations, update each iteration according to the above steps, and seek the optimal solution of the model. Thus, the steps of the WLRTC-P algorithm are as follows:Algorithm 2: WLRTC-P.
InputXi1, Mi1, PΩMil=PΩXil, ℳil=Til=0, αi, ρ, k(1)
fori=1 to k(2)
Update：Xil+1=Spαi/ρkℳil−1/ρTil(3)
Update：Mil+1=∑i=13Xil+1+1/pTilΩ¯+PΩ¯Mil(4)
Update：Til+1:=Til+ρXil+1−Mil+1(5)
if：Ml+1−MlF2/PΩMF2<EaOutputℳIn the whole solution process, the influence of noise is not considered, and the WRTC-P model only considers the reconstruction accuracy of the given data. Similar to most statistical learning models, it improves the recovery accuracy of traffic flow data by iterating the existing low-rank problem of observed data, so the problem of data crossing will not occur in the migration process, and it is suitable for most low-rank traffic flow data restoration problems.


## 4.1. P-Shrinkage Norm

First, describe the definition of the matrix p-shrinkage norm [28]given a matrix X∈ℝm×n, Gp as a map of ℝm×n⟶ℝ, the proximal operator of the p-shrinkage norm is Gpμ, and the p-shrinkage norm is defined as(17)GpμX=GpμσX=∑i=1ngpμσiX,where GpμσX means taking the p-shrinkage norm of the matrix, σiX means taking the i-th singular value of X, and the singular value size ordering of matrix X is σ1⩾σ2⩾⋯⩾σminm,n⩾0.Extending the p-shrinkage norm to the tensor model, for tensorX∈ℝn1×n2×n3, the p-shrinkage norm under the tensor singular value decomposition can be defined as follows:(18)Gpμx=1n3∑k=1n3GpμσX¯k=1n3∑k=1n3∑i=1minn1,n2gpμσX¯k.Corollary 1.
Optimal solution for any−∞<p<1, X∈ℝn1×n2×n3, and α, ρ>0 problem:(19)Spα/ρx=αGpμx+ρ2X−ℳF2.
Given the generalized singular value thresholding [29], defined the threshold operator for the p-shrinkage norm as Spμ, we have(20)Spα/ρx=USpα/ρΣijkV⊤,Spα/ρX is a threshold decomposition for each slice of the tensor, and the whole process can be done in the Fourier domain, so there have(21)Spα/ρΣijk=DiagσiX¯−αρσiX¯p−1,0+,Σijk denotes the tensor singular value of X, and for any a∈ℝ, the additive function .+ is defined as a+=maxa,0. P-shrinkage norm is meant to preserve the characteristics of large singular values, which shrinks less for larger singular values and penalizes excessively for smaller ones, and can be seen as another way of threshold shrinkage.
Based on the inference of [23], when −∞<p<1, GpμX is a nonconvex envelope of the average rank of the tensor in the unit sphere of the spectral norm, and it is tighter than the tensor nuclear norm, and we have pxp∗≤GpμX≤rankaX.


## 4.2. WLRTC-P Model

The original low-rank tensor completion framework is based on t-SVD as in (14), when there is data input, the nuclear norm of the tensor is only related to the first frontal slice of the tensor, and the lateral and horizontal slices are not processed; however, the frontal slice of the tensor contains data information in only two dimensions, which is contrary to our idea of using tensor modeling. Thus, to avoid the dependence of the model on the data input direction, and to solve the problem that one-way slicing does not fully utilize the correlation information within the data, inspired by the use of multiple TNN instead of the tensor rank minimization problem by LIU et al., take the weighted nuclear norm of the three-way slice to instead of the tensor average rank minimization problem, and as much as possible to reduce the influence of data input direction. Define the three tensors [30] size as shown in Figure 2, X1∈ℝn1×n2×n3, X2∈ℝn2×n3×n1, X3∈ℝn3×n1×n2, using three tensors weighting to replace the original tensor, weighting is desired to retain each direction data channel’s correlation information to improve the completion accuracy. So, we weighted the original tensor completion model and performed parameter optimization to recreate the model. Therefore, for the original problem (14), the model can be rewritten as follows:(22)minXi:∑i=13αixi∗,s.t:pΩxi=PΩxi.Figure 2
Three ways to slice a tensor (frontal slicing, lateral slicing, and horizontal slicing).The sum of∑i=13αi weights is 1. For the external three-way weighting, the tensor nuclear norm is taken in three directions, and the nuclear norm in each direction is given by(23)∑i=13αiXi∗=∑i=13αimibcircXi∗=∑i=13αimiX¯i∗,where m1=n3,m2=n1,m3=n2, X¯i are the corresponding block-diagonal matrices obtained from each direction tensor X by (2).Adding the idea of the p-shrinkage norm, the model is rebuilt and defined as WLRTC-P. To reduce the dependence effect, introduced three third-order auxiliary tensor variablesℳii=123 and an additional set of constraints ℳi=Xi, and transformed the above problem into a tractable form as follows:(24)minX,D,Z,W:αiGpμXi,s.tPΩMi=PΩYi,Mi=Xii=1,23.To solve this optimization problem, a straightforward and widely used approach is the Alternating Direction Method of Multipliers (ADMM [31]) framework. The ADMM framework is mainly applied to solve large-scale problems and to solve optimization problems with multiple nonsmooth terms in the objective, and the general idea is to decompose the large global problem into multiple smaller, easier-to-solve local subproblems by decomposing the coordination process, and to obtain the optimal solution of the global problem by coordinating the solutions of the subproblems. Defined the enhanced Lagrangian function, as follows:(25)LρXi,ℳi,Ti,ρ=αiGpXi+ρ2Xi−ℳiF2+Ti,Xi−ℳi,where Z∈ℝn1×n2×n3 is the Lagrange multiplier, ρ is the penalty parameter after adding the constraint, and the auxiliary variable Tili=1,2,3 is used for the dual update under ADMM. According to the solution of the ADMM framework, decompose xiMilTilp into local subproblems, respectively, and the optimal solution of the problem is approximated by alternating updates, fixing the other position elements in the process, and the order of iterative updates is as follows:(1)
First updateXil+1, fixing the parameters of ℳil and Til.(27)Xil+1=argminXℒXil,Mil,Til,ρ=argminXαiGpμXil+ρ2Xil−MilF2+Xil,Til=argminXαiGpμXil+ρ2Xil−Mil−1ρTilF2=argminXSpαi/ρkMil−1ρTil.Spαi/ρkℳil−1/ρTil is the tensor singular value thresholding decomposition, where μ=αi/ρk.(2)
Updateℳil+1, fixing the parameters of Xil+1 and Ti1.(28)Mil+1=argminMLMil,Xil+1,Til,ρ=argminM∑i=13ρ2Xil+1−MilF2+Mil,Til=argminM∑i=13ρ2MilF2−PXil+1,Mil+Mil,Til=∑i=13Xil+1+1PTil.Adding constraintPΩMi=PΩYi, we have(29)Mil+1=∑i=13Xil+1+1ρTilΩ¯+PΩ¯Mil.(3)
Finally solve forTil+1:(30)Til+1=Til+ρXil+1−ℳil+1.Set the number of iterations, update each iteration according to the above steps, and seek the optimal solution of the model. Thus, the steps of the WLRTC-P algorithm are as follows:Algorithm 2: WLRTC-P.
InputXi1, Mi1, PΩMil=PΩXil, ℳil=Til=0, αi, ρ, k(1)
fori=1 to k(2)
Update：Xil+1=Spαi/ρkℳil−1/ρTil(3)
Update：Mil+1=∑i=13Xil+1+1/pTilΩ¯+PΩ¯Mil(4)
Update：Til+1:=Til+ρXil+1−Mil+1(5)
if：Ml+1−MlF2/PΩMF2<EaOutputℳIn the whole solution process, the influence of noise is not considered, and the WRTC-P model only considers the reconstruction accuracy of the given data. Similar to most statistical learning models, it improves the recovery accuracy of traffic flow data by iterating the existing low-rank problem of observed data, so the problem of data crossing will not occur in the migration process, and it is suitable for most low-rank traffic flow data restoration problems.


## 5. Experiment

In this section, experiments are conducted on the WLRTC-P model using publicly available spatio-temporal traffic datasets, to verify the conjecture of the directional dependence influence of the data input, and to evaluate the algorithm for its efficacy and generalizability capability.


### 5.1. Data Set Settings

Experiments are conducted using two publicly available datasets, both collected from real spatio-temporal traffic systems, and experiments on real datasets can verify the applicability of the algorithm in real scenarios. Both datasets have the same tensor pattern “location/sensor × day × time” except for the difference in tensor size.(1)
An urban traffic dataset for Guangzhou, China, for 61 days from August 1, 2016, to September 30, 2016, at 10-minute intervals, consists of 214 anonymous road segments (mainly composed of urban highways and arterial roads). According to the spatio-temporal attributes, establish a third-order tensor with the size of 214 × 61 × 144, where the dimensions of the tensor are road sections, dates, and time windows, respectively. For simplicity, this dataset is referred to as “G” for short.(2)
Seattle freeway traffic speed data collected by Cui et al. [32]. This dataset contains the freeway traffic speed in Seattle, USA, for the whole year of 2015 from 323 loop detectors with a resolution of 5 minutes. Select a subset of January (4 weeks from January 1st to January 28th) as experimental data, and build tensor data with a size of 323 × 28 × 288 according to the space-time attributes, where each dimension of the tensor is road sections, dates, and time windows, respectively. For simplicity, this dataset is referred to as “G” for short.Constructed the traffic data as a tensor pattern of links/sensors, the dates, and time window ensures that the tensor data are sufficiently low rank on each slice, to maximize the use of internal correlation information and thus improve the data restoration accuracy where the link dimension represents the spatial-level features of the traffic data, and the date and time windows represent the time-level features.


### 5.2. Baseline Model

Compared the WLRTC-P model with the following baseline model:HaLRTC [14]high-precision low-rank tensor completion model, the most classical low-rank tensor completion model (LRTC), is based on the minimization of the nuclear norm of the tensor and is solved using the alternating multiplier method (ADMM) framework.BGCP [16]Bayesian Gaussian CP decomposition, a full Bayesian tensor decomposition model, uses Markov chain, Monte Carlo, to learn the potential factor matrix (low-rank structure).CP_ALS [33]CP decomposition model, the classical tensor completion model, is solved by least squares and gradient descent methods, which is more accurate than the original solution method.TRMF [34]temporal regularized matrix decomposition is a temporal matrix decomposition model applied to a multiple autoregressive (AR) process to model the underlying temporal factors.For the selection of baseline model, considered two data construction methods. For matrix-based completion, the TRMF model was selected for comparison, and the data were processed as “location/sensor × time”; for tensor-based completion, the tensor data model was constructed as “location/sensor × day × time,” and the classical CP_ALS, HaLRTC, and BGCP models are selected for comparison. In addition, we have other considerations, the WLRTC-P model does not consider the influence of the directionality of the input tensor slices, and the original model is compared with the proposed model to show the importance of spatio-temporal information.


### 5.3. Experimental Setup

The datasets given in the experiments all tend to be complete, and for the model’s completion accuracy comparison, the true and completion values of the datasets need to be used to calculate the metric, so when conducting the experiments, the data need to be lost according to the set scenario with the missing rate, and then, the missing values are filled by the completion algorithm. For accuracy comparison, the accuracy of completion is judged by comparing the value of MAPE (%) and RMSE. The definitions of MAPE (%) and RMSE are given as (29)(31)MAPE=1n∑i=1nyi−y^iyi×100,RMSE=1n∑i=1nyi−y^i2.In real life, traffic data are missing in two waysrandom missing (RM) and nonrandom missing (NM). The RM case is aimless, mostly caused by unstable signals during communication and transmission, while the NM case data are destroyed in a related way, mostly caused by sensor damage. The two missing scenarios set allowed for a better evaluation of the different model’s performance and effectiveness.Taking two publicly available missing datasets as an example, randomly select a link, to show the data sparsity under 30% and 60% missing spatio-temporal traffic data in a day, as shown in Figure3.Figure 3
Data missing rate is 30% and 60% (Green means missing, and white means no missing issues). (a) G data set RM scenario, (b) S data set RM scenario, (c) G data set NM scenario, (d) S data set NM scenario.
(a)(b)(c)(d)


#### 5.3.1. Parameter Setting

Refer to the way of original LRTC framework was set up. For the iteration stepρ, set the initial value to ρ  = 1e8/1e9, set the maximum value ρmax of ρ as 1e5, and update ρ by ρ=min1.05×ρ,ρmax in each iteration update; for the weighting parameter αi, set the sum of parameters to 1, where α1=α2=α3=1/3; μ values to αi/ρk; for p values, set to 0.2; for the convergence condition, use Ml+1−MlF2/PΩMF2<ε to determine whether the algorithm converges; for the convergence accuracy limit E, set the value of E to 1e-4; for the selection of the number of iterations k, set k to 200 can guarantee the convergence of the model.


#### 5.3.2. Experimental Scheme

According to the actual collection of traffic data, we set two groups of experiments in the RM scenarios and the NM scenarios. In two missing scenarios, the missing rates ofG and S are set to 20%–80%, respectively. The first experiment is to judge the influence problem of data input direction (slices in different directions) on the completion accuracy by comparing three data sets with different input directions in the same scenario to verify our conjecture. The second is to conduct accuracy comparison experiments, to compare with the baseline model by experiments in the same scenario, and to analyze the performance of the model. In the comparison of the accuracy values in both sets of experiments, the smaller the values of RMSE and MAPE (%), the higher the completion accuracy of the algorithm, and the higher the accuracy advantage and model generalization ability.


### 5.4. Experiment

Experiment 1.
Influence of data input direction (slicing) on completion accuracy:
In this experiment, according to the above parameter settings, we use the one way of the original model (14) when the data input directions are X1∈ℝn1×n2×n3, X2∈ℝn2×n3×n1, and X3∈ℝn3×n1×n2, and the completion accuracy of traffic data is related to frontal slicing, lateral slicing, and horizontal slicing, and in the experiment, use direction 1, direction 2, and direction 3 to replace them, respectively. The following Figure 4 shows the comparison of slices in three different directions under different missing rates and missing scenarios.
As we can see, in most cases, the accuracy is highest when the data input direction is direction 1, and the completion accuracy is close to direction 1 when the data input direction is direction 2. In the NM scenario of the S data set, the accuracy of direction 2 is better than that of direction 1 in some cases of missing rate, and the accuracy of direction 3 is the worst in comparison. For this reason, we applied the correlation information of other directions in the subsequent WLRTC-P model by way of direction weighting to improve the accuracy. The above experiments verified our conjecture that the different input directions of the traffic data will lead to a significant influence on the accuracy of the subsequent completions.Figure 4
Effect of data input direction (slicing) on completion accuracy (The missing rates of G and S data sets were 20%–80%, respectively). (a) G data set RM scenario, (b) G data set NM scenario, (c) S data set RM scenario, (d) S data set NM scenario.
(a)(b)(c)(d)Experiment 2.
Comparison of completion accuracy under different missing scenarios and missing degrees:
The missing rates ofG and S are set to 20%–80%, respectively. Through experiments, comparing the completion accuracy of the algorithm under different scenarios and different missing degrees, the performance of the WLRTC-P model is evaluated by comparing it with the baseline model.
For the baseline models, to represent the completion performance of the matrix decomposition model (TRMF), and tensor decomposition models (CP_ALS, BGCP, and HaLRTC) on the spatio-temporal dataset, the same settings as in the previous work are followed, and all matrix and tensor decomposition models are configured with the same number of rank functions and iterations, setting the rank of CP_ALS, BGCP, and TRMF to 50 in the random missing case and 10 in the nonrandom missing case. The comparison results are shown in Table1 and Table 2.Table 1
Comparison of accuracy (RMSE/MAPE%) onG and S data sets in the case of random missing.
CP_ALSHaLRTCBGCPTRMFWLRTC-PG20%3.59/8.333.33/8.133.57/8.283.14/7.472.61/5.9840%3.61/8.373.61/8.863.59/8.293.25/7.762.82/6.4550%3.66/8.493.77/9.303.77/9.313.34/80.022.88/6.5760%3.72/8.613.96/9.823.96/9.833.47/8.373.05/6.9670%3.82/8.844.18/10.454.18/10.453.70/89.743.18/7.2480%4.13/9.384.47/11.324.27/10.164.03/9.893.41/7.77S20%4.49/7.423.47/5.934.50/7.453.71/5.963.18/5.0640%4.55/7.583.83/6.764.54/7.583.79/6.163.39/5.4550%4.61/7.714.07/7.304.08/7.313.84/6.293.51/5.6760%4.63/7.774.34/7.904.34/7.913.91/6.463.63/5.9070%4.65/7.784.76/8.894.76/8.894.07/6.843.83/6.2580%4.81/8.125.31/10.254.68/7.844.38/5.164.10/6.78Table 2
Comparison of accuracy (RMSE/MAPE%) onG and S data sets in the case of non-random missing.
CP_ALSHaLRTCBGCPTRMFWLRTC-PG20%4.29/10.274.21/10.454.27/10.204.27/10.244.14/10.0140%4.32/10.284.38/10.884.32/10.254.37/10.374.23/10.1750%4.45/10.484.52/11.304.52/11.314.46/10.574.33/10.4460%4.58/10.734.69/11.804.69/11.814.53/10.804.39/10.5470%4.90/11.094.97/12.664.96/12.654.68/11.314.51/10.9180%5.24/11.645.77/14.615.32/11.884.87/11.854.62/11.08S20%5.63/9.954.69/8.795.65/9.935.26/9.124.50/7.9540%5.70/10.045.27/10.195.68/9.945.30/9.194.85/8.8450%5.78/10.345.64/11.195.64/11.205.38/9.445.01/9.1160%5.89/10.486.11/12.376.11/12.385.53/9.695.26/9.5970%9.31/20.826.94/14.346.91/14.305.72/10.055.69/10.5680%9.96/21.159.77/19.158.57/13.357.10/11.396.23/11.69


#### 5.4.1. Experimental Analysis in the RM Scenario

In the case of random missing, inG and S data sets, from the value of RMSE/MAPE (%), the completion performance of WLRTC-P is significantly better than other classical baseline models at different missing rates (20%–80%), and with the increase of the missing degree, for HaLRTC and BGCP, the accuracy of WLRTC-P model is more obvious, while CP_ALS and TRMF are relatively stable, and the decreasing trend of accuracy is similar to WLRTC-P model. Overall, the performance of the WLRTC-P model is optimal in the case of RM, and in extreme cases, its completion effect is more stable than the baseline models, with strong generalization and robustness.


#### 5.4.2. Experimental Analysis in the NM Scenario

In the case of nonrandom missing, inG and S data sets, the completion accuracy of WLRTC-P is also consistently better than other baseline models. When the missing rate is 20%–60%, the completion accuracy of WLRTC-P is not obvious at this stage, but as the degree of missing rate increases, when the missing is 70%–80%, the completion accuracy advantage of WLRTC-P began to show, and most of the baseline models (CP_ALS, HaLRTC, and BGCP) showed distortion in the case of extreme absence. The completion accuracy of TRMF is also high in the case of extreme missing in G dataset, but the generalization is relatively lower than the WLRTC-P model, and the effect is not good in the extreme case of S dataset. WLRTC-P uses the p-shrinkage norm to eliminate the singular values with less information, retain large data features, and compress enough correlation information. It can still ensure high accuracy in the extreme case of nonrandom shrinkage.Overall, under RM and NM scenarios in both datasets, the WLRTC-P model works optimally for the entire spatio-temporal traffic data recovery problem compared with the baseline model, has the highest data recovery accuracy and model generalizability capability, and performs more significantly in the case of extreme data deficiencies. The three-way weighting approach is able to extract correlation information for each slice direction of the tensor data and use the p-shrinkage norm in each direction to replace the nuclear norm, which is a tighter convex envelope for the rank-minimization problem and is able to compress to get more internal correlation information, and the two optimization methods make the completion accuracy of the WLRTC-P model higher than most baseline models.Here, in two missing scenarios, the completion renderings of the two data sets with the missing rate of 60% are shown in Figure5.Figure 5
Completion effect of G dataset and S dataset, red represents the true value, and blue represents the completion value (missing rate set to 60%). (a) G data set RM, (b) G data set NM, (c) S data set RM, (d) S data set NM.
(a)(b)(c)(d)


## 5.1. Data Set Settings

Experiments are conducted using two publicly available datasets, both collected from real spatio-temporal traffic systems, and experiments on real datasets can verify the applicability of the algorithm in real scenarios. Both datasets have the same tensor pattern “location/sensor × day × time” except for the difference in tensor size.(1)
An urban traffic dataset for Guangzhou, China, for 61 days from August 1, 2016, to September 30, 2016, at 10-minute intervals, consists of 214 anonymous road segments (mainly composed of urban highways and arterial roads). According to the spatio-temporal attributes, establish a third-order tensor with the size of 214 × 61 × 144, where the dimensions of the tensor are road sections, dates, and time windows, respectively. For simplicity, this dataset is referred to as “G” for short.(2)
Seattle freeway traffic speed data collected by Cui et al. [32]. This dataset contains the freeway traffic speed in Seattle, USA, for the whole year of 2015 from 323 loop detectors with a resolution of 5 minutes. Select a subset of January (4 weeks from January 1st to January 28th) as experimental data, and build tensor data with a size of 323 × 28 × 288 according to the space-time attributes, where each dimension of the tensor is road sections, dates, and time windows, respectively. For simplicity, this dataset is referred to as “G” for short.Constructed the traffic data as a tensor pattern of links/sensors, the dates, and time window ensures that the tensor data are sufficiently low rank on each slice, to maximize the use of internal correlation information and thus improve the data restoration accuracy where the link dimension represents the spatial-level features of the traffic data, and the date and time windows represent the time-level features.


## 5.2. Baseline Model

Compared the WLRTC-P model with the following baseline model:HaLRTC [14]high-precision low-rank tensor completion model, the most classical low-rank tensor completion model (LRTC), is based on the minimization of the nuclear norm of the tensor and is solved using the alternating multiplier method (ADMM) framework.BGCP [16]Bayesian Gaussian CP decomposition, a full Bayesian tensor decomposition model, uses Markov chain, Monte Carlo, to learn the potential factor matrix (low-rank structure).CP_ALS [33]CP decomposition model, the classical tensor completion model, is solved by least squares and gradient descent methods, which is more accurate than the original solution method.TRMF [34]temporal regularized matrix decomposition is a temporal matrix decomposition model applied to a multiple autoregressive (AR) process to model the underlying temporal factors.For the selection of baseline model, considered two data construction methods. For matrix-based completion, the TRMF model was selected for comparison, and the data were processed as “location/sensor × time”; for tensor-based completion, the tensor data model was constructed as “location/sensor × day × time,” and the classical CP_ALS, HaLRTC, and BGCP models are selected for comparison. In addition, we have other considerations, the WLRTC-P model does not consider the influence of the directionality of the input tensor slices, and the original model is compared with the proposed model to show the importance of spatio-temporal information.


## 5.3. Experimental Setup

The datasets given in the experiments all tend to be complete, and for the model’s completion accuracy comparison, the true and completion values of the datasets need to be used to calculate the metric, so when conducting the experiments, the data need to be lost according to the set scenario with the missing rate, and then, the missing values are filled by the completion algorithm. For accuracy comparison, the accuracy of completion is judged by comparing the value of MAPE (%) and RMSE. The definitions of MAPE (%) and RMSE are given as (29)(31)MAPE=1n∑i=1nyi−y^iyi×100,RMSE=1n∑i=1nyi−y^i2.In real life, traffic data are missing in two waysrandom missing (RM) and nonrandom missing (NM). The RM case is aimless, mostly caused by unstable signals during communication and transmission, while the NM case data are destroyed in a related way, mostly caused by sensor damage. The two missing scenarios set allowed for a better evaluation of the different model’s performance and effectiveness.Taking two publicly available missing datasets as an example, randomly select a link, to show the data sparsity under 30% and 60% missing spatio-temporal traffic data in a day, as shown in Figure3.Figure 3
Data missing rate is 30% and 60% (Green means missing, and white means no missing issues). (a) G data set RM scenario, (b) S data set RM scenario, (c) G data set NM scenario, (d) S data set NM scenario.
(a)(b)(c)(d)


### 5.3.1. Parameter Setting

Refer to the way of original LRTC framework was set up. For the iteration stepρ, set the initial value to ρ  = 1e8/1e9, set the maximum value ρmax of ρ as 1e5, and update ρ by ρ=min1.05×ρ,ρmax in each iteration update; for the weighting parameter αi, set the sum of parameters to 1, where α1=α2=α3=1/3; μ values to αi/ρk; for p values, set to 0.2; for the convergence condition, use Ml+1−MlF2/PΩMF2<ε to determine whether the algorithm converges; for the convergence accuracy limit E, set the value of E to 1e-4; for the selection of the number of iterations k, set k to 200 can guarantee the convergence of the model.


### 5.3.2. Experimental Scheme

According to the actual collection of traffic data, we set two groups of experiments in the RM scenarios and the NM scenarios. In two missing scenarios, the missing rates ofG and S are set to 20%–80%, respectively. The first experiment is to judge the influence problem of data input direction (slices in different directions) on the completion accuracy by comparing three data sets with different input directions in the same scenario to verify our conjecture. The second is to conduct accuracy comparison experiments, to compare with the baseline model by experiments in the same scenario, and to analyze the performance of the model. In the comparison of the accuracy values in both sets of experiments, the smaller the values of RMSE and MAPE (%), the higher the completion accuracy of the algorithm, and the higher the accuracy advantage and model generalization ability.


## 5.3.1. Parameter Setting

Refer to the way of original LRTC framework was set up. For the iteration stepρ, set the initial value to ρ  = 1e8/1e9, set the maximum value ρmax of ρ as 1e5, and update ρ by ρ=min1.05×ρ,ρmax in each iteration update; for the weighting parameter αi, set the sum of parameters to 1, where α1=α2=α3=1/3; μ values to αi/ρk; for p values, set to 0.2; for the convergence condition, use Ml+1−MlF2/PΩMF2<ε to determine whether the algorithm converges; for the convergence accuracy limit E, set the value of E to 1e-4; for the selection of the number of iterations k, set k to 200 can guarantee the convergence of the model.


## 5.3.2. Experimental Scheme

According to the actual collection of traffic data, we set two groups of experiments in the RM scenarios and the NM scenarios. In two missing scenarios, the missing rates ofG and S are set to 20%–80%, respectively. The first experiment is to judge the influence problem of data input direction (slices in different directions) on the completion accuracy by comparing three data sets with different input directions in the same scenario to verify our conjecture. The second is to conduct accuracy comparison experiments, to compare with the baseline model by experiments in the same scenario, and to analyze the performance of the model. In the comparison of the accuracy values in both sets of experiments, the smaller the values of RMSE and MAPE (%), the higher the completion accuracy of the algorithm, and the higher the accuracy advantage and model generalization ability.


## 5.4. Experiment

Experiment 1.
Influence of data input direction (slicing) on completion accuracy:
In this experiment, according to the above parameter settings, we use the one way of the original model (14) when the data input directions are X1∈ℝn1×n2×n3, X2∈ℝn2×n3×n1, and X3∈ℝn3×n1×n2, and the completion accuracy of traffic data is related to frontal slicing, lateral slicing, and horizontal slicing, and in the experiment, use direction 1, direction 2, and direction 3 to replace them, respectively. The following Figure 4 shows the comparison of slices in three different directions under different missing rates and missing scenarios.
As we can see, in most cases, the accuracy is highest when the data input direction is direction 1, and the completion accuracy is close to direction 1 when the data input direction is direction 2. In the NM scenario of the S data set, the accuracy of direction 2 is better than that of direction 1 in some cases of missing rate, and the accuracy of direction 3 is the worst in comparison. For this reason, we applied the correlation information of other directions in the subsequent WLRTC-P model by way of direction weighting to improve the accuracy. The above experiments verified our conjecture that the different input directions of the traffic data will lead to a significant influence on the accuracy of the subsequent completions.Figure 4
Effect of data input direction (slicing) on completion accuracy (The missing rates of G and S data sets were 20%–80%, respectively). (a) G data set RM scenario, (b) G data set NM scenario, (c) S data set RM scenario, (d) S data set NM scenario.
(a)(b)(c)(d)Experiment 2.
Comparison of completion accuracy under different missing scenarios and missing degrees:
The missing rates ofG and S are set to 20%–80%, respectively. Through experiments, comparing the completion accuracy of the algorithm under different scenarios and different missing degrees, the performance of the WLRTC-P model is evaluated by comparing it with the baseline model.
For the baseline models, to represent the completion performance of the matrix decomposition model (TRMF), and tensor decomposition models (CP_ALS, BGCP, and HaLRTC) on the spatio-temporal dataset, the same settings as in the previous work are followed, and all matrix and tensor decomposition models are configured with the same number of rank functions and iterations, setting the rank of CP_ALS, BGCP, and TRMF to 50 in the random missing case and 10 in the nonrandom missing case. The comparison results are shown in Table1 and Table 2.Table 1
Comparison of accuracy (RMSE/MAPE%) onG and S data sets in the case of random missing.
CP_ALSHaLRTCBGCPTRMFWLRTC-PG20%3.59/8.333.33/8.133.57/8.283.14/7.472.61/5.9840%3.61/8.373.61/8.863.59/8.293.25/7.762.82/6.4550%3.66/8.493.77/9.303.77/9.313.34/80.022.88/6.5760%3.72/8.613.96/9.823.96/9.833.47/8.373.05/6.9670%3.82/8.844.18/10.454.18/10.453.70/89.743.18/7.2480%4.13/9.384.47/11.324.27/10.164.03/9.893.41/7.77S20%4.49/7.423.47/5.934.50/7.453.71/5.963.18/5.0640%4.55/7.583.83/6.764.54/7.583.79/6.163.39/5.4550%4.61/7.714.07/7.304.08/7.313.84/6.293.51/5.6760%4.63/7.774.34/7.904.34/7.913.91/6.463.63/5.9070%4.65/7.784.76/8.894.76/8.894.07/6.843.83/6.2580%4.81/8.125.31/10.254.68/7.844.38/5.164.10/6.78Table 2
Comparison of accuracy (RMSE/MAPE%) onG and S data sets in the case of non-random missing.
CP_ALSHaLRTCBGCPTRMFWLRTC-PG20%4.29/10.274.21/10.454.27/10.204.27/10.244.14/10.0140%4.32/10.284.38/10.884.32/10.254.37/10.374.23/10.1750%4.45/10.484.52/11.304.52/11.314.46/10.574.33/10.4460%4.58/10.734.69/11.804.69/11.814.53/10.804.39/10.5470%4.90/11.094.97/12.664.96/12.654.68/11.314.51/10.9180%5.24/11.645.77/14.615.32/11.884.87/11.854.62/11.08S20%5.63/9.954.69/8.795.65/9.935.26/9.124.50/7.9540%5.70/10.045.27/10.195.68/9.945.30/9.194.85/8.8450%5.78/10.345.64/11.195.64/11.205.38/9.445.01/9.1160%5.89/10.486.11/12.376.11/12.385.53/9.695.26/9.5970%9.31/20.826.94/14.346.91/14.305.72/10.055.69/10.5680%9.96/21.159.77/19.158.57/13.357.10/11.396.23/11.69


### 5.4.1. Experimental Analysis in the RM Scenario

In the case of random missing, inG and S data sets, from the value of RMSE/MAPE (%), the completion performance of WLRTC-P is significantly better than other classical baseline models at different missing rates (20%–80%), and with the increase of the missing degree, for HaLRTC and BGCP, the accuracy of WLRTC-P model is more obvious, while CP_ALS and TRMF are relatively stable, and the decreasing trend of accuracy is similar to WLRTC-P model. Overall, the performance of the WLRTC-P model is optimal in the case of RM, and in extreme cases, its completion effect is more stable than the baseline models, with strong generalization and robustness.


### 5.4.2. Experimental Analysis in the NM Scenario

In the case of nonrandom missing, inG and S data sets, the completion accuracy of WLRTC-P is also consistently better than other baseline models. When the missing rate is 20%–60%, the completion accuracy of WLRTC-P is not obvious at this stage, but as the degree of missing rate increases, when the missing is 70%–80%, the completion accuracy advantage of WLRTC-P began to show, and most of the baseline models (CP_ALS, HaLRTC, and BGCP) showed distortion in the case of extreme absence. The completion accuracy of TRMF is also high in the case of extreme missing in G dataset, but the generalization is relatively lower than the WLRTC-P model, and the effect is not good in the extreme case of S dataset. WLRTC-P uses the p-shrinkage norm to eliminate the singular values with less information, retain large data features, and compress enough correlation information. It can still ensure high accuracy in the extreme case of nonrandom shrinkage.Overall, under RM and NM scenarios in both datasets, the WLRTC-P model works optimally for the entire spatio-temporal traffic data recovery problem compared with the baseline model, has the highest data recovery accuracy and model generalizability capability, and performs more significantly in the case of extreme data deficiencies. The three-way weighting approach is able to extract correlation information for each slice direction of the tensor data and use the p-shrinkage norm in each direction to replace the nuclear norm, which is a tighter convex envelope for the rank-minimization problem and is able to compress to get more internal correlation information, and the two optimization methods make the completion accuracy of the WLRTC-P model higher than most baseline models.Here, in two missing scenarios, the completion renderings of the two data sets with the missing rate of 60% are shown in Figure5.Figure 5
Completion effect of G dataset and S dataset, red represents the true value, and blue represents the completion value (missing rate set to 60%). (a) G data set RM, (b) G data set NM, (c) S data set RM, (d) S data set NM.
(a)(b)(c)(d)


## 5.4.1. Experimental Analysis in the RM Scenario

In the case of random missing, inG and S data sets, from the value of RMSE/MAPE (%), the completion performance of WLRTC-P is significantly better than other classical baseline models at different missing rates (20%–80%), and with the increase of the missing degree, for HaLRTC and BGCP, the accuracy of WLRTC-P model is more obvious, while CP_ALS and TRMF are relatively stable, and the decreasing trend of accuracy is similar to WLRTC-P model. Overall, the performance of the WLRTC-P model is optimal in the case of RM, and in extreme cases, its completion effect is more stable than the baseline models, with strong generalization and robustness.


## 5.4.2. Experimental Analysis in the NM Scenario

In the case of nonrandom missing, inG and S data sets, the completion accuracy of WLRTC-P is also consistently better than other baseline models. When the missing rate is 20%–60%, the completion accuracy of WLRTC-P is not obvious at this stage, but as the degree of missing rate increases, when the missing is 70%–80%, the completion accuracy advantage of WLRTC-P began to show, and most of the baseline models (CP_ALS, HaLRTC, and BGCP) showed distortion in the case of extreme absence. The completion accuracy of TRMF is also high in the case of extreme missing in G dataset, but the generalization is relatively lower than the WLRTC-P model, and the effect is not good in the extreme case of S dataset. WLRTC-P uses the p-shrinkage norm to eliminate the singular values with less information, retain large data features, and compress enough correlation information. It can still ensure high accuracy in the extreme case of nonrandom shrinkage.Overall, under RM and NM scenarios in both datasets, the WLRTC-P model works optimally for the entire spatio-temporal traffic data recovery problem compared with the baseline model, has the highest data recovery accuracy and model generalizability capability, and performs more significantly in the case of extreme data deficiencies. The three-way weighting approach is able to extract correlation information for each slice direction of the tensor data and use the p-shrinkage norm in each direction to replace the nuclear norm, which is a tighter convex envelope for the rank-minimization problem and is able to compress to get more internal correlation information, and the two optimization methods make the completion accuracy of the WLRTC-P model higher than most baseline models.Here, in two missing scenarios, the completion renderings of the two data sets with the missing rate of 60% are shown in Figure5.Figure 5
Completion effect of G dataset and S dataset, red represents the true value, and blue represents the completion value (missing rate set to 60%). (a) G data set RM, (b) G data set NM, (c) S data set RM, (d) S data set NM.
(a)(b)(c)(d)


## 6. Conclusions

In this paper, we proposed the WLRTC-P model, which reduces the dependence of the completion accuracy on the data input direction by direction weighting, extracts the data correlation information, improves the completion accuracy by using the p-shrinkage norm as a convex proxy of the tensor average rank, and solves it by the alternating direction method of multipliers. Experiments are conducted in two publicly available datasets with random missing scenarios and nonrandom missing scenarios. The data input directionality experiment shows the necessity of weighting. Subsequent comparison of completion accuracy with the current classical model shows that our proposed weighted model achieved good results in most scenarios.



---

*Source: 3304677-2022-10-13.xml*