Multi-Matrices Factorization with Application to Missing Sensor Data Imputation

Abstract

We formulate a multi-matrices factorization model (MMF) for the missing sensor data estimation problem. The estimation problem is adequately transformed into a matrix completion one. With MMF, an n-by-t real matrix, R, is adopted to represent the data collected by mobile sensors from n areas at the time, T₁, T₂, … , T_t, where the entry, R_i,j, is the aggregate value of the data collected in the ith area at T_j. We propose to approximate R by seeking a family of d-by-n probabilistic spatial feature matrices, U₍₁₎, U₍₂₎, … , U_(t), and a probabilistic temporal feature matrix, V ∈ ℝ^d×t, where $R_j\approx U_{(j)}^T{T}_j$. We also present a solution algorithm to the proposed model. We evaluate MMF with synthetic data and a real-world sensor dataset extensively. Experimental results demonstrate that our approach outperforms the state-of-the-art comparison algorithms.

1. Introduction

In this work, we study the following missing sensor data imputation problem: Let the matrix, R ∈ ℝ^n×t, consist of the data collected by a set of mobile sensors in spacial areas S₁, S₂, … S_n at time points T₁ < T₂ < … < T_t, where the entry, R_i,j, is the aggregate value collected by the sensors in S_i at T_j. In particular, if there is no sensor in S_i at time T_j, we denote the value of R_i,j as “ ⊥ ”, which indicates that it is missing. Our focus is to find the suitable estimations for the missing values in a given incomplete matrix, R. Results of this research could be helpful in recovering missing values in statistical analyses. For example, to predict floods, people usually place geographically distributed sensors in the water to continuously monitor the rising water levels. However, some data in a critical period of time might be corrupted, due to, e.g., sensor hardware failures. Such a kind of data needs be recovered to guarantee the prediction accuracy.

Many efforts have been devoted to the missing sensor data imputation problem. Typical examples include k nearest neighbor-based imputation [1], multiple imputation [2], hot/cold imputation [3], maximum likelihood and Bayesian estimation [4] and expectation maximization [5]. However, despite the various implementations of these methods, their main essence is based on the local consistency of the sensor data, i.e., the data collected at adjacent time points within the same spacial area should be close to each other, as well as the data collected at the same time from neighboring areas. We refer to them as local models. As is well known, these local models suffer from the cumulative error problem in scenarios where the missing ratio is high.

Matrix factorization (MF), as a global model, has caught substantial attention in recent years. Typically, in the Netflix rating matrix completion competition [6], some variations of the MF model, e.g., [7,8], achieved state-of-the-art performances, showing their potential to recover the missing data from highly incomplete matrices. On the other side, many well-studied MF models, such as non-negative matrix factorization [9], max margin matrix factorization [10,11], and probabilistic matrix factorization [7], are based on the i.i.d.assumption [12], which, in terms of our problem, implies that the neighborhood information among the data is disregarded and, hence, leaves vast room for improvement.

We in this paper, we propose a multi-matrices factorization model (MMF), which can be outlined as follows. For a matrix, X, denote X_j the jth column of X. Given a sensor data matrix, R, we seek a set of matrices, U₍₁₎, U₍₂₎, …, U_(t) ∈ ℝ^d×n, and a matrix, V ∈ ℝ ^d×t, such that for i = 1, 2, …,t, $U_{(i)}^T{V}_i\approx R_i$. Here, U_(i) is referred to as the spatial feature matrix, in which the jth column, U_(i),j, is the feature vector of area S_j at T_i. Similarly, V is referred to as the temporal feature matrix, in which the jth column, V_j, is the temporal feature vector of T_j. To predict the missing values in R, we first fit the matrices, U₍₁₎, U₍₂₎, …, U_(t) (single sub-indexes of matrix mean columns) and V, with the non-missing values in R; then, for each R_i,j = “ ⊥ ”, we take its estimation as ${\widehat{R}}_{i,j}=U_{(j),i}^T{V}_j$.

The remainder of the paper is organized as follows: Section 2 summarizes the notations used in the paper. Section 3 studies the related work on matrix factorization. In Section 4, we present our multi-matrices factorization model. The algorithm to solve the proposed model is outlined in Section 5. Section 6 is devoted to the experimental evaluations. Finally, our conclusions are presented in Section 7, followed by a presentation of future work, acknowledgments and a list of references.

2. Notations

For a vector V = [v₁, v₂, …v_n]′ ∈ ℝⁿ, we use ‖V‖₀, ‖V‖₁ and ‖V‖₂ to denote its 0-Norm, 1-Norm and 2-Norm, respectively, as follows:

-

${\Vert V\Vert }_0={\sum }_{i=1}^n\mathbb{\text{I}}({\mathit{\text{v}}}_i\ne 0)$, i.e., the number of nonzero entries of V;

-

${\Vert V\Vert }_1={\sum }_{i=1}^n\left|v_i\right|$;

-

${\Vert V\Vert }_2=\sqrt{{\sum }_{i=1}^n{v}_i^2}$.

For a matrix, X ∈ ℝ^n×m, we denote its Frobenius norm as ${\Vert X\Vert }_F=\sqrt{{\sum }_{i=1}^n{\sum }_{j=1}^m{X}_{i,j}^2}$.

3. Matrix Factorization

The essence of the Matrix Factorization problem is to find two factor matrices, U and V, such that their product can approximate the given matrix, R, i.e., R ≈ U^TV. As a fundamental model of machine learning and data mining, the MF method has achieved state-of-the-art performance in various applications, such as collaborative filtering [13], text analysis [14], image analysis [9,15] and biology analysis [16]. In principle, for a given matrix, R, the MF problem can be formulated as the optimization model below:

$$\left\{U^{*},V^{*}\right\}=\underset{U,V}{\text{min}}\mathit{\text{Loss}}(U^{T}V,R)$$

(1)

where the lose function, Loss, is used to measure the closeness of the approximation, U^TV, to the target, R. Usually, Loss(U^TV, R) can be decomposed into the sum of the pairwise loss between the entries of U^TV and R; that is, $\mathit{\text{Loss}}(U^{T}V,R)={\sum }_{i=1}^n{\sum }_{j=1}^m\mathit{\text{loss}}({(U^{T}V)}_{i,j},R_{i,j})$. Some of the most used forms include the square loss (loss(x, y) = (x - y)²) [7,8,17], the 0-1 loss (loss(x, y) = 𝕀(x = y)) [11] and the divergence loss $(\mathit{\text{loss}}(x,y)=\mathit{\text{xlog}}\frac{x}{y}-x+y)$ [9].

It is notable that for Equation (1), if {U*, V*} is a solution to it, then for any scalar, κ > 0, $\{\kappa U^{*},\frac{1}{\kappa }{V}^{*}\}$ is also another solution; hence, problem (1) is ill-posed. To overcome this obstacle, various constraints on U and V are introduced, such as constraints on the entries [15], constraints on the sparseness [18,19], constraints on the norms [7,20] and constraints on the ranks [21,22]. All these constraints, from the perspective of the statistical learning theory, can be regarded as the length of the model to be fitted. According to the minimum description length principle [23,24], a smaller length means a better model; hence, most of them can be incorporated into Model (1) as the additional regularized terms, that is:

$$\left\{U^{*},V^{*}\right\}=\underset{U,V}{\mathit{\text{min}}}\mathit{\text{Loss}}(U^{T}V,R)+P(U,V)$$

(2)

where the regularization factor, P(U, V), corresponds to the constraints on U and V .

As a transductive model, Model (2) has many nice mathematical properties, such as the generalization error bound [10] and the exactness [17,25]. However, as is well known, when compared with the generative model, one of the main restrictions of the transductive model is that it can hardly be used to describe the relations existing in the data. In particular, in terms of our problem, even though Model (2) may work well, it is laborious to express the dynamics of the data over time.

4. The Proposed Model

In this section, we elaborate on our multi-matrices factorization (MMF) approach. Given the sensor data matrix, R, in which the entry, R_i,j (1 ≤ i ≤ n and 1 ≤ j ≤ t), is collected from S_i at T_j, our goal is to find the factor matrices, U₍₁₎, U₍₂₎, …, U_(t) ∈ ℝ^d×n and V ∈ ℝ^d×t, such that: for j = 1, 2, …, t,

$$R_j\approx U_{(j)}^T{V}_j$$

(3)

where U_(j) is regarded to be composed of the spatial features of all areas at T_j and V is treated as consisting of the temporal features of all time points. We denote the i th column of U_j as U_(j),_i, which corresponds to the spatial feature value of S_i at T_j, and denote the j th column of V as V_j, which corresponds to the temporal feature value of T_j.

Taking advantage of the knowledge of the probability graph model, we assume that the dependent structure of the data in U₍₁₎, U₍₂₎, …, U_(t), V and R is as illustrated in Figure 1. More specifically, we have the following assumptions:

• Columns of U_(j) (1 ≤ j ≤ t) are linearly independent, i.e.,

  $$Pr(U_{(j)})=\prod _{i=1}^{n}Pr(U_{(j),i})$$

• U_(1),i (1 ≤ i ≤ n) follows the same Gaussian distribution with a mean of zero and a covariance matrix ${\sigma }_U^2$, i.e.,

  $$\mathit{\text{Pr}}(U_{(1),i}\mid {\sigma }_U)={(2\pi {\sigma }_U^2)}^{-\frac{d}{2}}\mathit{\text{exp}}\{-\frac{{\Vert U_{(1),i}\Vert }_2^2}{2{\sigma }_U^2}\}$$

• U_(j),i (1 ≤ i ≤ n) are dependent in time order with the pre-specified priors, ζ_U and σ_U, i.e.,

  $$\mathit{\text{Pr}}(U_{(j),i}\mid {\zeta }_U,{\sigma }_U)=\mathit{\text{Pr}}(U_{(1),i}\mid {\zeta }_U,{\sigma }_U)\times \prod _{j=2}^t\mathit{\text{Pr}}(U_{(j),i}\mid U_{(j-1),i},{\zeta }_U,{\sigma }_U)$$

  Moreover, for j > 1, we assume U_(j),i is a Laplace random vector with location parameter U_(j-1),i and scale parameter ζ_U, namely:

  $$\mathit{\text{Pr}}(U_{(j),i}\mid U_{(j-1),i},{\zeta }_U,{\sigma }_U)=\frac{1}{2{\zeta }_U}\text{exp}\{-\frac{\left|U_{(j),i}-U_{(j-1),i}\right|}{{\zeta }_U}\}$$

• The columns of V are linearly dependent in time order with the pre-specified priors, ζ_V and σ_V, i.e.,

  $$\mathit{\text{Pr}}(V\mid {\zeta }_V,{\sigma }_V)=\mathit{\text{Pr}}(V_1\mid {\zeta }_V,{\sigma }_V)\times \prod _{j=2}^t\mathit{\text{Pr}}(V_j\mid V_{j-1},{\zeta }_V,{\sigma }_V)$$

  We also assume that, for j > 1:

  $$\mathit{\text{Pr}}(V_j\mid V_{j-1},{\zeta }_V,{\sigma }_V)=\frac{1}{2{\zeta }_V}\mathit{\text{exp}}\{-\frac{\left|V_j-V_{j-1}\right|}{{\zeta }_V}\}$$

• The (i, j)th entry of R (1 ≤ i ≤ n, 1 ≤ j ≤ t) follows Gaussian distribution with a mean of $U_{(j),i}^T{V}_j$ and variance ${\sigma }_R^2$, i.e.,

  $$\mathit{\text{Pr}}(R_{i,j}\mid U_{(j),i}^T{V}_j,{\sigma }_R^2)={(2\pi {\sigma }_R^2)}^{-\frac{1}{2}}\mathit{\text{exp}}\{-\frac{{(R_{i,j}-U_{(j),i}^T{V}_i)}^2}{2{\sigma }_R^2}\}$$

Now, given R and the priors, σ_U, σ_V, σ_R, ζ_U, and ζ_V, let U = {U₍₁₎, U₍₂₎, …, U_(t)}; below, we find a solution to the following equation:

$$\{U^{*},V^{*}\}=\underset{U,V}{\text{arg}}\mathit{\text{maxPr}}(U,V\mid R,{\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)$$

(4)

First, applying Bayes' theorem, we have:

$$\mathit{\text{Pr}}(U,V\mid R,{\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)=\frac{\mathit{\text{Pr}}(U,V,R\mid {\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)}{\mathit{\text{Pr}}(R\mid {\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)}$$

Since R is observed and σ_U, σ_V, σ_R, ζ_U and ζ_V are pre-specified, the denominator, Pr(R∣σ_U, σ_V, σ_R, ζ_U, ζ_V), can be treated as a constant. Therefore:

$$\begin{array}{l}\text{Equation}(4)\iff \\ \{U^{*},V^{*}\}=\underset{U,V}{arg}\mathit{\text{maxPr}}(U,V,R\mid {\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)\end{array}$$

(5)

Combing Assumptions (I.) ∼ (V.) and the dependency structure illustrated in Figure 1, we have:

$$\begin{array}{l}\mathit{\text{Pr}}(R,U,\text{V}\mid {\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)\\ =\mathit{\text{Pr}}(R\mid U,V,{\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)\times \mathit{\text{Pr}}(U,V\mid {\sigma }_U,{\sigma }_V,{\sigma }_R,{\zeta }_U,{\zeta }_V)\\ =\mathit{\text{Pr}}(R\mid U,V,{\sigma }_R)\times \mathit{\text{Pr}}(U\mid {\sigma }_U,{\zeta }_U)\times \mathit{\text{Pr}}(V\mid {\sigma }_V,{\zeta }_V)\\ =\text{∏}_{i=1}^n\text{∏}_{j=1}^t\mathit{\text{Pr}}(R_{i,j}\mid U_{(j),i},V_j,{\sigma }_R)\times \text{∏}_{j=1}^t\mathit{\text{Pr}}(U_{(j)}\mid {\sigma }_U,{\zeta }_U)\times \mathit{\text{Pr}}(V\mid {\sigma }_V,{\zeta }_V)\\ =\text{∏}_{i=1}^n\text{∏}_{j=1}^t\mathit{\text{Pr}}(R_{i,j}\mid U_{(j),i},V_j,{\sigma }_R)\times \text{∏}_{j=1}^t\mathit{\text{Pr}}(U_{(j),1}\mid {\sigma }_U)\text{∏}_{i=1}^n\text{∏}_{j=2}^t\mathit{\text{Pr}}(U_{(j),i}\mid U_{(j-1),i},{\zeta }_U)\times \mathit{\text{Pr}}(V_1\mid {\sigma }_V)\times \text{∏}_{j=2}^t\mathit{\text{Pr}}(V_j\mid V_{j-1},{\zeta }_V)\\ \propto exp(-\frac{1}{2{\sigma }_R^2}\text{∑}_{i=1}^n\text{∑}_{j=1}^t{(U_{(j),i}^T{V}_i-R_{i,j})}^2))\times exp(-\frac{1}{2{\sigma }_U^2}\text{∑}_{i=1}^n{\Vert U_{(i),1}\Vert }_2^2)\times \text{∏}_{i=1}^{n}exp(-\text{∑}_{j=2}^t\frac{\left|U_{(j),i}-U_{(j-1),i}\right|}{{\zeta }_U})\times exp(-\frac{1}{2{\sigma }_V^2}{\Vert V_1\Vert }_2^2)\times exp(-\text{∑}_{j=2}^t\frac{\left|V_j-V_{j-1}\right|}{{\zeta }_V})\end{array}$$

Taking the logarithm on both sides and take the missing values into account, then we have:

$$\begin{array}{l}\text{Equation}(5)\iff \\ \{U^{*},V^{*}\}=\underset{U,V}{\mathit{\text{arg}}}min\{\text{∑}_{i=1}^n\text{∑}_{j=1}^t{({\mathit{\text{R}}}_{i,j}-U_{(j),i}^T{V}_j)}^2\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+\alpha \text{∑}_{j=1}^t{\Vert U_{(j),1}\Vert }_2^2+\gamma \text{∑}_{i=1}^n\text{∑}_{j=2}^t{\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1+\beta {\Vert V_1\Vert }_2^2+\mathrm{\lambda }\text{∑}_{j=2}^t{\Vert V_j-V_{j-1}\Vert }_1\}\end{array}$$

(6)

where $\alpha =\frac{{\sigma }_R^2}{{\sigma }_U^2}$, $\gamma =\frac{{\sigma }_R^2}{{\zeta }_U}$, $\beta =\frac{{\sigma }_R^2}{{\sigma }_V^2}$ and $\lambda =\frac{{\sigma }_R^2}{{\zeta }_V}$ are the regularization parameters.

As a supplement, we have the following comments on Model (6):

• On the selection of the Gaussian prior: In our model, since no prior information is available for the columns of the matrices, U₍₁₎ and V, hence, according to the max entropy principle [26], a reasonable choice for them is the Gaussian prior distribution.

• On the ability to formalize the dynamics of the sensor data: The ability to characterize the dynamics of the sensor data lies in the terms $\gamma {{\sum }_{i=1}^n{\sum }_{j=2}^t\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1$ and $\mathrm{\lambda }{{\sum }_{j=2}^t\Vert V_j-V_{j-1}\Vert }_1$. Obviously, for any two adjacent time points, T_j−1 and T_j, if the interval is small enough (namely, |T_j – T_j−1| → 0), then for any area, S_i, the values, R_i,j−1 and R_i,j, should be close to each other (namely, |R_i,j – R_i,j−1| → 0). This can been enforced by tuning the parameters, γ and λ (see the following elaboration).

  First of all, since for any x ∈ ℝⁿ, ‖x‖₂ ≤ ‖x‖₁, we have:

  $$\begin{array}{ll}{\Vert R_{i,j}-R_{i,j-1}\Vert }_2& ={\Vert U_{(j),i}^T{V}_j-U_{(j-1),i}^T{V}_{j-1}\Vert }_2\\ & ={\Vert {(U_{(j),i}-U_{(j-1),i})}^T{V}_j+U_{(j-1),i}^T(V_j-V_{j-1})\Vert }_2\\ & \le {\Vert U_{(j),i}-U_{(j-1),i}\Vert }_2{\Vert V_j\Vert }_2+{\Vert U_{(j-1),i}\Vert }_2{\Vert V_j-V_{j-1}\Vert }_2\\ & \le {\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1{\Vert V_j\Vert }_2+{\Vert U_{(j-1),i}\Vert }_2{\Vert V_j-V_{j-1}\Vert }_1\end{array}$$

  (7)
  Secondly, it is obvious that greater regularization parameters (i.e., α, β, γ and λ) result in smaller corresponding multipliers (i.e., ${{\sum }_{j=1}^t\Vert U_{(j),1}\Vert }_2^2$, ${\Vert V_1\Vert }_2^2$, ${\sum }_{i=1}^n{\sum }_{j=2}^t{\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1$ and ${\sum }_{j=2}^t{\Vert V_j-V_{j-1}\Vert }_1$. In particular:
  $$\gamma \to \infty \Rightarrow {\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1\to 0$$

  (8)
  and:
  $$\lambda \to \infty \Rightarrow {\Vert V_j-V_{j-1}\Vert }_1\to 0$$

  (9)
  Hence, combining Equations (7)-(9), when |T_j – T_j−1| → 0, we can simply take γ → ∞ and λ → ∞ and achieve ‖R_j – R_j−1‖₂ → 0.

  On the other side, when |T_j – T_j−1| → ∞, as is well known, the values in R_j and R_j−1 are regarded as being independent. In this case, we can take γ → 0 and λ → 0, allowing R_j to be irrelevant to R_j−1.

• On the ℓ₁ norm: It is straightforward to verify that, if we replace the ℓ₁ terms in Equation (6) with the ℓ₂ terms (equivalently, use the Gaussian distribution instead of the Laplace distribution in Assumptions (III.) and (IV.)), e.g., replacing ‖V_j – V_j−1‖₁ with ${\Vert V_j-V_{j-1}\Vert }_2^2$, ‖R_j – R_j−1‖₂ can still be bounded via tuning the regularization parameters, γ and λ. The reason for adopting the ℓ₁ norm here is two-fold: Firstly, as shown above, the ℓ₁ terms can lead to the bounded difference norm, ‖R_j – R_j−1‖₂, and hence, the proposed model accommodates the ability to characterize the dynamics of the sensor data; secondly, according to the recent emerging works on compressed sensing [27,28], under some settings, the behavior of the ℓ₁ norm is similar to that of the ℓ₀ norm. In terms of our model, this result indicates that the ℓ₁ terms can restrict not only the magnitudes of the dynamics happening to the features, but also the number of features that changed in adjacent time points. In other words, with ℓ₁ norms, our model gains more expressibility.

5. The Algorithm

Below, we present the algorithm to solve Model (6). We denote:

$$W=\underset{U,V}{\mathit{\text{arg}}}min\text{∑}_{i=1}^n\text{∑}_{j=1}^t{(R_{i,j}-U_{(J),i}^T{V}_j)}^2\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+\alpha \text{∑}_{j=1}^t{\Vert U_{(j),1}\Vert }_2^2+\beta {\Vert V_1\Vert }_2^2+\gamma \text{∑}_{i=1}^n\text{∑}_{j=2}^t{\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1+\lambda \text{∑}_{j=2}^t{\Vert V_j-V_{j-1}\Vert }_1$$

Apparently, W is convex with respect to U_(j),i and V_j (1 ≤ i ≤ n and 1 ≤ j ≤ t). Therefore, we can obtain the local minimum solution via coordinate descent [29].

First, we introduce the signum function, sgn, for a real variable, x:

$$\mathit{\text{sgn}}(x)=\{\begin{array}{ll}1& \text{if}x>0\\ -1& \text{if}x<0\hfill \\ 0& \text{if}x=0\end{array}$$

For X = [x₁, x₂, …x_n]′ ∈ ℝⁿ, we denote sgn(X) = [sgn(x₁), sgn(x₂), …, sgn(x_n)]′.

Then, we calculate the partial subgradient of W with regard to U_(j),i (1 ≤ i ≤ n and 1 ≤ j ≤ t) as follows:

• For j = 2, 3, …, t, define F_j,i,1 = γsgn(U_(j),i – U_(j−1),i).

• For j = 1, 2, …, t – 1, define F_j,i,2 = –γF_j+1,i,1.

• Let F_1,i,1 = F_t,i,2= 0, and we have:

• For i = 1, 2 …,n:

  $$\frac{\partial W}{\partial U_{(1),i}}=2\alpha U_{(1),i}-2(R_{i,1}-U_{(1),i}^T{V}_j)V_j\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+F_{1,i,1}+F_{1,i,2}$$

• For i = 1,2 …,n and j = 2, 3, …, t:

  $$\frac{\partial W}{\partial U_{(j),i}}=-2({\text{R}}_{i,j}-U_{(j),i}^T{V}_j)V_j\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime }{)+\text{F}}_{j,i,1}+F_{j,i,2}$$

Similarly, we calculate the partial subgradient of W with regard to V_j (1 ≤ j ≤ t):

• For 2 ≤ j ≤ t, define G_j,1 = λsgn(V_j – V_j−1).

• For 1 ≤ j ≤ t – 1, denote G_j,2 = – G_j+1,1.

• Let G_1,1 = G_t,2 = 0, and we have:

  $$\frac{\partial W}{\partial V_1}=2\beta V_1-2\text{∑}_{i=1}^n(R_{i,1}-U_{(1),i}^T{V}_1)U_{(1),i}\mathbb{\text{I}}(R_{i,1}{\ne }^{\prime }{\perp }^{\prime })+G_{1,1}+G_{1,2}$$
  and for j > 1:
  $$\frac{\partial W}{\partial V_j}=-2\text{∑}_{i=1}^n(R_{i,j}-U_{(j),i}^T{V}_j)U_{(j),i}\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+G_{j,1}+G_{j,2}$$

Finally, with the results above, we present the solution algorithm in Algorithm 1.

6. Applications on Missing Sensor Data Imputation

In this section, we evaluate our approach through two large-sized datasets and compare the results with two state-of-the-art algorithms in terms of parametric sensitivity, convergence and missing data recovery performance. The following paragraphs describe the set-up, evaluation methodology and the results obtained. To simplify the parameter tuning, we set α = β and λ = γ in the algorithm implementation.

6.1. Evaluation Methodology

Three state-of-the-art algorithms are selected for comparison to the proposed MMP model. The first one is the k-nearest neighbor-based imputation model [1]. As a local model, for every missing entry, R_i,j, the knnmethod takes the estimation, R̂_i,j, as the mean of the k nearest neighbors to it. Let (x) be the set consisting of the k non-empty entries to x; then:

$${\widehat{R}}_{i,j}=\frac{1}{k}\sum _{R_{i,l}\in \mathcal{N}(R_{i,j})}{R}_{i,l}$$

(10)

The second algorithm is the probabilistic principle components analysis model (PPCA) [30,31], which has achieved state-of-the-art performance in the missing traffic flow data imputation problem [31]. Denote the observations of the incomplete matrix, R, as R_o. Let x ∼ N(0, I); to estimate the missing values, PPCA first fits the parameters μ and C with:

$$\left\{\mu *,C*\right\}=\underset{\mu ,C}{\mathit{\text{arg}}}max\mathit{\text{Pr}}(x\mid R_o)~N(x{\sigma }^{-2}{(I+{\sigma }^{-2}C{C}^T)}^{-1}C{R}_o,{(I+{\sigma }^{-2}C{C}^T)}^{-1})$$

(11)

where σ is the tunable parameter. Then, with R_o, μ* and C*, it takes the estimation of the missing values (denoted as R_m) as:

$${\widehat{R}}_m=\underset{R_m}{\mathit{\text{arg}}}maxN(R_m\mid C^{T}x,{\sigma }^{2}I)$$

(12)

The third algorithm is the probabilistic matrix factorization model (PMF) [7], one of the most popular algorithms targeting the Netflix matrix completion problem. PMF first seeks the low rank matrices, U and V, that:

$$\left\{U*,V*\right\}=\underset{U,V}{\mathit{\text{arg}}}min\text{∑}_{i=1}^n\text{∑}_{j=1}^t{(R_{i,j}-U_i^T{V}_j)}^2\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+\alpha \text{∑}_{i=1}^n{\Vert U_i\Vert }_2^2+\beta \text{∑}_{j=1}^t{\Vert V_j\Vert }_2^2$$

(13)

then for the missing entry R_i,j, it takes the estimation as

$${\hat{R}}_{i,j}=U_i^T{V}_j$$

(14)

Algorithm 1: The Multi-Matrices Factorization Algorithm.
	Input: matrix R; number of the latent features, d; learning rates, η1, η2, η3 and η4; regularization parameters, α and λ; threshold ϵ.
	Output: the estimated matrix, R̂.
	// Initialize U and V.
1	Draw random vectors, U(1),1, U(1),2, …, U(1),n, V1 ∼ N(0, I);
2	for j = 2; j ≤ t; j + + do
3	Let Vj = Vj−1 + Z; here Z ∼ Laplace(0, 1);
4	end
5	for j = 2; j ≤ t; j + + do
6	for i = 1; i ≤ n; i + + do
7	Let U(j),i = U(j-1),i + Z; here Z ∼ Laplace(0, 1);
8	end
9	end
	// Coordinate descent.
10	$W_1={\sum }_{i=1}^n{\sum }_{j=1}^t{(R_{i,j}-U_{(j),i}^T{V}_j)}^2\mathbb{\text{I}}(R_{i,j}{\ne }^{\prime }{\perp }^{\prime })+\alpha {\sum }_{j=1}^t{\Vert U_{(j),1}\Vert }_2^2+\beta {\Vert V_1\Vert }_2^2+\gamma {\sum }_{i=1}^n{\sum }_{j=2}^t{\Vert U_{(j),i}-U_{(j-1),i}\Vert }_1+\lambda {\sum }_{j=2}^t{\Vert V_j-V_{j-1}\Vert }_1$;
11
12	W2 = inf;
13	while |W2 – W1| > ϵ do
14	W2 = W1;
15	for i= 1, 2 …, n do
16	Let $U_{(1),i}^{\mathit{\text{new}}}=U_{(1),i}-{\eta }_1\frac{\partial W_1}{\partial U_{(1),i}}$;
17	end
18	for j > 1 and i = 1, 2 …, n do
19	Let $U_{(j),i}^{\mathit{\text{new}}}=U_{(j),i}-{\eta }_2\frac{\partial W_1}{\partial U_{(j),i}}$
20	end
21	$V_1^{\mathit{\text{new}}}=V_1-{\eta }_3\frac{\partial W_1}{\partial V_1}$;
22	for j = 2 …, t do
23	Let $V_j^{\mathit{\text{new}}}=V_j-{\eta }_4\frac{\partial W_1}{\partial V_j}$;
24	end
25	Replace all U(j),is with $V_j^{\mathit{\text{new}}}s$ and Vjs with , recompute W1;
26	end
27	return R̂, where ${\widehat{R}}_j=U_{(j)}^T{V}_j$;

These three algorithms, as well as the proposed MMF, are employed to perform missing imputations for the incomplete matrix, R, on the same datasets.

The testing protocol adopted here is the Given X (0 < X < 1) protocol [32], i.e., given a matrix, R, only X percent of its observed entries are revealed, while the remaining observations are concealed to evaluate the trained model. For example, a setting with X = 10% means that the algorithm is trained with 10% of the non-missing entries, and the rest of the 90% non-missing ones are held and are to be recovered. In both of the experiments on synthetic and real datasets, the data partition is repeated five times, and the average results, as well as the standard deviations over the five repetitions are recorded.

Similar to many other missing imputation problems [1,3-5,7,13,33-35], we employ the root mean square error (RMSE) to depict the distance between the real values and the estimations: Let S = {s₁, s₂, … s_n} be the test dataset and Ŝ = {ŝ₁, ŝ₂, … ŝ_n} be the estimated set; here, ŝ_i is the estimation of s_i. Then, the RMSE of the estimation is given by $\sqrt{\frac{1}{n}{\sum }_{k=1}^n{(s_k-{\widehat{s}}_k)}^2}$.

6.2. Synthetic Validation

To conduct a synthetic validation of the studied approaches, we randomly draw a 100 × 10, 000 matrix R using the procedure detailed in Algorithm 2. The rows in R correspond to the areas, S₁, S₂, …, S_n, and the columns correspond to the time. Thus, R_i,j represents the data collected in S_i at time T_j. Notably, the parameter, r_i,j, in Algorithm 2 is used to control the magnitude of the variation happening to S_i from T_j−1 to T_j. Combining lines 4 and 5, we have: for i ∈ [1, n] and j ∈ [1, t], $|\frac{R_{i,j}-R_{i,j-1}}{R_{i,j-1}}|=|r_{i,j}|\le 0.1$. This constraint ensures that the data collected in S_i does not change too much over time T_j−1 to T_j.

Algorithm 2: Synthetic Data Generating Procedure.
1	for i = 1, 2, …, n do
2	Draw Ri,1 ∼ N(0, 1);
3	for j = 2, 3, …, t do
4	Let ri,j ∼ Uniform(−0.1, 0.1);
5	Ri,j ← Ri,j−1 + ri,jRi,j−1;
6	end
7	end

We first evaluate the sensitivity of the proposed algorithm to the regularization parameters, α and λ. Half of the entries in R are randomly selected as testing data and recovered using the remaining 50% as the training data. Namely, we take X = 50% in the Given X protocol. In the experiments, we first fix α = 0.01, tune λ via λ = 0.01 × 2ⁿ (n = 0, 1, … 7) and, then, do the reverse by changing α via α = 0.01 × 2ⁿ (n = 0, 1, … 7), but setting λ = 0.01. The average RMSEs with the same parameter settings on different data partitions are summarized in Figure 2.

In Figure 2, the RMSE-1 curve represents the recovery errors obtained by fixing α and changing the value of λ. The RMSE-2 curve corresponds to the errors with different α values and fixed λ. We can see that even when λ is expanded by more than 100 times (2⁷ = 128), the RMSE still remains stable. A similar result also appears in the experiments on the parameter, α, where a significant change of the RMSE only occurs when n is greater than six, i.e., when α is expanded more than 60 times.

The second experiment we conduct is to study the prediction ability of the proposed algorithm, as well as that of the comparison algorithms. In the Given X protocol, we set X = 10%, 20%, …, 90% in sequence. Then, for each X value, we perform missing imputations via our algorithm and the comparison algorithms. In the all implementations, we set k = 5 for the knn model. As with the MF-based algorithms, we examine their performance with respect to the latent feature dimension d = 10 and d = 30 respectively. Furthermore, in the implementation of MMF, we fix α = λ = 0.1, η₁ = η₂ = η₃ = η₄ = 0.01. All results are summarized in

Table 1. Recovery errors on the synthetic dataset (mean ± std).

		10%	20%	30%	40%	50%	60%	70%	80%	90%
	knn	13.41 ± 0.62	6.80 ± 0.19	4.44 ± 0.06	3.12 ± 0.05	2.27 ± 0.02	1.81 ± 0.09	1.72 ± 0.01	1.69 ± 0.05	1.62 ± 0.06
d = 10	PPCA	17.09 ± 2.08	20.24 ± 1.30	22.75 ± 1.32	23.96 ± 2.19	20.79 ± 1.00	16.67 ± 1.34	18.68 ± 0.81	11.98 ± 0.70	5.11 ± 3.10
	PMF	3.23 ± 0.23	3.33 ± 0.19	3.29 ± 0.12	3.34 ± 0.07	3.29 ± 0.09	1.69 ± 0.04	1.83 ± 0.02	1.81 ± 0.03	1.85 ± 0.04
	MMF	3.07 ± 0.07	2.21 ± 0.10	2.14 ± 0.09	1.98 ± 0.06	1.93 ± 0.06	1.92 ± 0.04	1.75 ± 0.03	1.84 ± 0.02	1.80 ± 0.03
d = 30	PPCA	19.65±3.01	22.48 ± 0.86	24.86 ± 1.54	23.99 ± 0.61	22.67 ± 0.49	20.22 ± 1.46	17.53 ± 0.70	14.23 ± 1.70	11.14 ± 1.72
	PMF	3.20 ± 0.21	3.32 ± 0.13	3.35 ± 0.09	3.36 ± 0.11	3.31 ± 0.07	1.78 ± 0.04	1.81 ± 0.02	1.79 ± 0.02	1.86 ± 0.11
	MMF	3.06 ± 0.05	2.17 ± 0.08	2.07 ± 0.05	1.94 ± 0.03	1.74 ± 0.03	1.62 ± 0.02	1.64 ± 0.03	1.65 ± 0.01	1.69 ± 0.03

.

As shown in Table 1, when X is large, e.g., X ≥ 50%, knn is competitive with the matrix factorization methods, while in the other situations, the MF methods outperform it significantly. In terms of the MF-based methods, we find that our algorithm outperforms PPCA significantly in all settings. The RMSEs of our algorithm are at most roughly 20% of that of PPCA. Specifically, for X ∈ [30%; 80%], the RMSEs of the proposed algorithm are even only 10% of that of PPCA. We can also observe that the parameter, d, has a different impact on the performances of our algorithm and PPCA: When d changes from 10 to 30, most of the RMSEs of PPCA increase evidently, while for our algorithm, the RMSEs are reduced by roughly 5%.

When compared with PMF, our algorithm also performs better in most of the settings: PMF achieves lower RMSEs than MMF only in two cases, in which d = 10, X = 60% and d = 10; X = 80% respectively. Another interesting finding is that the promotion of the feature number, d (from 10 to 30), has little impact on the performance of PMF.

We also exam the convergence speed of the proposed algorithm. In the missing recovery experiments conducted above, for each X setting, we record the average RMSE of the recovered results after every 10 iterations of all data partitions. We can see from Figure 3 that, for all X values, the errors drop dramatically in the first 20 iterations and remain stable after the first 100 iterations. We can conclude that the proposed algorithm converges to the local optimization solutions after around 100 iterations.

6.3. Application to Impute the Missing Traffic Speed Values

To evaluate the feasibility of the proposed approach on real-world applications, in this section, we conduct another experiment on a traffic speed dataset, which was collected in the urban road network of Zhuhai City [36], China, from April 1, 2011 to April 30, 2011. The data matrix, R, consists of 1,853 rows and 8,729 columns. Each row corresponds to a road, and each column corresponds to a five minute-length time interval. All columns are arranged in ascending order of time. An entry, R_i,j (1 ≤ i ≤ 1, 853,1 ≤ j ≤ 8729), in R is the aggregate mean traffic speed of the ith road in the jth interval. Since all the data in R are collected by floating cars [37], the value of R_i,j could be missing if there is no car on the i th road in the j th time interval. Our statistics show that in R, there are nearly half of the entries, i.e., eight million entries are missing values.

We perform missing imputation on matrix R using the studied algorithms with parameter settings k = 5 and d = 10. In the implementation of MMF, we fix α = 0.25, λ = 0.5, η₁ = η₂ = η₃ = η₄ = 0.5. We summarize all results in

Table 2. Recovery errors on the transportation dataset (mean ± std).

	10%	20%	30%	40%	50%	60%	70%	80%	90%
knn	40.47 ± 0.02	31.79 ± 0.01	25.41 ± 0.02	21.35 ± 0.00	18.33 ± 0.00	15.89 ± 0.01	13.73 ± 0.01	11.67 ± 0.00	9.45 ± 0.02
PPCA	17.90 ± 0.01	17.36 ± 0.01	13.00 ± 0.02	12.25 ± 0.01	11.47 ± 0.01	11.31 ± 0.03	11.19 ± 0.02	11.14 ± 0.04	11.16 ± 0.10
PMF	14.41 ± 0.01	12.83 ± 0.03	12.43 ± 0.01	12.33 ± 0.02	12.36 ± 0.01	12.35 ± 0.01	12.13 ± 0.00	11.99 ± 0.03	11.96 ± 0.02
MMF	11.79 ± 0.02	11.51 ± 0.01	11.43 ± 0.01	11.05 ± 0.02	11.05 ± 0.01	11.01 ± 0.00	10.83 ± 0.01	10.69 ± 0.01	10.70 ± 0.02

, from which both the feasibility and effectiveness of MMF are well verified. In detail, when X is large enough, e.g., X ≥ 80%, knn is competitive, while in the other cases, knn cannot work as well as the MF based algorithms. As for the MF algorithms, we see that the proposed MMF outperforms PPCA and PMF in all X settings. Particularly, when the observations are few (X = 10% and X = 20%), the errors of our algorithm reduce by 33% compared to those of PPCA and by 10% compared to those of PMF, respectively. When X > 20%, the RMSE differences between PPCA and our algorithm tend to be slight, but the overall errors of PPCA are roughly 3% ∼ 5% higher than those of MMF. For PMF, the RMSEs remain about 10% higher than MMF in all settings.

7. Conclusion

Missing estimation is one of the main concerns in current studies on sensor data-based applications. In this work, we formulate the estimation problem as a matrix completion one and present a multi-matrices factorization model to address it. In our model, each column, R_j, of the target matrix, R, is approximated by the product of a spatial feature matrix, U_(j), and a temporal feature vector, V_j. Both U_j and V_j are time dependent, and hence, their product accommodates the ability to describe the time variant sensor data. We also present a solution algorithm to the factorization model. Empirical studies on a synthetic dataset and real sensor data show that our approach outperforms the comparison algorithms.

Reviewing the present work, it is notable that the proposed model only incorporates the temporal structure information, while the information on the spatial structure is disregarded, e.g., the data collected in two adjacent areas, S_k and S_l, should be close to each other. Hence, our next step is to extend our model with more complex structured data.