1. Introduction

In the past few decades, due to the importance of position based services in various practical applications such as smart home, robot navigation, library route guidance [1], large-scale supermarket shopping [2] and so on, the demand for indoor position based services has been increasing. Many positioning systems have been studied to estimate the position of objects in the environment. For outdoor environments, global positioning system(GPS) can achieve high accuracy positioning. However, GPS cannot work reliably indoors due to the blockage of buildings. Accordingly, some indoor positioning systems have been developed, including technologies such as Wi-Fi [3], radio-frequency identification [4], ultra-wideband [5], and ZigBee [6]. Unfortunately, these technologies have limitations in terms of positioning accuracy and stability due to multipath effect.

In recent years, visible light communication(VLC) technology has shown great potential for high-precision indoor positioning [7]. VLC uses light in the visible band as the information carrier for communication, which has the advantages of low energy consumption and high security [8]. For visible light positioning(VLP), the most commonly used signal characteristic is the received signal strength (RSS), correspondingly, photodiode (PD) is usually used as the optical signal detector. Visible light-based localization algorithms mainly include proximity-, fingerprinting-, and geometry-based methods [9]. Fingerprint-based positioning systems have been widely studied for the feasibility and low complexity [10], which can obtain accurate positions by simple calibration.

Owing to factors such as obstacle occlusions and sunlight influences in the environment, it is actually difficult to obtain an accurate channel model. Therefore, it is hard to effectively construct the corresponding fingerprint database. An alternative way is to arrange dense PDs to make detections without the knowledge of the channel model. However, due to the space-limited scenario and high cost, sparse PD arrangement is indispensable. Hence, how to recover the fingerprint database from the sparse PD measurements is a challenge.

Recovering the fingerprint database through sparse PD measurements can be regarded as a sparse data filling problem. At present, some data recovery algorithms have been widely studied. For each sample to be determined, K-nearest neighbors (KNN) algorithm [11,12] achieves classification or numerical prediction by searching for the $K$ nearest samples in the feature space. Accordingly, this method requires extensive computation. Singular value threshold(SVT) algorithms [13–15] perform singular value decomposition on the matrix to restore the sparse matrix iteratively. It can not fully utilize the multidimensional characteristics of data, which results in unsatisfactory accuracy. Compressed sensing(CS) algorithms [16–18] creatively combine the norm minimization with undersampling to reconstruct sparse signals. But its performance is limited to the used random sampling matrix. Matrix filling algorithms [19] and tensor completion algorithms [20,21] recover data based on the low-rank feature of matrices or tensors. However, the accuracy of these algorithms could be further improved, as they fail to fully exploit the potential spatio-temporal correlations of the data. More importantly, the PDs in the environment are usually sparsely arranged beforehand, which results in missing data in the locations without PDs. Unfortunately, the data recovery algorithms in the recent literature focus more on random missing, rather than the missing form caused by sparsely arranged PDs. To solve this problem, a new tensor completion algorithm based on CANDECOMP/PARAFAC(CP) decomposition is proposed in this paper. It is worth noting that Abou-Shehada et al. [22] proposed a weighted KNN approach for fingerprint recognition in VLP. They used the sparse RSS dataset to obtained the fingerprint database by pathloss exponent calibration and finally achieved accurate positioning. However, the pathloss calibration is still based on the availability of the channel model. The advantage of our paper is that no calibration is required and the fingerprint database can be directly and accurately restored by the proposed algorithm, thus allowing it to be applied in a wide range of scenarios. The main contributions in this paper are summarized as follows:

The remainder of the paper is organized as follows. The preliminaries are presented in Section 2. Section 3 introduces the system model and the low-rank feature of the data tensor. Section 4 details the proposed algorithm and the process of solving optimization problems. In section 5, we evaluate the performance of the proposed algorithm via experimental data. Finally, we summarized the work in section 6.

2. Preliminaries

In this section, some basic notations and necessary definitions are presented. The following notations are used throughout this paper. Tensors are denoted by uppercase calligraphic letters (${\cal X,Y,Z\ldots }$), matrices are written as uppercase boldfaced letters(X,Y,Z…), vectors are represented by lowercase boldfaced letters (x,y,z…). Symbol $\lVert. \rVert _0$, $\lVert. \rVert _1$ and $\lVert. \rVert _2$ respectively stands for the L$_0$, L$_1$ and L$_2$ norm. Subsequently, symbol $\lVert. \rVert _F$ and $\lVert. \rVert _*$ denotes the Frobenius norm and nuclear norm respectively. The relevant variables used in this paper are summarized in Table 1.

Table 1. Key variables in this paper

View Table | View all tables in this article

Definition 1 (Kronecker product of matrix [23]): The Kronecker product of matrix A$\in \mathbb {R}^{m \times n}$ and B$\in \mathbb {R}^{M \times N}$ is defined as A$\otimes$B=[$a_{ij}\cdot$B]$\in \mathbb {R}^{mM \times nN}$, in which $a_{ij}\in$A.

Definition 2 (Khatri-Rao product of matrix): The Khatri-Rao product of matrix A$\in \mathbb {R}^{I \times K}$ and B$\in \mathbb {R}^{J \times K}$ is defined as A$\odot$B=[a$_{1}\otimes$b$_{1}$ a$_{2}\otimes$b$_{2}$ …a$_{K}\otimes$b$_{K}$]$\in \mathbb {R}^{IJ \times K}$, in which ${\textbf {a}}_{k}$ and ${\textbf {b}}_{k}$ respectively denotes a column of matrix A and B, ${k}\;\in \;\left \{ {1,2 \ldots,K} \right \}$.

Definition 3 (Slice of tensor [24]): A $3$-order tensor can be represented by a set of two-dimensional slices. For tensor ${\cal X}\;\in \;\mathbb {R}^{L \times W \times T}$, slices in various directions can be obtained by observing from different dimensions, denoted by $\textbf{X}_{l,:,:}$, $\textbf{X}_{:,w,:}$, and $\textbf{X}_{:,:,t}$ respectively, where $l\in \left \{ {1,2,\ldots,L} \right \}$, $w\in \left \{ {1,2,\ldots,W} \right \}$ and $t\in \left \{ {1,2,\ldots,T} \right \}$.

Definition 4 (Tensor expansion by mode-$n$ [25]): Mode-$n$ expansion is to unfold tensor into matrix format according to a predefined rule. Generally, mode-$n$ expansion of tensor ${\cal X}_N\;\in \;\mathbb {R}^{I_1 \times I_2 \times \ldots \times I_N}$ along its $n$-th mode into a matrix is defined as X$_{N(n)}\;\in \;\mathbb {R}^{I_n \times \prod _{j\not =n} I_j}$, $n\in \left \{ {1,2,\ldots,N} \right \}$. As shown in Fig. 1, we illustrate the mode-$1$, mode-$2$ and mode-$3$ expansion of the $3$-order tensor ${\cal X}\;\in \;\mathbb {R}^{L \times W \times T}$.

 

Fig. 1. Schematic diagram of mode-$1$, $2$, $3$ expansion.

Download Full Size | PDF

Definition 5 (Rank-1 tensor): If tensor ${\cal X}_N\;\in \;\mathbb {R}^{I_1 \times I_2 \times \ldots \times I_N}$ can be written as the outer product of $N$ vectors, it is called rank-1 tensor. Accordingly, rank-1 tensor ${\cal X}_N$ can be denoted as ${\cal X}_N$=${\textbf{v}^{(1)} \circ \textbf{v}^{(2)} \circ \ldots \circ \textbf{v}^{(N)}}$, in which ${\textbf{v}^{(n)}}\in \mathbb {R}^{I_n}$ and $n\in \left \{ {1,2,\ldots,N} \right \}$.

Definition 6 (Tensor rank): The rank of tensor is defined as the minimum number of rank-1 tensor required to restore the current tensor. A tensor ${\cal X}_N$ with rank $R$ can be expressed as ${\cal X}_N$=$\sum\limits _{r = 1}^{R}{\textbf{v}_r^{(1)} \circ \textbf{v}_r^{(2)} \circ \ldots \circ \textbf{v}_r^{(N)}}$, in which ${\textbf{v}_r^{(n)}}\in \mathbb {R}^{I_n}$ and $n\in \left \{ {1,2,\ldots,N} \right \}$.

Definition 7 (CP decomposition [26]): CP decomposition is to decompose a tensor into the sum of a series of rank-1 tensors, by this way, the $N$-order tensor ${\cal X}_N\;\in \;\mathbb {R}^{I_1 \times I_2 \times \ldots \times I_N}$ with rank $R$ can be decomposed to ${\cal X}_N$=$\sum\limits _{r = 1}^{R}{\textbf{v}_r^{(1)} \circ \textbf{v}_r^{(2)} \circ \ldots \circ \textbf{v}_r^{(N)}}$. By introducing the factor matrices V$_{(n)}$=$[\textbf{v}_1^{(n)},\textbf{v}_2^{(n)},\ldots,\textbf{v}_R^{(n)}]\in \mathbb {R}^{I_n\times R}$, tensor ${\cal X}$ can be further denoted as ${\cal X}_N$=$\left [\kern -0.15em\left [ {\textbf {V}^{(1)},\textbf {V}^{(2)},\ldots,\textbf {V}^{(N)}} \right ]\kern -0.15em\right ]$.

3. System model and problem formulation

In this section, we illustrate the system configuration diagram, including the layout of LED luminaires and PDs. The illuminance received by PDs is modeled as a tensor. Furthermore, we analyze the low-rank feature of the tensor. Based on that, the corresponding optimization problem is formulated.

3.1 System configuration

Considering an indoor scenario, the LED is mounted in the geometric centre of the ceiling. The PDs are placed on the floor to measure illuminance. The scenario diagram of the indoor VLC system is shown in Fig. 2.

 

Fig. 2. The VLC system model.

Download Full Size | PDF

In the VLC network, it is assumed that the receiver plane of each PD is paralleled to the ceilling. The vertical distance $h$ from the LED to the ground is constant. Since the light signal intensity from reflection and diffusion is much less than that from the line of sight(LOS) transmission in the VLC environment [27], the direct current (DC) channel gain of the LOS link between the LED and PD$_i$ is given as [28],

The gain of optical concentrator ${g({\psi _{i}})}$ is given as [29]:

The affecting factors of the missing data in the fingerprint database mainly include: (1) the size of the grids division, (2) the sparseness of the sensors arrangement, (3) the sampling time, etc. The illuminance in the environment varies continuously with space and time. However, we only obtain the illuminance by arranging sparse PDs on the discrete grids, which results in missing data for grids without PD. Consequently, the incomplete data are definitely introduced. The location set of all grids is defined as:

To distinguish the different scenarios that our algorithm and most other algorithms focus on, here we define two forms of data missing, including random missing (RM) and sensor time series missing (STSM). Arranging sensors at all grids to sample randomly at certain slots results in RM. However, it is difficult to achieve in practice. If sparse sensors are placed in the environment to sample data for a period of time simultaneously, it results in STSM. We construct the measured data as a third-order tensor $\cal X\;\in \mathbb {R} ^{L,W,T}$, where $T$ denotes the number of slots.

Schematic diagrams of RM and STSM are shown in Fig. 3(a) and Fig. 3(b) respectively. The orange cubes represent the measured data and the white cubes denote the missing data. In the practical scenario, only sparse PDs can be arranged. Hence, we focus on the scenario of STSM.

 

Fig. 3. Two missing forms of the data tensor: RM and STSM.

Download Full Size | PDF

3.2 Low rank analysis

The low-rank feature of the tensor provides the theoretical basis for data recovery. Tensor singular value decomposition(t-SVD) is a new decomposition framework based on the tensor product operation [30], mainly for the decomposition of 3-D tensor. Unlike CP decomposition or Tucker decomposition, the advantage of t-SVD is that the singular values of the tensor can be obtained intuitively. The t-SVD decomposition model of a third-order tensor is shown in Fig. 4.

 

Fig. 4. Schematic diagram of t-SVD.

Download Full Size | PDF

Accordingly, the t-SVD decomposition of $\cal X$ can be expressed as follows:

The singular values have been sorted from the largest to the smallest in the process of t-SVD. The sum of the first $n$ singular values is $E_n$, the sum of all singular values is $E_{total}$, and we define the ratio $E_n/E_{total}$ as the energy accumulation ratio. Accordingly, to illustrate the effect of missing data on the low-rank feature of the tensor, we further analyze the variation of the energy accumulation ratio with the number of singular values in both RM and STSM scenarios.

As shown in Fig. 5, only a few singular values are needed to take the energy accumulation rate close to 1, which reflects the low-rank feature of the tensor. Moreover, comparing Fig. 5(a) and (b), the energy accumulation rate in the STSM scenario rises more slowly than in the RM scenario. It can be concluded that the low-rank feature of the tensor is certainly destroyed in the STSM scenario.

 

Fig. 5. The change of energy accumulation rate under different data missing ratios.

Download Full Size | PDF

3.3 Problem formulation

The fingerprint database is recovered by measured data of sparse PDs, which can actually be transformed into a tensor completion problem. Based on the low-rank feature of the data tensor, we propose the corresponding optimization problem for recovering data tensor.

Before introducing the low rank tensor completion problem, let us start with the well-known low rank matrix completion [31]:

Inspired by the low rank matrix completion problem, the low rank tensor completion problem can be expressed as:

4. Proposed SCTC algorithm

In this section, we propose the SCTC algorithm based on CP decomposition to recover sparse PDs measured data. Then the procedure of constructing the spatial constraint matrix based on graph theory is given. Subsequently, the Toeplitz matrix is used as the temporal constraint matrix.

4.1 Spatio-temporal constraint tensor comlpetion based on CP decomposition

Note that although the rank of matrix in optimization problem (7) is NP-hard [32] and difficult to solve, since the nuclear norm is the tightest convex envelope of the matrix rank, some literature proves that currently using nuclear norm as a substitution is still the optimal solution [33,34]. Accordingly, applying this method to tensors, we deduce the optimization problem (7) as follows:

Since the nuclear norm is proposed for matrices, here we give the definition of the nuclear norm of the $n$-th order tensor ${\cal X}_n$:

In essence, we can understand the nuclear norm of the tensor as some linear convex combination of the nuclear norm of all matrices expanded along each mode. It should be noted that, when the mode number $n$=2, the nuclear norm of tensor is completely consistent with the definition of the matrix. After introducing the definition of nuclear norm, we can further rewrite the optimization problem as:

The tensor obtained from the total $T$ slots measured data is ${\cal X}_T$. Accordingly, we introduce a sampling operator ${\cal P}_{\Omega }$ to limit the consistency of elements at the same positions in different tensors, which can be denoted as:

By introducing the constraint term as a regularization into the optimization problem, we obtain the new optimization problem as:

The rank minimization constraint only ensures that the tensor is recovered with a minimal number of rank-one tensor to approximate the original tensor. However, it only obtains coarse optimization results. It is indispensable to further impose constraint on the optimization problem. However, imposing constraints on the nuclear norm is not easy. Therefore, the CP decomposition of tensor ${\cal X}$ is introduced.

By introducing the rank adjustment factor $\lambda$, the previous regularization factor $\rho$ can be substituted. We separate the three factor matrices and use the Frobenius norm to replace the nuclear norm. The optimization problem can be further transformed into:

Such an optimization problem can recover RM data. However, the recovery accuracy is limited for STSM data. Accordingly, we introduce spatial and temporal constraint matrices for the optimization problem. The spatio-temporal correlation of the tensor and the construction of the spatial and temporal constraint matrices would be described in the next subsection. Thus, we obtain the new optimization problem $f$ as:

The optimization problem in Eq. (17) is nonconvex and difficult to solve for the joint variables $\textbf{I},J$ and $\textbf{K}$. Therefore, we decompose the original optimization problem into three sub-optimization problems $f_1,f_2,f_3$ to solve $\bf {I,J,K}$ separately:

For Eq. (18), we solve it by alternating optimization. Specifically, the $\textbf{J}$ and $\textbf{K}$ matrices are fixed when solving for $\textbf{I}$. Since the optimization problem $f_1$ is convex and derivable with respect to matrix $\textbf{I}$, we directly derive $f_1$ as follows:

Then, we let the derivative of $\frac {\partial f_1}{\partial \textbf{I}} = 0$ to obtain:

Similarly, the optimization problems $f_2$ and $f_3$ are convex and derivable with respect to matrix $\textbf{I}$ and $\textbf{K}$ respectively, we directly derive $f_2$ and $f_3$ as follows:

Let the derivative of $\frac {\partial f_2}{\partial \textbf{J}} = 0$ and $\frac {\partial f_3}{\partial \textbf{K}} = 0$, we obtain:

Algorithm 1. Spatio-temporal Constraint Tensor Completion (SCTC)

View Table | View all tables in this article

Due to the Eq. (20),(23) and (24) are complex and difficult to solve, we use the conjugate gradient method and solve alternately to obtain $\bf {I},\bf {J}$ and $\bf {K}$. Finally, the estimated tensor $\hat {\cal X}$ is recovered through the factor matrices $\bf {I}$, $\bf {J}$ and $\bf {K}$:

The pseudocode of the proposed SCTC algorithm is illustrated in Algorithm  1.

4.2 Construction of constraint matrix

In this subsection, we firstly analyze the spatio-temporal correlation of data tensor. Subsequently, we illustrate the construction of the spatial and temporal constraint matrices.

4.2.1 Temporal constraint matrix

We calculate the spatial and temporal correlation for the data tensor. Specifically, we divide the data tensor into slice matrices and calculate the correlation from three dimensions of the tensor separately. The Pearson correlation coefficient is used to measure the correlation between the two slice matrices. For matrices $\textbf{{P,Q}} \in \mathbb {R}^{M \times N}$, the Pearson correlation coefficient is defined as follows:

For example, we divide the tensor into $T$ time slice matrices on the time dimension. We calculate the correlation coefficients for any two time slice matrices to obtain the correlation schematic in Fig. 6(a). It can be seen that there is a high correlation between the two adjacent time slice matrices. This is because the illuminance is relatively stable over time in the absence of strong fluctuations from the outside.

 

Fig. 6. Temporal and spatial correlations of the considered tensor.

Download Full Size | PDF

Similarly, Fig. 6(b) and (c) show the high correlation among the slice matrices on spatial dimensions including horizontal and vertical directions. Accordingly, due to the symmetry of the data in the vertical and horizontal directions, a strong correlation is shown near the main diagonal in Fig. 6(b) and (c). Based on the high temporal and spatial correlation of the data tensor, we further design the corresponding constraint matrices.

According to the temporal correlation between any two time slice matrices in Fig. 6(a), the Toeplitz matrix is designed according to the temporal correlation of adjacent moments.

For $f_3$ in the optimization problem Eq. (18), the value of the constraint term should be as small as possible to avoid additional errors. The data measured by a PD corresponds to a column in the mode-$3$ expansion of the data tensor. Therefore, the Toeplitz matrix is used as the temporal constraint matrix multiplied with the mode-$3$ expansion, which maintains the temporal correlation among the data while solving the optimization problem.

4.2.2 Spatial constraint matrix

As shown in Fig. 6, there is a high correlation among the data tensor on the spatial dimension. We utilize the idea of undirected graph [35] to describe the relationship of spatial grids. The undirected graph can be expressed as $G = (V, E)$, where $V$ is a set of vertices representing grids and $E$ is a set of edges. And the correlation can be described by Laplace matrix $\textbf{L}$, which is defined as $\textbf{L} = D-\textbf{A}_d$, where $\textbf{D}$ and $\textbf{A}_d$ are the degree matrix and adjacency matrix of the graph respectively. Thus, we construct spatial constraint matrices based on Laplace transform.

For the data tensor, the rows of mode-$1$ expansion matrix correspond to the rows of grids, and the columns of mode-$2$ expansion matrix correspond to the columns of grids. Hence, we utilize $\textbf{X}_{(1)}$ and $\textbf{X}_{(2)}$ to construct the spatial constraint matrices $\textbf{S}_1$ and $\textbf{S}_2$, respectively. The process of constructing $\textbf{S}_1$ is as follows.

Firstly, We calculate the correlation coefficient between any two rows of $\textbf{X}_{(1)}$. The correlation coefficient $\rho _{ij}$ between two vectors ${{\textbf{x}_i}}$ and ${{\textbf{x}_j}}$ is defined as:

Referring to the undirected graph theory, the adjacency matrix $\textbf{A}_d\in \mathbb {R}^{L\times L}$ can be expressed as:

We set the threshold $u$ of the correlation coefficient to 0.9. If the correlation coefficient $\rho _{ij}$ between ${{\textbf{x}_i}}$ and ${{\textbf{x}_j}}$ is greater than $u$, let $a_{ij}=1$, otherwise $a_{ij}=0$, $i,j \in \left \{ {1,2, \ldots,L} \right \}$. It can be expressed as:

In particular, since we do not consider the autocorrelation of each row vector, all elements on the diagonal of the matrix $\textbf{A}_d$ are 0, i.e., $a_{ii}=0$.

Next we calculate the sum of the elements of each column of the adjacency matrix $\textbf{A}_d$. Accordingly, the sum are placed on the diagonal of the corresponding column separately, and the elements on the other positions are zero. The matrix obtained is called the degree matrix $\bf {D}$.

Obviously, $\bf {D}$ is a diagonal matrix. Finally, we obtain the spatial constraint matrix $\textbf{S}_1$ by subtracting the degree matrix $\bf {D}$ from the adjacency matrix $\textbf{A}_d$.

The process of constructing the spatial constraint matrix $\textbf{S}_1$ is illustrated in Algorithm  2.

Algorithm 2. The label and the caption are preserved here, and the body of the algorithm can be found below.

View Table | View all tables in this article

The construction method of the spatial constraint matrix $\textbf{S}_2$ is similar to $\textbf{S}_1$, and it is sufficient to use the mode-$2$ expansion of the tensor $\cal X$.

5. Experiment results

In this section, we validate the performance of the proposed algorithm using the experimental measured data.

5.1 Experimental setup

5.1.1 Dataset

A $1.2\times 0.6\times 1.5 m^3$ VLC system is shown in Fig. 7. The LED luminaire is embedded in the center of the ceiling, and the floor is divided into square grids with $0.1m$ side length. Each PD is placed in the center of the grid. The received signal is the illuminance in lux. The entire experimental environment was covered with a black cloth to reduce the impact of light and pedestrians in the environment.

 

Fig. 7. Experimental setup.

Download Full Size | PDF

Accordingly, the relevant parameters are given in Table 2.

Table 2. Experimental parameters

View Table | View all tables in this article

We model the measured data as tensor ${\cal X}\in \;\mathbb {R}^{L\times W\times T}$. $L$=12 and $W$=6 are the number of horizontal and vertical grids respectively. $T$ = 120 represents the number of time slots. The element ${\cal X}_{l,w,t}$, $l \in \left \{ {1,2, \ldots,L} \right \}$, $w \in \left \{ {1,2, \ldots,W} \right \}$, $t \in \left \{ {1,2, \ldots,T} \right \}$ in tensor represents the illuminance of the $t$-th time slot at the $l$-th and $w$-th grid.

5.1.2 Parameter setting

Considering the optimization problem of Eq. (18), the values of each regulator are empirically set as $\beta _1$=$\beta _2$=$\beta _3=6\times 10^{1}$, $\lambda =0.3$. Meanwhile, we will analyze the errors of different ranks to determine the appropriate rank value.

5.1.3 Error evaluation criteria

In order to quantitatively evaluate the performance of the proposed algorithm, we use relative mean square error (RMSE) to measure the accuracy of data recovery:

5.2 Performance analysis

We present the recovery results of the proposed SCTC algorithm using the measured data. Subsequently, we analyze the performance from various aspects, such as algorithm convergence and error.

5.2.1 Robustness analysis

Figure 8 shows the error analysis of $R$=4,6 and 8 in the scenario of STSM. It can be seen that for a certain rank, the RMSE increases with the increase of the missing proportion. Moreover, when the missing proportion is fixed, the RMSE decreases as the rank increases. This is because when the rank increases, more rank-one tensors can be used to construct the original tensor, more accurate recovery results are obtained.

 

Fig. 8. Variation of error with missing rate at $R$=4, 6 and 8.

Download Full Size | PDF

In general, the RMSE variation among the different ranks is small, and all are within an acceptable range. Since spatio-temporal constraint matrices are introduced in the algorithm, which can effectively constrain the data recovery process, the algorithm is robust for different ranks. Accordingly, the results illustrate that the proposed SCTC algorithm is robust and recovers well for different ranks.

5.2.2 Convergence analysis

The convergence analysis of the proposed algorithm is given in Fig. 9. It can be seen that as the missing rate increases, the RMSE increases and the convergence speed slows down. When the missing rate is fixed, the RMSE decreases as the rank increases. Because when more rank-one tensors are used to restore the original tensor, more accurate results can be obtained with less error. Moreover, the convergence curves show that the algorithm converges basically within 5 iterations. Hence, we can conclude that the proposed SCTC algorithm has good convergence and robustness. Considering the errors and computational complexity under different ranks, we set $R$=6 for subsequent experiments.

 

Fig. 9. The convergence of the algorithm when R=4,6 and 8 using experimental data.

Download Full Size | PDF

5.2.3 Recovery error comparison

In this section, we compare the performance of the proposed SCTC algorithm with several commonly used data recovery algorithms. The baseline data recovery algorithms used here are mainly divided into the following four categories: (1)SVT decomposition algorithm [36]; (2) CS algorithm [37]; (3)High accuracy low rank tensor completion(HaLRTC) [31]. (4)Low-Rank High-Order Tensor Completion algorithm With Applications in Visual Data (HTNN-FFT) [38]. For comparison, we also used the algorithm based on CP decomposition, which does not involve utilizing the spatio-temporal correlation of the data tensor.

Through experiment, we compare the recovery performance of the proposed SCTC algorithm and several other baseline data recovery algorithms in the STSM scenario, where the missing ratio varies from 0.1 to 0.9. Experimental results of data recovery errors are given in Fig. 10. Obviously, the RMSE of various recovery algorithms increases as the missing rate increases. Since the proposed SCTC algorithm introduces spatio-temporal constraint matrices, which can provide effective guidance in the process of data recovery to avoid serious distortion of the recovered data. In general, the proposed SCTC algorithm has better stability and a smaller RMSE, especially when the proportion of missing data increases.

 

Fig. 10. Performance comparison of various recovery algorithms under STSM scenario.

Download Full Size | PDF

5.2.4 Time complexity analysis

The time complexity of various algorithms is shown in Fig. 11. Overall, the time required for various algorithms is not significant due to the small amount of data used. Since the proposed SCTC algorithm converges only within a few iterations, which reduces a large amount of computational and time complexity. Moreover, compared with the recovery algorithm using only CP decomposition, faster convergence can be obtained because the proposed algorithm imposes spatial and temporal constraint matrices. In general, the time complexity of our algorithm is significantly lower than other algorithms.

 

Fig. 11. Comparison of the calculation time of different data recovery algorithms.

Download Full Size | PDF

6. Conclusion

In this paper, we proposed a SCTC algorithm to solve the problem of missing data caused by sparsely arranged PDs. Reasonable constraint for the optimization problem in this algorithm was constructed by introducing spatial and temporal constraint matrices. Then we validated the performance of the proposed algorithm using an experimental measured dataset. The results showed that the proposed algorithm outperformed other algorithms in terms of RMSE and time complexity.