# Matrix Factorization for Evolution Data

Authors: Xiao-Yu Huang; Xian-Hong Xiang; Wubin Li; Kang Chen; Wen-Xue Cai; Lei Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We study a matrix factorization problem, that is, to find two factor matricesU and V such that R≈UT×V, where R is a matrix composed of the values of the objects O1,O2,…,On at consecutive time points T1,T2,…,Tt. We first present MAFED, a constrained optimization model for this problem, which straightforwardly performs factorization on R. Then based on the interplay of the data in U, V, and R, a probabilistic graphical model using the same optimization objects is constructed, in which structural dependencies of the data in these matrices are revealed. Finally, we present a fitting algorithm to solve the proposed MAFED model, which produces the desired factorization. Empirical studies on real-world datasets demonstrate that our approach outperforms the state-of-the-art comparison algorithms.


---

## Body

## 1. Introduction

The essence of thematrix factorization (MF) problem is to find two factor matrices U and V, such that their product can approximate a given matrix R; that is, R≈UTV. As a fundamental model in machine learning and data mining, MF methods have been widely used in various applications, such as collaborative filtering [1, 2], social network analysis [3], text analysis [4], image analysis [5, 6], and biology analysis [7].In this work, we study a special variety of the MF problem. Let the matrixR∈Rn×t consist of the values of the objects O1,O2,…,On at a serial of time points T1,T2,…,Tt, where the entry Ri,j is the value of Oi at Tj. Our goal is to find the suitable factor matrices U and V, which can not only approximate R, but also depict theevolution of the objects over time.A typical application of this work is to estimate the missing historical traffic speed data, which is necessary for many transportation information systems [8–10]. Given an urban road network with n roads L1,L2,…,Ln, we let Oi(1≤i≤n) correspond to the traffic speed of Li and let entry Ri,j be the speed of Li at Tj; then R is composed of speed values of the roads L1~Ln at the time T1~Tt. Especially, if the value of Ri,j was not collected, we say it ismissing and denote as Ri,j= “⊥”. With this representation, to estimate the missing values, we can first fit U and V with the nonmissing entries of R, and then for each Ri,j= “⊥”, take its estimation as R^i,j=UiTVj, where Ui/Vj is the ith/jth column of U/V.A main difference between this work and the other existing MF models is that we take the time evolution effects into account. Well-studied MF models, such as nonnegative matrix factorization [5], max margin matrix factorization [11, 12], and probabilistic matrix factorization [13], are based on the i.i.d assumption, which implicates that the behavior of the objects evolved with time is ignored and treated as being independent of time, and thus their abilities on describing time-varying data are poor. To tackle this issue, the factors that affect the evolution of the objects are explicitly modeled in our contribution, and that enables our model to interplay with the time dependent information. Furthermore, our empirical studies confirm that the proposed model can scale well with size of the data.The remainder of this paper is organized as follows. Section2 briefly summarizes the notations used in the paper. Section 3 reviews the works on matrix factorization. Section 4 presents our matrix factorization model as well as the statistical mechanism of the model. The proposed algorithm is introduced in Section 5. Section 6 is devoted to the analysis of the experiments. Our conclusions and future works are presented in Section 7.


## 2. Notations

For a vectorV=[v1,v2,…,vn]′∈Rn, we denote its 2-Norm as ∥V∥2, where
(1)∥V∥2=∑i=1nvi2‍.
For a matrix X∈Rn×m, we let Xk be the kth column of X and denote itsFrobenius norm as
(2)∥X∥F=∑i=1n∑j=1mXi,j2.


## 3. Matrix Factorization

The essence of matrix factorization is to find the suitable factor matricesU and V, such that their product can approximate a given matrix R. In principle, the MF problem can be formulated as the optimization model below:
(M1){U*,V*}=min⁡U,V⁡Loss(UTV,R),
where the function Loss is used to measure thecloseness of the approximation UTV to the target R. In general, Loss(UTV,R) can be decomposed into the sum of the pairwise loss between the entries of UTV and R; that is, Loss(UTV,R)=∑i=1n‍∑j=1m‍Loss((UTV)i,j,Ri,j). Most common used forms of the Loss function include the square loss (Loss(x,y)=(x-y)2) [1, 13, 14], the 0-1 Loss (Loss(x,y)=I(x=y)) [11], and the divergence loss (Loss(x,y)=xlog⁡(x/y)-x+y) [6].Notably, if{U*,V*} is a solution of (M1), then for any scalar ζ>0, {ζU*,(1/ζ)V*} is another solution, and hence the problem (M1) is ill posed. To overcome this obstacle, various constraints on U and V are introduced, such as constraints on the entries [5], constraints on the sparseness [15, 16], constraints on the norms [13, 17], and constraints on the ranks [18, 19]. All these constraints, from the perspective of the statistical learning theory, can be regarded as thelength of the model to be fitted. According to the minimum description length principle [20, 21], smaller length means better model, and thus most of them can be incorporated into (M1) as additional regularized terms; that is,
(M2){U*,V*}=min⁡U,V⁡Loss(UTV,R)+P(U,V),
where the regularization factor P(U,V) corresponds to the constraints on U and V.As a transductive model,(M2) has many appealing mathematical properties, such as the generalization error bound [22] and the exactness [14, 23]. However, as well known, when compared with the generative model, a main weakness of the transductive model is that it can hardly be used to describe the relation in the data. In particular, for our studied problem, even though the model (M2) may work well, it is difficult to express the interplay between the model and the evolutions of the objects.


## 4. The Proposed Model

In this section the proposed model for matrix factorization for evolution data (MAFED) is presented. We first formalize MAFED as a constraint optimization model and then present a probabilistic graphical model, the solution of which is equivalent to the MAFED model. Finally we attain the generative interpretation of MAFED, with which we elaborate the ability of MAFED to describe the data relation. The last subsection is devoted to the solution algorithm for MAFED.


### 4.1. The Model

LetR∈Rn×t be the matrix generated by n independent objects O1,O2,O3,…,On at t consecutive time points T1,T2,…,Tt; then our MAFED model can be formulated as follows: (3a){U*,V*}=argU,V⁡min⁡∑i=1n∑j=1t(Ri,j-UiTVj)2‍‍(3b)s.t.∥Ui∥2⩽A,i=1,2,…,n(3c)s.t.∥V1∥2⩽B(3d)s.t.∥Vj+1-Vj∥2⩽Cj,j=1,2,…,t-1.Here Ui and Vj correspond to the ith column of U and the jth column of V, respectively.From constraints (3b)–(3d), we can see that the roles of U and V in (3a) are asymptotic. We call U the object matrix, of which the ith column (Ui) is the invariant feature vector of Oi. Similarly, we call V the environment matrix, of which the jth column (Vj) is the time varying environment feature vector at Tj. Note that for each Vj, the range it takes effects in is global; that is, for each time point Tj, all the objects O1,O2,…,On share the same environment feature vector Vj. As a result, every entry Ri,j(1≤i≤n,1≤j≤t) is composed by the object’s intrinsic feature vector Ui and the environment feature vector Vj, and hence for each object Oi, its evolutions over the time T1,T2,…,Tt are controlled by T1,T2,…,Tt. To be more illustrative, let us consider such an example. Let R be the speed matrix of the roads L1,L2,…,Ln at the time points T1,T2,…,Tt, let U be the feature matrix of the roads, and let V be the feature matrix of the environment; then Ui corresponds to the feature description of Li and the entries of Ui may consist of the intrinsic features of the road, such as the surface rough, the number of lanes, and the role in the traffic network; similarly, for V, every Vj corresponds to the description of the environment at Tj, and the entries of Vj can be the time, the weather, the visibility, and so on. Therefore, for the product Ri,j=UiTVj, it can be interpreted as the speed of Li at Tj is composed by the road feature Ui and the environment feature Vj.The ability of MAFED to describe the time evolution effect lies in the constraint (3d). As well known, a fundamental characteristic of the evolution is for the all objects; when the interval between two adjacent time points, for example, Tj and Tj+1, is small enough, the corresponding values of the objects, that is, Rj and Rj+1, should tend to remain the same. To satisfy this constraint, in (3d) we introduce the tunable values C1,C2,…,Ct-1. When |tj-tj+1|→0, we let Cj→0. Now since
(4)∥Rj-Rj+1∥2=∥UTVj-UTVj+1∥2≤∥U∥2∥Vj-Vj+1∥2,∥U∥2≤∥U∥F≤∑i=1n∥Ui∥2‍≤n×A,
we have
(5)∥Rj-Rj+1∥2≤n×A×∥Vj-Vj+1∥2≤n×A×Cj.
When n and A are fixed, Cj→0 leads to n×A×Cj→0, and so via (3d) it is guaranteed that |tj-tj+1|→0⇒||Rj-Rj+1||2→0.On the other hand, when the time interval|tj-tj+1|→+∞, the values that Rj and Rj+1 take will tend to be independently. For this case, we let Cj→+∞, indicating that Rj+1 can take its value regardless of Rj.With the constraints (3c) and (3d), we can control the value of ∥Vj∥2(j=2,…,t) via
(6)∥Vj∥22=∥V1+V2-V1+⋯+Vj-Vj-1∥22≤∥V1∥22+∥V2-V1∥22+⋯+∥Vj-Vj-1∥22≤B2+∑k=1j-1Ck2.The result above shows that we can bound the sum∑j=1t‍∥Vj∥22 by selecting appropriate parameters B and C. Besides, from the constraint (3b), the sum ∑i=1n‍∥Ui∥22 is bounded by nA2. According to [12, 24], under some suitable situations, the sum 1/2(∑j=1t‍∥Vj∥22+∑i=1n‍∥Ui∥22) is the convex envelop of rank⁡ (UTV), and hence MAFED accommodates the ability to control the rank of the model. On the other side, as shown by Candès and Recht [14] and Candès and Tao [23], when the target matrix islow rank, it is possible to exactly recover it from only a few of its observations. This declaration shows that MAFED is possible to achieve high accuracy in the applications of missing imputations.


### 4.2. The Generative Interpretation

In this section, to investigate how MAFED interplays with the evolution of the data from the view of generative modeling, we first present a probabilistic graphical model (PGM [25]) and then show that the maximum likelihood solution to the proposed PGM is exactly the same as the solution of (3a).First of all, the Lagrange dual of (3a) is
(7){U*,V*}=argU,V⁡min⁡∑i=1n∑j=1t(Ri,j-UiTVj)2+α∑i=1n‍∥Ui∥22+β∥V1∥22+∑j=1t-1‍λj∥Vj+1-Vj∥22s.t.α,β,λ1,λ2,…,λt-1≥0.Hereα, β, and λi(i=1,2,…,t-1) are tunable parameters, corresponding to the upper bounds A, B, and Ci(i=1,2,…,t-1) of constraints (3b)–(3d), respectively, where greater bounds correspond to smaller parameters, and vise versa.Considering the PGM presented in Figure1, we have the following assumptions for U, V, and R.(1)
The columns ofU are from the same Gaussian distribution with mean 0 and covariance matrix σU2I; that is, for 1≤i≤n, we have (here we use d to denote the dimension of Ui)
(8)Pr(Ui∣σU)=(2πσU2)-d/2exp⁡⁡{-∥Ui∥222σU2}.(2)
The columns ofV are linearly dependent in the order of their subscripts with respect to prespecified priors σP and σQ1,σQ2,…,σQT-1; that is,
(9)Pr(V∣σP,σΛ)=Pr(V1∣σP,σΛ)×∏j=2TPr⁡(Vj∣Vj-1,σP,σΛ)‍.Here for clarity, we letΛ={Q1,Q2,…,Qt-1} and σΛ={σQ1,σQ2,…,σQT-1}. We also assume V1 is a Gaussian random vector with mean 0 and covariance σP2I and Vj (j>1) is a Gaussian random vector with mean Vj-1 and covariance σQj-12I; that is,
(10)Pr(V1∣σP,σΛ)=(2πσP2)-d/2exp⁡⁡{-∥V1∥222σP2},Pr(Vj∣Vj-1,σP,σΛ)=(2πσQj-1)-d/2exp⁡⁡{-∥Vj-Vj-1∥222σQj-12}.(3)
The(i,j)th entry of R (1≤i≤N, 1≤j≤T) is a Gaussian random variable with mean UiTVj and variance σR2; that is,
(11)Pr(Ri,j∣UiTVj,σR2)=(2πσR2)-1/2exp⁡⁡{-(Ri,j-UiTVj)22σR2}.With these assumptions, we in the following part of this section show that, for a given matrix R and priors σU, σP, σΛ,  σR, model (7) is equivalent to the following maximum likelihood fitting problem:
(12){U*,V*}=argU,V⁡max⁡Pr(U,V∣R,σU,σP,σΛ,σR).
In fact, according to the Bayesian theorem, we have
(13)Pr(U,V∣R,σU,σP,σΛ,σR)=Pr(U,V,R∣σU,σP,σΛ,σR)Pr(R∣σU,σP,σΛ,σR).
Since R is observed and σU, σP, σΛ, and σR are prespecified, the denominator Pr(R|σU,σP,σΛ,σR) can be treated as a constant. Therefore, we get
(14)(12)⟺{U*,V*}=argU,V⁡max⁡Pr(U,V,R∣σU,σP,σΛ,σR).
Combining Figure 1 and the assumptions (1)~(3), we have
(15)Pr(R,U,V∣σU,σP,σΛ,σR)=Pr⁡(R∣U,V,σU,σP,σΛ,σR)×Pr(U,V∣σU,σP,σΛ,σR)=Pr(R∣U,V,σR)×Pr(U∣σU)×Pr(V∣σP,σΛ)=∏i=1n∏j=1tPr(Ri,j∣Ui,Vj,σR)⁡‍‍×∏i=1nPr⁡(Ui∣σU)×Pr⁡(V1∣σP)‍×∏j=2tPr⁡(Vj∣Vj-1,σΛ)‍∝exp⁡(-12σR2∑i=1n∑j=1t(UiTVj-Ri,j)2)×exp⁡⁡(-12σU2∑i=1t‍∥Ui∥22)×exp⁡⁡(-12σP2∥V1∥22)×exp⁡(-∑j=1t-1‍12σQj2∥Vj+1-Vj∥22).
Taking the logarithm on the two sides of the equation, we get
(16)(14)⟺{U*,V*}=argU,V⁡min⁡1σR2×∑i=1n∑j=1t(Ri,j-UiTVj)2‍+1σU2∑i=1N∥Ui∥22‍+1σP2∥V1∥22+∑j=1t-11σQj2∥Vj+1-Vj∥22‍.
Comparing (7) and (16), the PGM model is equivalent to the optimization model (3a) if we let α=σR2/σU2, β=σR2/σP2, and λj=σR2/σQj2(j=1,2,…,t-1).Figure 1
The generative model for MAFED.


## 4.1. The Model

LetR∈Rn×t be the matrix generated by n independent objects O1,O2,O3,…,On at t consecutive time points T1,T2,…,Tt; then our MAFED model can be formulated as follows: (3a){U*,V*}=argU,V⁡min⁡∑i=1n∑j=1t(Ri,j-UiTVj)2‍‍(3b)s.t.∥Ui∥2⩽A,i=1,2,…,n(3c)s.t.∥V1∥2⩽B(3d)s.t.∥Vj+1-Vj∥2⩽Cj,j=1,2,…,t-1.Here Ui and Vj correspond to the ith column of U and the jth column of V, respectively.From constraints (3b)–(3d), we can see that the roles of U and V in (3a) are asymptotic. We call U the object matrix, of which the ith column (Ui) is the invariant feature vector of Oi. Similarly, we call V the environment matrix, of which the jth column (Vj) is the time varying environment feature vector at Tj. Note that for each Vj, the range it takes effects in is global; that is, for each time point Tj, all the objects O1,O2,…,On share the same environment feature vector Vj. As a result, every entry Ri,j(1≤i≤n,1≤j≤t) is composed by the object’s intrinsic feature vector Ui and the environment feature vector Vj, and hence for each object Oi, its evolutions over the time T1,T2,…,Tt are controlled by T1,T2,…,Tt. To be more illustrative, let us consider such an example. Let R be the speed matrix of the roads L1,L2,…,Ln at the time points T1,T2,…,Tt, let U be the feature matrix of the roads, and let V be the feature matrix of the environment; then Ui corresponds to the feature description of Li and the entries of Ui may consist of the intrinsic features of the road, such as the surface rough, the number of lanes, and the role in the traffic network; similarly, for V, every Vj corresponds to the description of the environment at Tj, and the entries of Vj can be the time, the weather, the visibility, and so on. Therefore, for the product Ri,j=UiTVj, it can be interpreted as the speed of Li at Tj is composed by the road feature Ui and the environment feature Vj.The ability of MAFED to describe the time evolution effect lies in the constraint (3d). As well known, a fundamental characteristic of the evolution is for the all objects; when the interval between two adjacent time points, for example, Tj and Tj+1, is small enough, the corresponding values of the objects, that is, Rj and Rj+1, should tend to remain the same. To satisfy this constraint, in (3d) we introduce the tunable values C1,C2,…,Ct-1. When |tj-tj+1|→0, we let Cj→0. Now since
(4)∥Rj-Rj+1∥2=∥UTVj-UTVj+1∥2≤∥U∥2∥Vj-Vj+1∥2,∥U∥2≤∥U∥F≤∑i=1n∥Ui∥2‍≤n×A,
we have
(5)∥Rj-Rj+1∥2≤n×A×∥Vj-Vj+1∥2≤n×A×Cj.
When n and A are fixed, Cj→0 leads to n×A×Cj→0, and so via (3d) it is guaranteed that |tj-tj+1|→0⇒||Rj-Rj+1||2→0.On the other hand, when the time interval|tj-tj+1|→+∞, the values that Rj and Rj+1 take will tend to be independently. For this case, we let Cj→+∞, indicating that Rj+1 can take its value regardless of Rj.With the constraints (3c) and (3d), we can control the value of ∥Vj∥2(j=2,…,t) via
(6)∥Vj∥22=∥V1+V2-V1+⋯+Vj-Vj-1∥22≤∥V1∥22+∥V2-V1∥22+⋯+∥Vj-Vj-1∥22≤B2+∑k=1j-1Ck2.The result above shows that we can bound the sum∑j=1t‍∥Vj∥22 by selecting appropriate parameters B and C. Besides, from the constraint (3b), the sum ∑i=1n‍∥Ui∥22 is bounded by nA2. According to [12, 24], under some suitable situations, the sum 1/2(∑j=1t‍∥Vj∥22+∑i=1n‍∥Ui∥22) is the convex envelop of rank⁡ (UTV), and hence MAFED accommodates the ability to control the rank of the model. On the other side, as shown by Candès and Recht [14] and Candès and Tao [23], when the target matrix islow rank, it is possible to exactly recover it from only a few of its observations. This declaration shows that MAFED is possible to achieve high accuracy in the applications of missing imputations.


## 4.2. The Generative Interpretation

In this section, to investigate how MAFED interplays with the evolution of the data from the view of generative modeling, we first present a probabilistic graphical model (PGM [25]) and then show that the maximum likelihood solution to the proposed PGM is exactly the same as the solution of (3a).First of all, the Lagrange dual of (3a) is
(7){U*,V*}=argU,V⁡min⁡∑i=1n∑j=1t(Ri,j-UiTVj)2+α∑i=1n‍∥Ui∥22+β∥V1∥22+∑j=1t-1‍λj∥Vj+1-Vj∥22s.t.α,β,λ1,λ2,…,λt-1≥0.Hereα, β, and λi(i=1,2,…,t-1) are tunable parameters, corresponding to the upper bounds A, B, and Ci(i=1,2,…,t-1) of constraints (3b)–(3d), respectively, where greater bounds correspond to smaller parameters, and vise versa.Considering the PGM presented in Figure1, we have the following assumptions for U, V, and R.(1)
The columns ofU are from the same Gaussian distribution with mean 0 and covariance matrix σU2I; that is, for 1≤i≤n, we have (here we use d to denote the dimension of Ui)
(8)Pr(Ui∣σU)=(2πσU2)-d/2exp⁡⁡{-∥Ui∥222σU2}.(2)
The columns ofV are linearly dependent in the order of their subscripts with respect to prespecified priors σP and σQ1,σQ2,…,σQT-1; that is,
(9)Pr(V∣σP,σΛ)=Pr(V1∣σP,σΛ)×∏j=2TPr⁡(Vj∣Vj-1,σP,σΛ)‍.Here for clarity, we letΛ={Q1,Q2,…,Qt-1} and σΛ={σQ1,σQ2,…,σQT-1}. We also assume V1 is a Gaussian random vector with mean 0 and covariance σP2I and Vj (j>1) is a Gaussian random vector with mean Vj-1 and covariance σQj-12I; that is,
(10)Pr(V1∣σP,σΛ)=(2πσP2)-d/2exp⁡⁡{-∥V1∥222σP2},Pr(Vj∣Vj-1,σP,σΛ)=(2πσQj-1)-d/2exp⁡⁡{-∥Vj-Vj-1∥222σQj-12}.(3)
The(i,j)th entry of R (1≤i≤N, 1≤j≤T) is a Gaussian random variable with mean UiTVj and variance σR2; that is,
(11)Pr(Ri,j∣UiTVj,σR2)=(2πσR2)-1/2exp⁡⁡{-(Ri,j-UiTVj)22σR2}.With these assumptions, we in the following part of this section show that, for a given matrix R and priors σU, σP, σΛ,  σR, model (7) is equivalent to the following maximum likelihood fitting problem:
(12){U*,V*}=argU,V⁡max⁡Pr(U,V∣R,σU,σP,σΛ,σR).
In fact, according to the Bayesian theorem, we have
(13)Pr(U,V∣R,σU,σP,σΛ,σR)=Pr(U,V,R∣σU,σP,σΛ,σR)Pr(R∣σU,σP,σΛ,σR).
Since R is observed and σU, σP, σΛ, and σR are prespecified, the denominator Pr(R|σU,σP,σΛ,σR) can be treated as a constant. Therefore, we get
(14)(12)⟺{U*,V*}=argU,V⁡max⁡Pr(U,V,R∣σU,σP,σΛ,σR).
Combining Figure 1 and the assumptions (1)~(3), we have
(15)Pr(R,U,V∣σU,σP,σΛ,σR)=Pr⁡(R∣U,V,σU,σP,σΛ,σR)×Pr(U,V∣σU,σP,σΛ,σR)=Pr(R∣U,V,σR)×Pr(U∣σU)×Pr(V∣σP,σΛ)=∏i=1n∏j=1tPr(Ri,j∣Ui,Vj,σR)⁡‍‍×∏i=1nPr⁡(Ui∣σU)×Pr⁡(V1∣σP)‍×∏j=2tPr⁡(Vj∣Vj-1,σΛ)‍∝exp⁡(-12σR2∑i=1n∑j=1t(UiTVj-Ri,j)2)×exp⁡⁡(-12σU2∑i=1t‍∥Ui∥22)×exp⁡⁡(-12σP2∥V1∥22)×exp⁡(-∑j=1t-1‍12σQj2∥Vj+1-Vj∥22).
Taking the logarithm on the two sides of the equation, we get
(16)(14)⟺{U*,V*}=argU,V⁡min⁡1σR2×∑i=1n∑j=1t(Ri,j-UiTVj)2‍+1σU2∑i=1N∥Ui∥22‍+1σP2∥V1∥22+∑j=1t-11σQj2∥Vj+1-Vj∥22‍.
Comparing (7) and (16), the PGM model is equivalent to the optimization model (3a) if we let α=σR2/σU2, β=σR2/σP2, and λj=σR2/σQj2(j=1,2,…,t-1).Figure 1
The generative model for MAFED.


## 5. The Algorithm

We in this section present a fitting algorithm to solve (7). For clarity, we only discuss the case that all the time intervals are of equal length; that is, T2-T1=T3-T2=⋯=Tt-Tt-1. Under this hypothesis, we take λ1=λ2=⋯=λt-1 in model (7). Cases with unequal length time intervals can be handled by removing the constraint λi=λj if Ti+1-Ti=Tj+1-Tj.Let(17)S=argU,V⁡min⁡∑i=1n∑j=1t(Ri,j-UiTVj)2‍‍+α∑i=1n∥Ui∥22+β∥V1∥22+λ∑j=1t-1∥Vj+1-Vj∥22.
It is straightforward to verify that S is convex with respect to Ui(1≤i≤n) and Vj(1≤j≤t), respectively. Therefore we can achieve the local minimum solution of S via coordinate descent [26]. We first calculate the partial derivative with respect to U1,U2,…,Un and V1,V2,…,Vt, respectively. Then we have the following.Fori=1,2,…,n,
(18)∂S∂Ui=2αUi-2∑j=1t(Rij-UiTVj)‍Vj
and for j=2,3,…,t-1,
(19)∂S∂Vj=2λ(2Vj-Vj-1-Vj+1)-2∑i=1n(Ri,j-UiTVj)Ui.
Similarly,
(20)∂S∂V1=2βV1-2∑i=1n(Ri,1-UiTV1)‍Ui-2λ(V2-V1),∂S∂Vt=2λ(Vt-Vt-1)-2∑i=1n(Ri,t-UiTVt)‍Ui.
With the equations above, Algorithm 1 is presented.Algorithm 1:MAFED(R,d,η1,η2,η3,α,β,λ,ϵ).
Input: matrix R; number of the latent features d; learning rates η1, η2 and η3;
regularization parameters α, β and λ; threshold ϵ
Output: the estimated matrix UTV.
//Initialize  U  and  V.
(1) Generate random vectorsU1,U2,…,Un,V1~N(0,I);
(2)forj=2; j≤t; j++do
(3)Generate Z with N(0,I)
(4)Let Vj=Vj-1+Z
(5)end
//Coordinate descent.
(6)S1=∑i=1,j=1i=n,j=t‍(Ri,j-UitVj)2+α∑i=1n‍∥Ui∥22+β∥V1∥22+λ∑j=1t-1‍∥Vj+1-Vj∥22;
(7)S2=inf⁡;
(8)while|S2-S1|>ϵdo
(9)S2=S1;
(10)fori=1,2,…,ndo
(11) Let Uinew=Ui-η1(∂S/∂Ui);
(12)end
(13)V1new=V1-η2∂S∂V1;
(14)forj=2,…,tdo
(15) Let Vjnew=Vj-η3(∂S/∂Vj);
(16)end
(17)Replace the all Uis with Uinews, and Vjs with Vjnews, recompute S1;
(18)end
(19)returnUTV;


## 6. Experiments

In this section, experimental evaluations against a finance dataset and a traffic speed dataset are presented. The proposed MAFED algorithm is compared with the other two state-of-the-art approaches.


### 6.1. Evaluation Methodology

To evaluate the algorithms, we perform missing imputation on an incomplete matrixR. Similar to many other matrix competition problems, the testing protocol adopted is theGiven X(0<X<1) protocol [27]; that is, for R, we only show X of its observed entries to the users, while holding the remaining (1-X) observations to evaluate the trained model. For example, when X=10%,Given X means that the algorithm is trained with 10% of the nonmissing entries and the rest of 90% nonmissing ones are held and to be recovered.With theGiven X setting, (16) can be rewritten as
(21){U*,V*}=argU,V⁡min⁡∑i=1N‍∑j=1T‍(Ri,j-UiTVj)2I(Ri,j≠“⊥”)+α∑i=1N∥Ui∥22‍+β∥V1∥22+λ∑j=2T∥Vj-Vj-1∥22.
The first selected comparison algorithm is the probabilistic principle component analysis model (PPCA [28]), which achieves the state-of-the-art performance in the missing traffic flow data imputation problem [29]. The second is the probabilistic matrix factorization model (PMF), which is one of the most popular algorithms in the Netflix matrix completion problem [13].The evaluation criterion employed is the root mean square error (RMSE). More formally, for a test datasetS={s1,s2,…,sn} and the estimated set S^={s^1,s^2,…,s^n}  (s^i is the estimation of si), the RMSE of the estimation is given by (1/n)∑k=1n‍(sk-s^k)2.In all our experiments, for everyGivenX, the data partition is repeated for 5 times, and the average results as well as standard deviations are recorded.


### 6.2. Experiments on the Finance Dataset

The finance dataset is a524-by-245 matrix R, which consists of the opening prices of 524 leading U.S. companies in 245 consecutive trading days from August 21, 2009 to August 20, 2010. Each row in R corresponds to a company and each column corresponds to a day.To characterize the evolution of the opening prices, we introduce the termincremental rate; for the ith company, theincremental rate of its opening price in the jth day is given by
(22)Inci,j=Ri,j-Ri,j-1Ri,j-1.
Intuitively, theincremental rate quantifies howstrong the evolution is. A smaller |Inci,j| implies a closer connection between Ri,j-1 and Ri,j, and hence we say the evolution from Ri,j-1 to Ri,j isgradual; while a larger |Inci,j| indicates a looser relation between Ri,j and Ri,j-1, and thus the evolution is considered to besaltatory.We calculate the incremental rates for all companies in the days2~245. Figure 2 illustrates the cumulative probability distribution of these rates. We can observe that almost all the incremental rates locate in a very sharp interval (−10%, 10%). This implies that the changes of the opening prices are very slight. In other words, the evolutions of the opening prices aregradualism.Figure 2
The cumulative probability distribution curve of the incremental rates of the finance dataset.We evaluate the parametric sensitivity of the algorithm by tuning the parametersα and λ. In our experiments, we first fix α=0.01 and vary λ via an assignment expression λ=0.01×2n, where n=0,1,…,9. Then we do the reverse by fixing λ=0.01 and changing α via α=0.01×2n. For theGiven X protocol setting, we take X=50%; that is, 50% of the data in R is randomly selected as training data, with which the algorithm is trained and recovers the remaining 50% as test data. The average RMSEs of the experiments are shown in Figure 3.Figure 3
Empirical results on parameter sensitivity.As shown in Figure3, the RMSE values remain stable even when λ is expanded by more than 200 times (28=256). Similar result also appears in the experiments on the parameter α. We can observe that significant changes of the RMSE value only occur in cases where n>7 (i.e., α is expanded by more than 128 times).To study the prediction ability of the proposed algorithm, we increase theX value in theGiven X protocol from 10% to 50%. Then for each X setting, missing imputations are performed using MAFED with α=β=λ=0.5 and the comparison algorithms with the latent feature number d=10 and d=30, respectively. As illustrated in Table 1, MAFED outperforms PPCA significantly in all settings. For any X value, the RMSE of MAFED is at most 20% of that of PPCA. Specifically, for X∈{30%,40%,50%}, the RMSE of MAFED is even only 10% of PPCA. More notably, the value d has strikingly different impacts on MAFED and PPCA. When d changes from 10 to 30, almost all the RMSEs of PPCA dramatically increase, while for MAFED, RMSEs are decreased by nearly 5% except for the case with X=10%. This observation suggests that, with carefully setting of the parameter d, it is possible to improve the predicting accuracy of MAFED. We can also observe that, compared with PMF, MAFED achieves smaller RMSEs in the all settings. Especially, when X∈{30%,40%,50%}, the RMSE of MAFED is only half of that of PMF. Another interesting finding is that, for PMF, the value of d has little impact on predictions..Table 1
Imputation results on the finance dataset.
10%
20%
30%
40%
50%
D
=
10
PPCA
18.52
±
0.73
20.84
±
0.81
24.18
±
0.65
22.82
±
0.59
19.57
±
0.33
PMF
3.33
±
0.03
3.29
±
0.03
3.30
±
0.02
3.30
±
0.00
3.28
±
0.02
MAFED
3.14
±
0.01
2.92
±
0.01
2.84
±
0.02
2.32
±
0.01
2.05
±
0.02
D
=
30
PPCA
24.22
±
0.61
21.84
±
0.93
24.51
±
0.51
23.22
±
0.37
22.61
±
0.35
PMF
3.33
±
0.03
3.29
±
0.02
3.30
±
0.01
3.30
±
0.02
3.29
±
0.02
MAFED
3.11
±
0.03
2.87
±
0.02
2.70
±
0.00
2.21
±
0.01
2.09
±
0.00The convergence rate of MAFED is also explored in this work. We record the RMSEs of data recovered by MAFED every10 iterations and plot them on Figure 4, where the x-axis is the number of iterations and y-axis represents the RMSEs. We can observe that, for all X values, corresponding curves drop dramatically in the first 20 iterations and remain stable afterwards. From the result presented in Figure 3 and Table 1, we can conclude that, in our experiments, MAFED converges to the local optimization solutions within 100 iterations.Figure 4
Empirical results on the algorithm convergence.


### 6.3. Experiments on the Traffic Speed Dataset

To reconfirm the recovery performance of MAFED, we in this section conduct another evaluation on a traffic speed dataset collected in the urban road network of Zhuhai City [30], China, from April 1, 2011 to April 30, 2011. Again we adopt a matrix R to represent the data, where R consists of 1853 rows and 8729 columns. Each row corresponds to a road and each column corresponds to a 5-minute-length time interval. All columns are arranged in ascending order of time. The entry Ri,j(1≤i≤1853,1≤j≤8729) in R is the aggregate mean traffic speed of the ith road in the jth interval. As the data in R are collected by probing vehicles, the value of Ri,j might bemissing if there is no probing vehicle on the ith road during the jth time interval. In the data set, nearly half of the data, that is, 8 million entries in R, are such missing values.We calculate the incremental rates of the nonmissing data inR and present the cumulative density distribution for the incremental rates in Figure 5. The evolution of the traffic speeds differs from that of the stock opening prices remarkably. As aforementioned, the incremental rates of the finance data set mainly concentrate on the sharp interval [-10%, 10%]. While for the traffic speed data, as demonstrated in Figure 5, the range of the incremental rates spreads from -1 to 10. This indicates that there are plenty ofsudden changes of the road speeds occuring in adjacent time intervals. In other words, many of the speed evolutions aresaltation. A possible cause of this observation could be the effects of the traffic lights in the urban road network. When the traffic light of a road turns from green to yellow (and then red), traffic speed on the road immediately reduces to 0. On the other side, when the light turns from red to green, vehicles are restarted instantly, resulting in a significant acceleration of the speed on the road.Figure 5
The cumulative probability distribution curve of the incremental rates of the traffic speed dataset.This experiment is conducted with settingsd=10, α=0.5, and λ=10. Table 2 shows that in the all X settings, MAFED also outperforms PPCA and PMF. In particular, when there are few observations (e.g., X=10% and X=20%), RMSEs of MAFED are 33% lower than those of PPCA and 10% lower than those of PMF. When X>20%, the RMSE differences between PPCA and MAFED tend to be slight. Despite that, the overall errors of PPCA are roughly 3%~5% higher than those of MAFED, and for PMF, the RMSEs remain roughly 10% higher than those of MAFED.Table 2
Imputation results on the traffic dataset.
10%
20%
30%
40%
50%
PPCA
17.88
±
0.33
17.36
±
0.35
12.00
±
0.26
12.26
±
0.12
11.47
±
0.06
PMF
14.44
±
0.10
12.75
±
0.11
12.49
±
0.07
12.39
±
0.06
12.36
±
0.05
MAFED
13.30
±
0.09
12.17
±
0.10
11.92
±
0.07
11.84
±
0.06
11.04
±
0.08


## 6.1. Evaluation Methodology

To evaluate the algorithms, we perform missing imputation on an incomplete matrixR. Similar to many other matrix competition problems, the testing protocol adopted is theGiven X(0<X<1) protocol [27]; that is, for R, we only show X of its observed entries to the users, while holding the remaining (1-X) observations to evaluate the trained model. For example, when X=10%,Given X means that the algorithm is trained with 10% of the nonmissing entries and the rest of 90% nonmissing ones are held and to be recovered.With theGiven X setting, (16) can be rewritten as
(21){U*,V*}=argU,V⁡min⁡∑i=1N‍∑j=1T‍(Ri,j-UiTVj)2I(Ri,j≠“⊥”)+α∑i=1N∥Ui∥22‍+β∥V1∥22+λ∑j=2T∥Vj-Vj-1∥22.
The first selected comparison algorithm is the probabilistic principle component analysis model (PPCA [28]), which achieves the state-of-the-art performance in the missing traffic flow data imputation problem [29]. The second is the probabilistic matrix factorization model (PMF), which is one of the most popular algorithms in the Netflix matrix completion problem [13].The evaluation criterion employed is the root mean square error (RMSE). More formally, for a test datasetS={s1,s2,…,sn} and the estimated set S^={s^1,s^2,…,s^n}  (s^i is the estimation of si), the RMSE of the estimation is given by (1/n)∑k=1n‍(sk-s^k)2.In all our experiments, for everyGivenX, the data partition is repeated for 5 times, and the average results as well as standard deviations are recorded.


## 6.2. Experiments on the Finance Dataset

The finance dataset is a524-by-245 matrix R, which consists of the opening prices of 524 leading U.S. companies in 245 consecutive trading days from August 21, 2009 to August 20, 2010. Each row in R corresponds to a company and each column corresponds to a day.To characterize the evolution of the opening prices, we introduce the termincremental rate; for the ith company, theincremental rate of its opening price in the jth day is given by
(22)Inci,j=Ri,j-Ri,j-1Ri,j-1.
Intuitively, theincremental rate quantifies howstrong the evolution is. A smaller |Inci,j| implies a closer connection between Ri,j-1 and Ri,j, and hence we say the evolution from Ri,j-1 to Ri,j isgradual; while a larger |Inci,j| indicates a looser relation between Ri,j and Ri,j-1, and thus the evolution is considered to besaltatory.We calculate the incremental rates for all companies in the days2~245. Figure 2 illustrates the cumulative probability distribution of these rates. We can observe that almost all the incremental rates locate in a very sharp interval (−10%, 10%). This implies that the changes of the opening prices are very slight. In other words, the evolutions of the opening prices aregradualism.Figure 2
The cumulative probability distribution curve of the incremental rates of the finance dataset.We evaluate the parametric sensitivity of the algorithm by tuning the parametersα and λ. In our experiments, we first fix α=0.01 and vary λ via an assignment expression λ=0.01×2n, where n=0,1,…,9. Then we do the reverse by fixing λ=0.01 and changing α via α=0.01×2n. For theGiven X protocol setting, we take X=50%; that is, 50% of the data in R is randomly selected as training data, with which the algorithm is trained and recovers the remaining 50% as test data. The average RMSEs of the experiments are shown in Figure 3.Figure 3
Empirical results on parameter sensitivity.As shown in Figure3, the RMSE values remain stable even when λ is expanded by more than 200 times (28=256). Similar result also appears in the experiments on the parameter α. We can observe that significant changes of the RMSE value only occur in cases where n>7 (i.e., α is expanded by more than 128 times).To study the prediction ability of the proposed algorithm, we increase theX value in theGiven X protocol from 10% to 50%. Then for each X setting, missing imputations are performed using MAFED with α=β=λ=0.5 and the comparison algorithms with the latent feature number d=10 and d=30, respectively. As illustrated in Table 1, MAFED outperforms PPCA significantly in all settings. For any X value, the RMSE of MAFED is at most 20% of that of PPCA. Specifically, for X∈{30%,40%,50%}, the RMSE of MAFED is even only 10% of PPCA. More notably, the value d has strikingly different impacts on MAFED and PPCA. When d changes from 10 to 30, almost all the RMSEs of PPCA dramatically increase, while for MAFED, RMSEs are decreased by nearly 5% except for the case with X=10%. This observation suggests that, with carefully setting of the parameter d, it is possible to improve the predicting accuracy of MAFED. We can also observe that, compared with PMF, MAFED achieves smaller RMSEs in the all settings. Especially, when X∈{30%,40%,50%}, the RMSE of MAFED is only half of that of PMF. Another interesting finding is that, for PMF, the value of d has little impact on predictions..Table 1
Imputation results on the finance dataset.
10%
20%
30%
40%
50%
D
=
10
PPCA
18.52
±
0.73
20.84
±
0.81
24.18
±
0.65
22.82
±
0.59
19.57
±
0.33
PMF
3.33
±
0.03
3.29
±
0.03
3.30
±
0.02
3.30
±
0.00
3.28
±
0.02
MAFED
3.14
±
0.01
2.92
±
0.01
2.84
±
0.02
2.32
±
0.01
2.05
±
0.02
D
=
30
PPCA
24.22
±
0.61
21.84
±
0.93
24.51
±
0.51
23.22
±
0.37
22.61
±
0.35
PMF
3.33
±
0.03
3.29
±
0.02
3.30
±
0.01
3.30
±
0.02
3.29
±
0.02
MAFED
3.11
±
0.03
2.87
±
0.02
2.70
±
0.00
2.21
±
0.01
2.09
±
0.00The convergence rate of MAFED is also explored in this work. We record the RMSEs of data recovered by MAFED every10 iterations and plot them on Figure 4, where the x-axis is the number of iterations and y-axis represents the RMSEs. We can observe that, for all X values, corresponding curves drop dramatically in the first 20 iterations and remain stable afterwards. From the result presented in Figure 3 and Table 1, we can conclude that, in our experiments, MAFED converges to the local optimization solutions within 100 iterations.Figure 4
Empirical results on the algorithm convergence.


## 6.3. Experiments on the Traffic Speed Dataset

To reconfirm the recovery performance of MAFED, we in this section conduct another evaluation on a traffic speed dataset collected in the urban road network of Zhuhai City [30], China, from April 1, 2011 to April 30, 2011. Again we adopt a matrix R to represent the data, where R consists of 1853 rows and 8729 columns. Each row corresponds to a road and each column corresponds to a 5-minute-length time interval. All columns are arranged in ascending order of time. The entry Ri,j(1≤i≤1853,1≤j≤8729) in R is the aggregate mean traffic speed of the ith road in the jth interval. As the data in R are collected by probing vehicles, the value of Ri,j might bemissing if there is no probing vehicle on the ith road during the jth time interval. In the data set, nearly half of the data, that is, 8 million entries in R, are such missing values.We calculate the incremental rates of the nonmissing data inR and present the cumulative density distribution for the incremental rates in Figure 5. The evolution of the traffic speeds differs from that of the stock opening prices remarkably. As aforementioned, the incremental rates of the finance data set mainly concentrate on the sharp interval [-10%, 10%]. While for the traffic speed data, as demonstrated in Figure 5, the range of the incremental rates spreads from -1 to 10. This indicates that there are plenty ofsudden changes of the road speeds occuring in adjacent time intervals. In other words, many of the speed evolutions aresaltation. A possible cause of this observation could be the effects of the traffic lights in the urban road network. When the traffic light of a road turns from green to yellow (and then red), traffic speed on the road immediately reduces to 0. On the other side, when the light turns from red to green, vehicles are restarted instantly, resulting in a significant acceleration of the speed on the road.Figure 5
The cumulative probability distribution curve of the incremental rates of the traffic speed dataset.This experiment is conducted with settingsd=10, α=0.5, and λ=10. Table 2 shows that in the all X settings, MAFED also outperforms PPCA and PMF. In particular, when there are few observations (e.g., X=10% and X=20%), RMSEs of MAFED are 33% lower than those of PPCA and 10% lower than those of PMF. When X>20%, the RMSE differences between PPCA and MAFED tend to be slight. Despite that, the overall errors of PPCA are roughly 3%~5% higher than those of MAFED, and for PMF, the RMSEs remain roughly 10% higher than those of MAFED.Table 2
Imputation results on the traffic dataset.
10%
20%
30%
40%
50%
PPCA
17.88
±
0.33
17.36
±
0.35
12.00
±
0.26
12.26
±
0.12
11.47
±
0.06
PMF
14.44
±
0.10
12.75
±
0.11
12.49
±
0.07
12.39
±
0.06
12.36
±
0.05
MAFED
13.30
±
0.09
12.17
±
0.10
11.92
±
0.07
11.84
±
0.06
11.04
±
0.08


## 7. Conclusion and Future Works

Matrix factorization models are fundamental in machine learning and data mining. In this paper, we present MAFED, a matrix factorization model for evolution data. In MAFED, the two factor matrices are treated as composed of intrinsic features of the objects and time-varying features of the environment, respectively. Hence, their product accommodates the ability to describe the evolution of the objects over time. Besides, we construct a probabilistic graphical model, with which the statistical natural of MAFED is elaborated. We also present a fitting algorithm to find the solution of MAFED. Finally, we evaluate MAFED through missing data imputation on two real-world datasets. Experimental results indicate that the proposed model outperforms the comparison algorithms in both gradualism cases and saltation cases.According to the generative interpretation of MAFED, its success mainly attributes to the introduction of the linear chain structural prior; as a result, a natural extension of the present work is to deal with the factorizations on matrices with more complex structural priors. In particular, as have been shown in the empirical studies section, many of the pieces of data in real-world scenarios are highly incomplete, so it is important and interesting to study the factorization on the incomplete matrix with no prespecified structural information.



---

*Source: 525398-2014-05-25.xml*