Efficient Data Gathering Methods in Wireless Sensor Networks Using GBTR Matrix Completion

To obtain efficient data gathering methods for wireless sensor networks (WSNs), a novel graph based transform regularized (GBTR) matrix completion algorithm is proposed. The graph based transform sparsity of the sensed data is explored, which is also considered as a penalty term in the matrix completion problem. The proposed GBTR-ADMM algorithm utilizes the alternating direction method of multipliers (ADMM) in an iterative procedure to solve the constrained optimization problem. Since the performance of the ADMM method is sensitive to the number of constraints, the GBTR-A2DM2 algorithm obtained to accelerate the convergence of GBTR-ADMM. GBTR-A2DM2 benefits from merging two constraint conditions into one as well as using a restart rule. The theoretical analysis shows the proposed algorithms obtain satisfactory time complexity. Extensive simulation results verify that our proposed algorithms outperform the state of the art algorithms for data collection problems in WSNs in respect to recovery accuracy, convergence rate, and energy consumption.


Introduction
Wireless sensor networks (WSNs) are composed of large-scale, self-organized sensor nodes, which are capable of sensing, data storage, and communication. WSNs have lots of applications, such as remote environment sensing, industrial automation, smart city, and military monitoring. In practical applications, lots of ordinary nodes are deployed in an unattended mode. These ordinary nodes perform data collection tasks individually and transmit the raw data to the sink node in multi-hop access. Since it is difficult to recharge or replace the limited power supply of ordinary nodes, developing energy efficient data gathering methods is becoming crucial.
A large number of data collection methods have been proposed to reduce the energy consumption with different levels of data reconstruction precision in the literature [1][2][3]. These obtained data in WSNs possess spatial and temporal correlations, which are intrinsic characteristics of a physical environment. A previous article [1] proposed a clustered aggregation (CAG) algorithm for data collection, which utilizes the spatial and temporal correlations of the sensed data. Pham et al. [2] presented a divide and conquer approximating (DCA) algorithm to reduce power consumption. Since the sensed data require to be transmitted to the sink node in multihop communication, Rosana et al. [3] proposed a novel algorithm to construct spanning trees for efficient data gathering in wireless sensor networks. Unfortunately, data gathering methods in traditional mode have limitations. Firstly, the clustering methods (or the spanning tree construction methods) reflect high computational (1) The features of sensor datasets are analyzed in consideration of their topology information, which reveals that the data matrix is sparse under the graph based transform. (2) The graph based transform regularized (GBTR) Matrix Completion problem is formulated.
To reconstruct the missing values efficiently, the GBTR by Alternating Direction Method of Multipliers (GBTR-ADMM) algorithm is proposed. Simulation results reveal that GBTR-ADMM outperforms the state of art algorithms in view of the recovery accuracy and the energy consumption. (3) To accelerate the convergence of GBTR-ADMM, GBTR-A2DM2 algorithm is proposed, which benefits from a restart rule and the fusion of multiple constraints. (4) The time complexity of our proposed algorithms is analyzed, which shows that the complexity is low.
The structure of the paper is concluded as follows: In Section 2, the problem formulation about matrix completion is given. Section 3 presents the features of the real datasets and the synthesized dataset. The proposed GBTR-ADMM and GBTR-A2DM2 algorithms are expatiated individually in Sections 4 and 5. Section 6 shows the time complexity of the proposed algorithms. In Section 7, the performances of the proposed algorithms are studied. The conclusions and the future works are summarized in Section 8.

Problem Formulation
In this section, we introduce the related issues in respect to matrix completion theory. The main notations of the paper are summarized in Table 1. The observed ratio r The matrix rank λ The GBT sparsity regularization parameter β The Lagrange penalty parameter X The original data matrix X The reconstructed data matrix M The observed data matrix D The degree matrix A The adjacency matrix L The Laplacian matrix Ψ The GBT matrix W The introduced auxiliary variable Z The Lagrange multiplier Suppose there are N sensor nodes in the WSNs. Using x i N i=1 denotes the sensor data, where x i ∈ R M represents a data vector collected by node i in time slot t 1 , t 2 , · · · , t m . The sample interval is assumed to be equal. Thus, the data matrix X ∈ R M×N can be used to represent the sensor data gathered by N sensor nodes in M time slots.
In order to reduce energy consumption in resource-constrained WSNs, only a small amount of sensor data is transmitted to the sink node. Let Ω ⊂ {1, 2, · · · , M} × {1, 2, · · · , N} denotes the indices of the corresponding observed data of X. Similarly, let Ω c denotes the indices of omitted value. Let π Ω Sensors 2016, 16, 1532 4 of 18 be the linear projection operator that keeps the entries in Ω invariant and adjusts the entries in Ω c to zero, that is: Suppose matrix M ∈ R M×N is the observed data, which is the incomplete version of matrix X with entries those outside Ω zeros. That is π Ω X ij = π Ω M ij , ∀ (i, j) ∈ Ω.
Our goal is to reduce the amount of data transmission to the sink node, and to design relevant matrix completion algorithm to reconstruct the original data matrix X as closely as possible. The observed ratio is defined as: Next, the features of the datasets are studied in detail, which would be utilized in our designed algorithms.

Exploring the Features of Datasets
In this subsection, three datasets are utilized for analysis. The first two datasets are gathered by the GreenOrbs [16] system, which is deployed in the forest environment with up to 330 nodes. Its topology is exhibited in Figure 1. We mainly consider the temperature and the humidity data collected between 3 and 5 August 2011. The third dataset is the smooth data generated with a second order Autoregressive (AR) model. In detail, the AR filter H(z) = 1 1+a 1 z −1 +a 2 z −2 is used, where a 1 and a 2 is predefined as −0.1 and −0.8 individually. Five hundred nodes assigned with the generated data are randomly deployed in a 1000 m × 1000 m area, which is shown in Figure 2. The detailed information about these datasets is given in Table 2.
is the observed data, which is the incomplete version of matrix X with entries those outside Ω zeros. That is Our goal is to reduce the amount of data transmission to the sink node, and to design relevant matrix completion algorithm to reconstruct the original data matrix X as closely as possible. The observed ratio is defined as: Next, the features of the datasets are studied in detail, which would be utilized in our designed algorithms.

Exploring the Features of Datasets
In this subsection, three datasets are utilized for analysis. The first two datasets are gathered by the GreenOrbs [16] system, which is deployed in the forest environment with up to 330 nodes. Its topology is exhibited in Figure 1. We mainly consider the temperature and the humidity data collected between 3 and 5 August 2011. The third dataset is the smooth data generated with a second order Autoregressive (AR) model. In detail, the AR filter and 2 a is predefined as −0.1 and −0.8 individually. Five hundred nodes assigned with the generated data are randomly deployed in a 1000 m × 1000 m area, which is shown in Figure 2. The detailed information about these datasets is given in Table 2.

Low-Rank Property
Since sensor readings have spatial and temporal correlations, the rank r of data matrix X would be small, such as . This low-rank property of the sensor data has been studied in previous papers [14,15,17]. Let   rank X denote the rank of matrix X. Candes et al. [18] proposed to solve the matrix completion problem by minimizing   rank X under suitable constraints.
However, the minimization problem of rank function cannot be figured out with a global solution in polynomial time because of its non-convexity. Fortunately, the nuclear norm  X , which can be solved using various convex programming methods, is the tightest convex relaxation of the rank function [12]. Also, the relationship between rank function and nuclear norm in matrix completion is similar to the relation between 0 l norm and the convex 1 l norm in compressive sensing.

GBT Sparsity
Since these datasets are coupled with their topologies, we construct the graph-based transform (GBT) to sparsely represent them. The sensor network is represented by a graph consists of a vertex set V (sensor nodes) and an edge set E (sensor links). The link ( , ) e i j is supposed to exist if the Euclidean distance between any two nodes (node i and node j ) is smaller than the communication range. The topological graph is mathematically denoted with the adjacent matrix A: The degree matrix D is a diagonal matrix, where the diagonal elements   . Thus, the Laplacian matrix can be induced as: Since matrix L is symmetric, the eigenvalue decomposition can be obtained. Then, the eigenvector of the Laplacian matrix L constitutes the columns of the graph-based transform matrix, which is denoted as Ψ. Clearly, the GBT matrix Ψ is orthogonal. Detailed information about GBT basis can be found in [19].
The performance of the GBT matrix as a sparse basis is demonstrated in Figure 3. As can be seen, nearly 10% of the sorted GBT coefficients assemble more the 99% of the energy. Therefore, the matrix

Low-Rank Property
Since sensor readings have spatial and temporal correlations, the rank r of data matrix X would be small, such as r min (M, N). This low-rank property of the sensor data has been studied in previous papers [14,15,17]. Let rank (X) denote the rank of matrix X. Candes et al. [18] proposed to solve the matrix completion problem by minimizing rank (X) under suitable constraints. However, the minimization problem of rank function cannot be figured out with a global solution in polynomial time because of its non-convexity. Fortunately, the nuclear norm X * , which can be solved using various convex programming methods, is the tightest convex relaxation of the rank function [12]. Also, the relationship between rank function and nuclear norm in matrix completion is similar to the relation between l 0 norm and the convex l 1 norm in compressive sensing.

GBT Sparsity
Since these datasets are coupled with their topologies, we construct the graph-based transform (GBT) to sparsely represent them. The sensor network is represented by a graph G = (V, E), which consists of a vertex set V (sensor nodes) and an edge set E (sensor links). The link e(i, j) is supposed to exist if the Euclidean distance between any two nodes (node i and node j) is smaller than the communication range. The topological graph is mathematically denoted with the adjacent matrix A: The degree matrix D is a diagonal matrix, where the diagonal elements D (i, i) denote the number of links connected to node , and D (i, j) = 0, ∀i = j. Thus, the Laplacian matrix can be induced as: Since matrix L is symmetric, the eigenvalue decomposition can be obtained. Then, the eigenvector of the Laplacian matrix L constitutes the columns of the graph-based transform matrix, which is denoted as Ψ. Clearly, the GBT matrix Ψ is orthogonal. Detailed information about GBT basis can be found in [19].
The performance of the GBT matrix as a sparse basis is demonstrated in Figure 3. As can be seen, nearly 10% of the sorted GBT coefficients assemble more the 99% of the energy. Therefore, the matrix Ψ −1 X is extremely sparse. In other words, the l 1 norm of matrix Ψ −1 X, represented as Ψ −1 X 1 , is very small.

The Proposed Optimization Algorithm
Previous matrix completion algorithm is not realistic, as the overfitting problem cannot be avoided when the sampling rate is low. Thus, the recovery error could be large due to the overfitting problem. To obtain exact reconstruction of the missing values, the GBT sparsity and the low-rank property of the data matrix are utilized. Finally, the following optimization problem is formulated: where  is the GBT sparsity regularization parameter.
ADMM [20] algorithm can blend the decomposability of dual ascent with an extra variable. Benefitted from the method of multiplier, the algorithm has superior convergence. With an introduced auxiliary variable , the above problem is rewritten as follows: In the following, some prerequisite properties are presented to sever as the foundation to solve the above problem.

Definition 1. Define the soft-thresholding operator is:
, , where 0   , and the operator can be applied to vectors or matrices in an element-wise manner.
In consideration of Definition 1, the following helpful theorem is given, as stated in [13].

The Proposed Optimization Algorithm
Previous matrix completion algorithm is not realistic, as the overfitting problem cannot be avoided when the sampling rate is low. Thus, the recovery error could be large due to the overfitting problem. To obtain exact reconstruction of the missing values, the GBT sparsity and the low-rank property of the data matrix are utilized. Finally, the following optimization problem is formulated: where λ is the GBT sparsity regularization parameter. ADMM [20] algorithm can blend the decomposability of dual ascent with an extra variable. Benefitted from the method of multiplier, the algorithm has superior convergence. With an introduced auxiliary variable W ∈ R M×N , the above problem is rewritten as follows: In the following, some prerequisite properties are presented to sever as the foundation to solve the above problem.

Definition 1.
Define the soft-thresholding operator is: where ε > 0, and the operator can be applied to vectors or matrices in an element-wise manner.
In consideration of Definition 1, the following helpful theorem is given, as stated in [13].
, where U ∈ R M×r and V ∈ R N×r are orthonormal matrix, and r = rank(Y) . Then, for any matrix Y ∈ R M×N and ∀λ ≥ 0, the following equations are available: and where UΣV T is the SVD of matrix Y.

Lemma 1.
Suppose Ψ ∈ R N×N is the unit orthogonal real matrix, and Ψ T denotes the transpose matrix of Ψ . Then, the Frobenius norm of any matrix A is invariant under a unitary transformation, that is: The inverse of matrix Ψ is equal to Ψ T . Then, following the definition of trace, we have: and Then, the GBT Regularization by Alternating Direction Method of Multipliers (GBTR-ADMM) is proposed to solve Problem (6). To get rid of the linear constraint in Problem (6), the augmented Lagrangian function is formulated as: where Z is the Lagrange multiplier and β > 0 is the penalty parameter. For a comparative analysis, a fixed parameter β and an adaptive update strategy for optimal β are all studied in the experimental analysis section. GBTR-ADMM updates the variables in a three-step iterative approach under the constraint of fixed β. The augmented Lagrangian function L (X, Z, W, β) is minimized in respect of the variables in a Gauss-Seidel manner. In each step, a single variable is updated by fixing the rest of the variables. By updating the variables alternately, each subproblem is solved with a closed form solution. More specifically, the iterations of GBTR-ADMM are formulated as follows: Firstly, the variable X k+1 is computed with fixed value of Z k and W k . Then L (X, Z, W, β) is minimized as follows: Removing the constant term in Function (14), this function can be rewritten as: Obviously, the above optimization problem has the same form as defined in Theorem 1. Thus, a closed form solution is obtained as follows: where U, V , and Σ P are obtained from the SVD of matrix P, and P equals to W k − 1 β Z k . Secondly, the variable W k+1 is updated in the choice of default values X k+1 and Z k . The minimization of L (X, Z, W, β) goes as follows: Ignoring the constant term in this step, we can obtain the following optimization problem: Taking into consideration of the orthogonal invariance of the Forbenius norm, which is defined in Lemma 1, we obtain the following theorem. (18) is defined as follows:

Theorem 2. The closed form solution of Problem
Proof. Since matrix Ψ −1 is orthogonal, the following equation is obvious in combination with Lemma 1.
and defining Q = Ψ −1 W, we have: By Theorem 1, the closed form solution of the above Problem (21) is obtained as follows: Substituting Q = Ψ −1 W to Problem (22), then the following closed form solution is available: In view of the second constraint term in Problem (6), the final form of W k+1 is defined as: Thirdly, with the derived value of X k+1 and W k+1 in the above two steps, the calculation of Lagrange multiplier is updated as: The main procedure of GBTR-ADMM is shown in Algorithm 1. Note that the choice of penalty parameter β has high influence on the performance of ADMM algorithm. As it is difficult to choose an optimal value, the adaptive renewal mechanism is preferred in practical application. The performance difference of GBTR-ADMM with varying penalty value is studied in Section 7.1. What is more, the convergence of the ADMM based method is theoretically demonstrated in [20].

The Proposed Method for Accelerated Convergence
The performance of ADMM is highly sensitive to the number of variables and the number of constraints. As is stated in [21,22], more memory is required, and the rate of convergence is reduced with multiple variables constraints. What is more, the convergence property is not proven theoretically when the number of variables is greater than or equal to 3. In optimization Problem (6), these two constraints are considered separately, as shown in Algorithm 1. This may slow down the convergence speed.
In this section, a new approach is proposed to solve Problem (6) with fast constringency speed. Firstly, the two constraints in Problem (6) are merged together in a linear operator. Thus, the convergence rate is accelerated with only one constraint. Then, we introduce the GBT Regularization by accelerated alternating direction method of multipliers (GBTR-A2DM2) As we know, the convergence rate of A2DM2 [23,24] algorithm is O 1 k 2 while the convergence rate of ADMM (as Algorithm 1) is O 1 k .

The Fusion of Two Constraints
In consideration of the two constraints in Problem (6), X = W and π Ω (W ) = π Ω (M), two linear operators, which are represented as A and B: R M×N → R 2M×2N , are defined as follows: where C ∈ R 2M×2N is a constant matrix. Thus, Problem (6) is reformulated as follows: Also, the Lagrange function for the above optimization problem is: Similar to Algorithm 1, GBTR-A2DM2 decomposes the minimization of L (X, Z, W, β) into several subproblems. In each subproblem, GBTR-A2DM2 updates a variable keeping in mind that the other variables are fixed. Specifically, the optimization scheme of GBTR-A2DM2 for Problem (28) is resolved in the following steps: The pseudocode of GBTR-A2DM2 algorithm is shown in Algorithm 2. Next, we will discuss the accelerated technique of GBTR-A2DM2 with a restarting rule.

The Accelerated Technique
Since the objective function in Problem (27) is not very convex, the accelerated ADMM method with a restart rule is employed. To determine when to restart the value assignment, the primal error and the dual error are combined: whereẐ k andŴ k represent the second updated step in iteration steps 6-7 of Algorithm 2. For each iteration, m k is compared with m k−1 and if m k < ηm k−1 , where η is defined equal to 0.999, the algorithm is accelerated with steps 5-7. Otherwise, the method is restarted in process of steps 8-9. In comparison with GBTR-ADMM, Algorithm 2 has a higher convergence rate. Also, the convergence property of Algorithm 2 is guaranteed by A2DM2 with a restarting rule [23].

Time Complexity Analysis
In this part, the computational complexity of the proposed algorithms is discussed. The calculation of an inverse matrix cost much, which has the time complexity of O(n 3 ) (n is the dimension of an invertible matrix). Since matrix Ψ is orthogonal, the expensive computation of matrix inversion in our implementation can be substituted by its transposition. Thus, the dominated computational cost of GBTR-ADMM and GBTR-A2DM2 is the execution of matrix SVD in each iteration. As pointed out in [25], the time complexity of SVD operation is O(MN 2 ). In our implementation, the famous PROPACK [26] is utilized to perform partial SVD for the proposed algorithms. Since the low-rank property of the objective matrix, it is inefficient to compute the full SVD. To obtain the dominated energy of the objective matrix, only those singular values exceeding than a certain threshold are necessary. The limitation of PROPACK is that it cannot automatically determine the necessary calculations, except for a predefined number. Thus we are supposed to estimate the number of singular values and assign the number to PROPACK in each iteration.
Suppose svp k is the number of positive singular values of X k , and sv k is the number of singular value to be measured at k-th iteration. Then, the following updated strategy [27] is used, where the initial estimated value of sv 0 is 10. Benefiting from the software package PROPACK, the time complexity for a M × N matrix with rank of r is O(rMN). Hence, the total time complexity of our proposed algorithms is O(rMN). Nevertheless, the state of art algorithms for matrix completion problem [15,17] demand a complexity of O(r 2 MN) for each iteration.

Performance Evaluation
In this section, we evaluate the performances of GBTR-ADMM and GBTR-A2DM2. The experimental datasets and their topological structures are described in Section 3. Since the proposed algorithms are heavily influenced by several input parameters, it is necessary to choose the optimal parameters to maximize the algorithm performance. With the optimal parameters for GBTR-ADMM and GBTR-A2DM2, the recovery accuracy and the convergence properties are compared with the state of art algorithms. At last, the energy consumption of the proposed algorithms are compared with the state of art data gathering methods for WSNs. Simulation results show that GBTR-ADMM and GBTR-A2DM2 can highly reduce energy consumption in WSNs. Thus, the network lifetime is prolonged.
To measure the performance of the proposed algorithm, the reconstructed data matrixX is achieved. Thus, the recovery performance is estimated by the Normalized Mean Absolute Error (NMAE):

Parameter Setting
In this subsection, the choice of optimal parameters for GBTR-ADMM is discussed. Previous studies have shown that global convergence for ADMM algorithm holds for any fixed β > 0. However, different parameter values result in various convergence speeds. Thus, the input values, as listed in Table 3, to Algorithm 1 are selected by experience to obtain the best performance. The performance of GBTR-ADMM is also studied with different parameter values of β. The variation of the objective function values of Problem (13) with the increase of iteration numbers is shown in Figure 4. As we can see, β = β 0 achieves the best performance. Meanwhile, the descending speed becomes slower when the choice of parameter β is too large or too small. This is because the penalty parameter β trades off between minimizing the primal residual and the residual of the dual problem. A large penalty value may drop the primal residual, but at the expense of an increase of the dual residual, and vice versa. Table 3. Input values to Algorithm 1.   Figure 4 demonstrates the results in the synthesized datasets, and the optimal chosen value of β changes randomly in other datasets. Instead, an adaptive penalty update method is preferred, which is based on previous study [20]. The update strategy is formulated as follows:

Parameter Name
where βmax denotes the maximum value of the βk. The value of variable ρ is updated as: where 0 0   is a constant and  is the predefined threshold value. Obviously, when the residual value between is less than the threshold, the value of 1 k   increases to 0 k   . Thus, the convergence speed is improved in this way.
The effect of the sparsity regularization parameter λ is also analyzed. Figure 5 shows the variation of recovery error with different parameter value. As can be seen, the recovery error is quite large with small value of λ. With the increase of λ, recovery error declines rapidly, and remains stable as λ > 0.01. So, the optimal value for the sparsity regularization parameter λ is set as 0.01 in our experiments. Since similar trends are obtained for GBTR-A2DM2, we just omit it here.  Figure 4 demonstrates the results in the synthesized datasets, and the optimal chosen value of β changes randomly in other datasets. Instead, an adaptive penalty update method is preferred, which is based on previous study [20]. The update strategy is formulated as follows: where β max denotes the maximum value of the β k . The value of variable ρ is updated as: where ρ 0 > 0 is a constant and κ is the predefined threshold value. Obviously, when the residual value between X k+1 − X k 2 F and W k+1 − W k 2 F is less than the threshold, the value of β k+1 increases to ρ 0 β k . Thus, the convergence speed is improved in this way.
The effect of the sparsity regularization parameter λ is also analyzed. Figure 5 shows the variation of recovery error with different parameter value. As can be seen, the recovery error is quite large with small value of λ. With the increase of λ, recovery error declines rapidly, and remains stable as λ > 0.01. So, the optimal value for the sparsity regularization parameter λ is set as 0.01 in our experiments. Since similar trends are obtained for GBTR-A2DM2, we just omit it here. to 0 k   . Thus, the convergence speed is improved in this way.
The effect of the sparsity regularization parameter λ is also analyzed. Figure 5 shows the variation of recovery error with different parameter value. As can be seen, the recovery error is quite large with small value of λ. With the increase of λ, recovery error declines rapidly, and remains stable as λ > 0.01. So, the optimal value for the sparsity regularization parameter λ is set as 0.01 in our experiments. Since similar trends are obtained for GBTR-A2DM2, we just omit it here.

Recovery Accuracy
In this subsection, we compare the recovery accuracy of GBTR-ADMM and GBTR-A2DM2 with the state of art algorithms for matrix completion. The first chosen method is the Spatio-Temporal Compressive Data Collection (STCDG) [15]. The second method is the Compressive Data Collection (CDC) [11]. GBTR-ADMM is considered to solve the optimization Problem (6), while GBTR-A2DM2 is used for the optimization Problem (27) with only one constraint.
Simulations are executed on both the real datasets and the synthesized dataset, which are exploited in detail in Section 3. For each parameter setting of the simulation, the results are averaged for 50 independent trials. Figures 6-8 show that our GBTR based methods can reconstruct the missing values with high accuracy. In general, the recovery errors of all reconstruction algorithms decrease rapidly with the increase of the sampling ratio. When the sampling ratio is high enough, all reconstruction methods achieve smaller recovery error. Since our proposed two GBTR based methods are used to solve the same matrix completion problems, their performances are nearly the same. Figure 6 shows the recovery errors on GreenOrbs temperature data. As can been seen, our proposed methods achieve about 25% recovery error while the error of other two algorithms is more than 80%, when the sampling ratio is 1%.

Recovery Accuracy
In this subsection, we compare the recovery accuracy of GBTR-ADMM and GBTR-A2DM2 with the state of art algorithms for matrix completion. The first chosen method is the Spatio-Temporal Compressive Data Collection (STCDG) [15]. The second method is the Compressive Data Collection (CDC) [11]. GBTR-ADMM is considered to solve the optimization Problem (6), while GBTR-A2DM2 is used for the optimization Problem (27) with only one constraint.
Simulations are executed on both the real datasets and the synthesized dataset, which are exploited in detail in Section 3. For each parameter setting of the simulation, the results are averaged for 50 independent trials. Figures 6-8 show that our GBTR based methods can reconstruct the missing values with high accuracy. In general, the recovery errors of all reconstruction algorithms decrease rapidly with the increase of the sampling ratio. When the sampling ratio is high enough, all reconstruction methods achieve smaller recovery error. Since our proposed two GBTR based methods are used to solve the same matrix completion problems, their performances are nearly the same. Figure 6 shows the recovery errors on GreenOrbs temperature data. As can been seen, our proposed methods achieve about 25% recovery error while the error of other two algorithms is more than 80%, when the sampling ratio is 1%. Similar results can be obtained from Figure 7, which is simulated on the humidity dataset. The recovery error of GBTR based methods is still much less than STCDG and CDC when the sampling ratio is small. When the sampling is 1%, GBTR-ADMM and GBTR-A2DM2 can reconstruct the original missing values with recovery errors of less than 20%. Meanwhile, the recovery errors of STCDG and CDC are nearly 100%.
In Figure 8, the experiment results show that GBTR based methods outperform STCDG and Similar results can be obtained from Figure 7, which is simulated on the humidity dataset. The recovery error of GBTR based methods is still much less than STCDG and CDC when the sampling ratio is small. When the sampling is 1%, GBTR-ADMM and GBTR-A2DM2 can reconstruct the original missing values with recovery errors of less than 20%. Meanwhile, the recovery errors of STCDG and CDC are nearly 100%.
In Figure 8, the experiment results show that GBTR based methods outperform STCDG and CDC by a larger margin. Compared with Figures 6 and 7, Figure 8 shows that GBTR-ADMM and GBTR-A2DM2 achieve the best performance on the synthesized dataset. The recovery error on synthesized dataset is smaller than the real datasets at the same sampling ratio. The reason is that the synthesized data have much better sparsity than the other two real datasets under the GBT basis.

Convergence Behavior Analysis
In this subsection, the convergence performances are studied in the synthesized dataset. The compared methods are SVD and STCDG. For each method, we set the same stop conditions, where the tolerance error  is 4 10  . Figure 9 shows the necessary number of iterations to obtain accurate reconstruction at different sampling ratios. As can be seen from Figure 9, the convergence speed of our proposed two methods surpasses the SVD and STCDG. SVD has the slowest convergence rate of the four methods. Also, STCDG converges faster than SVD. Although the recovery accuracy of GBTR-A2DM2 and GBTR-ADMM behaves similar, as shown in Figures 6-8, GBTR-A2DM2 converges much faster than GBTR-ADMM. Note that GBTR-ADMM needs nearly 160 iterations to converge at the sampling ratio 0.9, while only about 40 iterations leads to the convergence of GBTR-A2DM2. Even when the sampling ratio is 0.5, the necessary number of iterations for GBTR-A2DM2 is about 125,

Convergence Behavior Analysis
In this subsection, the convergence performances are studied in the synthesized dataset. The compared methods are SVD and STCDG. For each method, we set the same stop conditions, where the tolerance error  is 4 10  . Figure 9 shows the necessary number of iterations to obtain accurate reconstruction at different sampling ratios. As can be seen from Figure 9, the convergence speed of our proposed two methods surpasses the SVD and STCDG. SVD has the slowest convergence rate of the four methods. Also, STCDG converges faster than SVD. Although the recovery accuracy of GBTR-A2DM2 and GBTR-ADMM behaves similar, as shown in Figures 6-8, GBTR-A2DM2 converges much faster than GBTR-ADMM. Note that GBTR-ADMM needs nearly 160 iterations to converge at the sampling ratio 0.9, while only about 40 iterations leads to the convergence of GBTR-A2DM2. Even when the sampling ratio is 0.5, the necessary number of iterations for GBTR-A2DM2 is about 125,

Convergence Behavior Analysis
In this subsection, the convergence performances are studied in the synthesized dataset. The compared methods are SVD and STCDG. For each method, we set the same stop conditions, where the tolerance error ξ is 10 −4 . Figure 9 shows the necessary number of iterations to obtain accurate reconstruction at different sampling ratios. As can be seen from Figure 9, the convergence speed of our proposed two methods surpasses the SVD and STCDG. SVD has the slowest convergence rate of the four methods. Also, STCDG converges faster than SVD. Although the recovery accuracy of GBTR-A2DM2 and GBTR-ADMM behaves similar, as shown in Figures 6-8, GBTR-A2DM2 converges much faster than GBTR-ADMM. Note that GBTR-ADMM needs nearly 160 iterations to converge at the sampling ratio 0.9, while only about 40 iterations leads to the convergence of GBTR-A2DM2. Even when the sampling ratio is 0.5, the necessary number of iterations for GBTR-A2DM2 is about 125, which is less than half of GBTR-ADMM. Next, when the sampling ratio is fixed as 0.6, the relative recovery errors of all the compared methods are analyzed. As can be seen from Figure 10, Compared to SVD and STCDG, both GBTR-ADMM and GBTR-A2DM2 gain less error in several iteration processes. Clearly, GBTR-A2DM2 converges much faster, which converges in about 100 iterations. Also, with no more than 250 iterations, both GBTR-ADMM and GBTR-A2DM2 terminate as the relative recovery errors drop below the tolerance error. Meanwhile, the relative errors of SVD and STCDG are about one order of magnitude larger than GBTR-ADMM and GBTR-A2DM2. In general, compared with the state of the art methods for matrix completion, our proposed methods can achieve smaller recovery error at the same number of iterations.

Energy Consumption and Network Lifetime
In this subsection, the energy consumption of the proposed algorithms for data gathering problems is analyzed. Five nodes are randomly deployed in a 1000 m × 1000 m area. The topology is shown in Figure 2, where the sink node is deployed in the center. The data transmission and recovery process are fulfilled in three steps. First, the sink node broadcasts the sampling ratio through the whole network. In our setting, the sampling ratio determines the probability of node to gather data. In the second step, the selected nodes transmit the gathering data to the sink node. Finally, these Next, when the sampling ratio is fixed as 0.6, the relative recovery errors of all the compared methods are analyzed. As can be seen from Figure 10, Compared to SVD and STCDG, both GBTR-ADMM and GBTR-A2DM2 gain less error in several iteration processes. Clearly, GBTR-A2DM2 converges much faster, which converges in about 100 iterations. Also, with no more than 250 iterations, both GBTR-ADMM and GBTR-A2DM2 terminate as the relative recovery errors drop below the tolerance error. Meanwhile, the relative errors of SVD and STCDG are about one order of magnitude larger than GBTR-ADMM and GBTR-A2DM2. In general, compared with the state of the art methods for matrix completion, our proposed methods can achieve smaller recovery error at the same number of iterations. Next, when the sampling ratio is fixed as 0.6, the relative recovery errors of all the compared methods are analyzed. As can be seen from Figure 10, Compared to SVD and STCDG, both GBTR-ADMM and GBTR-A2DM2 gain less error in several iteration processes. Clearly, GBTR-A2DM2 converges much faster, which converges in about 100 iterations. Also, with no more than 250 iterations, both GBTR-ADMM and GBTR-A2DM2 terminate as the relative recovery errors drop below the tolerance error. Meanwhile, the relative errors of SVD and STCDG are about one order of magnitude larger than GBTR-ADMM and GBTR-A2DM2. In general, compared with the state of the art methods for matrix completion, our proposed methods can achieve smaller recovery error at the same number of iterations.

Energy Consumption and Network Lifetime
In this subsection, the energy consumption of the proposed algorithms for data gathering problems is analyzed. Five nodes are randomly deployed in a 1000 m × 1000 m area. The topology is shown in Figure 2, where the sink node is deployed in the center. The data transmission and recovery process are fulfilled in three steps. First, the sink node broadcasts the sampling ratio through the whole network. In our setting, the sampling ratio determines the probability of node to gather data. In the second step, the selected nodes transmit the gathering data to the sink node. Finally, these

Energy Consumption and Network Lifetime
In this subsection, the energy consumption of the proposed algorithms for data gathering problems is analyzed. Five nodes are randomly deployed in a 1000 m × 1000 m area. The topology is shown in Figure 2, where the sink node is deployed in the center. The data transmission and recovery process are fulfilled in three steps. First, the sink node broadcasts the sampling ratio through the whole network. In our setting, the sampling ratio determines the probability of node to gather data. In the second step, the selected nodes transmit the gathering data to the sink node. Finally, these missing data are reconstructed by implementing GBTR-ADMM and GBTR-A2DM2 at the sink node. The compared methods are CDC and STCDG. In the traditional data gathering method, all sensor nodes are required to transmit their sampling data to the sink node. Thus, the traditional method is selected as a baseline method for comparison. The energy consumption model in paper [28] is employed in our simulation. The detailed simulation parameters are presented in Table 4. The initial energy for every sensor node is 2 J. Each packet contains 64 bits. The network is supposed to be symmetrical. The energy consumption for one bit transmission, defined as E Tx , is 100 nJ. Meanwhile, the reception of a packet consumes E Rx = 120 nJ. E Amp is the unit energy consumption of the power amplifying circuit. The synthesized dataset is exploited in our experiment. The network lifetimes of CDC, STCDG, and our proposed methods are evaluated. Detailed information is revealed in Figure 11. Note that the total energy consumption of the baseline method is relatively constant. That is because the baseline method transmits all the sensor data to the sink node no matter how the sampling ratio varies. To be different, the total power consumptions of CDC, STCDG, and GBTR based methods increase with enlargement of the sampling ratio. The reason is that more sensor data are needed to be transmitted when the sampling ratio increases. Thus, the network lifetime decreases. GBTR-ADMM and GBTR-A2DM2 outperform CDC at the same setting value of sampling ratio. However, the energy consumption of our proposed methods is equal to the baseline method when the sampling radio is exactly 1. Note that the lifetime of CDC is smaller than the baseline when the sampling ratios are higher than 75%. This phenomenon could be explained as per below. In CDC, all sensor nodes transmit the data M times, which is the necessary number to reconstruct original signals. At the same time, the ordinary nodes need only one data transmission and necessary relay transmissions for each sampling in the baseline method. Thus, when the sampling ratio is higher than a specific threshold, the total numbers of transmission for CDC is larger than that in the baseline method. In addition, since STCDG and our proposed methods are all based on the matrix completion theory, the curve variations over their lifetime coincide with each other.
transmit the data M times, which is the necessary number to reconstruct original signals. At the same time, the ordinary nodes need only one data transmission and necessary relay transmissions for each sampling in the baseline method. Thus, when the sampling ratio is higher than a specific threshold, the total numbers of transmission for CDC is larger than that in the baseline method. In addition, since STCDG and our proposed methods are all based on the matrix completion theory, the curve variations over their lifetime coincide with each other. Figure 11. Network lifetime comparison. Figure 11. Network lifetime comparison.

Conclusions and Future Works
In this paper, the data gathering problem based on Matrix Completion theory is studied. Except for the low-rank property, the sensed data are observed to be sparse under the graph based transform. By taking full advantage of these features, two novel reconstruction algorithms (named GBTR-ADMM and GBTR-A2DM2) are proposed. The time complexity is also analyzed, which shows their complexity is low. Several experiments on both real datasets and synthesized datasets are carried out. The experiment results show that our proposed methods outperform the state of the art algorithm for data gathering problems in WSNs Furthermore, it is observed that GBTR-A2DM2 converges much faster than GBTR-ADMM. For future works, will focus on applying our proposed algorithms to other datasets of real networks, which may exhibit complex topological information other than random networks, such as the scale-free or the small-world networks.