DEIM Forum 2019 H1-1

Time Series Link Prediction using NMF
Faith MUTINDA† , Atsuhiro NAKASHIMA† , Koh TAKEUCHI†† , Yuya SASAKI† , and Makoto
ONIZUKA†† Graduate School of Information Science and Technology, Osaka University
1–5 Yamadaoka, Suita, Osaka, Japan†† NTT Communication Science Laboratories, 2–4 Hikaridai, Seika, Soraku, Kyoto, Japan
E-mail: †{wavinya.faith,nakashima.atsuhiro,sasaki,onizuka}@ist.osaka-u.ac.jp, ††koh.t@acm.org
Abstract Data in many fields such as e-commerce, social networks, and web data can be modeled as graphs, where a node
represents a person and/or an object and a link represents the relationship between people and/or objects. Since the relationships change with time, data mining techniques for time series graphs have been actively studied. In this paper, we study the
problem of predicting the links in the future graph from historical graphs. Although various studies have been carried out onlink prediction, the prediction accuracy of existing methods is still low because it is difficult to capture continuous change withtime. Therefore, we propose a new method that combines non-negative matrix factorization (NMF) and Holt-Winters method.
NMF extracts the latent features while the Holt-Winters method captures the changes of features with time. Our method canpredict hidden links that do not appear in historical graphs. Our experiments with real dataset show that our method has ahigher prediction accuracy compared to existing methods.
Key words Link prediction, Holt-Winters forecasting, Non-negative matrix factorization (NMF)

1 Introduction
The amount of data in many applications such as e-commerce,social networks, and the web is increasing rapidly. Extracting useful information from these big data can improve the quality of applications or services and can create new profits. Data from theseapplications can be modeled as graphs, where nodes represent objects and links represent the relationship between the objects. For
example, if a customer purchases an item from an e-commerce site,
we can represent this relationship using a link between two nodes,
those represent the customer and the item. Graphs can be used invarious data mining tasks, such as detecting hidden groups, detecting missing links, and ranking objects [9].
Real world data such as product sales is dynamic because their relationships change with time. Therefore, it is useful to predict how
the links in the graphs change with time. By predicting the futurelink structure, we can predict future trends and behaviors. Therefore, e-commerce sites such as amazon can recommend products
which match customers tastes and increase their sales and customer
satisfaction [8, 12, 13, 17, 18]. We formally define our problem asfollows.
Problem definition: (Time Series Link Prediction). Given thetask is to predict the set of binary links ET +1 in the future graph
GT +1 = (V1 ,V2 , ET +1 ) at time T + 1.
Various link prediction methods have been proposed, but the prediction accuracy is still low [8, 12] for time series graphs. In orderto improve the prediction accuracy, it is important to model the features of the graph properly and capture how these features changewith time. Non-negative matrix factorization (NMF) method can beused in modelling graph features because it extracts latent features
effectively and therefore helps in understanding the underlying link
structure [19]. In addition, Holt-Winters method is commonly usedin time series data forecasting because it effectively captures periodicity or seasonal fluctuations in the data [3, 10].
In this paper, we propose a new method that combines NMF and
Holt-Winters method for predicting the links in the future graphfrom historical graphs. Our method first applies NMF to extract
the latent features from the matrix representing the structure of eachhistorical graph at time t (1 <= T ) and, then, applies Holt-Winters=t <to predict the latent features at time T + 1 from the time series ofthe decomposed matrices from 1 to T . In addition, we propose two
extensions to our method so as to improve the preciseness of theprediction; matrix decomposition by leveraging time invariant features and future prediction by ensemble learning. As for the former

set of time series bipartite graphs G1 = (V1 ,V2 , E1 ), · · · , GT (V1 ,V2 , ET ) extension, we assume that given a bipartite graph with two sets offrom time 1 to T , where V1 and V2 represent the sets of nodes

nodes, the latent features in one set are time invariant, while theand Et represents the set of weighted links at each time t, thelatent features in the other set change with time. Then, we apply

NMF and Holt-Winters method to an adjacency matrix of the bi-

all the matrices from time 1 to time T , and Θ ∈ (0, 1) is a parame-partite graph, where the column or row features are held constant.

ter that is chosen by the user depending on the experiments on the

As for the latter, ensemble learning improves the accuracy of pre-training data. As shown by (3), CWT assigns more weight to thedictions [6]. By employing ensemble learning, multiple models canmost recent links. TSVD matrix decomposition method can be used

be created by changing the parameters and they are combined towith the matrices resulting from CT and CWT, we refer to them as

improve the prediction accuracy and avoid overfitting. We evaluate

TSVD CT and TSVD CWT respectively.

our method by comparing it with other methods using a real dataset

In addition, Canonical Polyadic (CP) decomposition is a common

and our method showed an improvement in the prediction accuracy.

tensor matrix decomposition method which decomposes a three di-

Further, Vector Auto Regressive (VAR) model has been success-mensional matrix to three rank-one tensors, just like SVD decom-

fully used in multivariate time series data forecasting [21, 22]. We

poses a matrix to three rank-one matrices [8, 11]. However unlike

compare the performance of Holt-Winters and VAR forecasting

SVD, the resultant tensors from CP decomposition are not orthogo-

methods, and the results show that Holt-Winters results are morenal, but they are unique and hence they can be used for forecasting.

accurate and therefore it outperforms VAR on our dataset. In spite

The CP decomposition of a tensor X of size M × N × T is definedof the popularity of VAR models, the results indicate that in someas follows:K

domains VAR does not perform well.

X ≈ ∑ λk ak ◦ bk ◦ ck ,

(4)

k=1

The remainder of this paper is organized as follows. Section 2
describes related work. Section 3 overviews the preliminaries ofwhere K is the number of components and the symbol ◦ represents

this work. Section 4 describes the details of our proposed method.

the outer product 1 . ak and bk extract the column and row features

Section 5 reviews the results of our experiments. Section 6 provides

respectively, while ck extracts the temporal components. The three

our brief conclusion.

matrices correspond to the axes of the original three-dimensional
tensor.

2 Related work

2.2 Forecasting approach

In this section, we describe existing methods for the link pre-Various time series data forecasting methods have been stud-diction problem using matrix decomposition and forecasting tech-ied [20]. Prediction of how a link changes between certain nodes

niques.

is a univariate time series prediction task, and therefore future links
can be predicted by solving all the node combinations.

2.1 Matrix decomposition approach
Truncated singular value decomposition (TSVD) is a low-rank

Exponential smoothing is a simple univariate time series predic-

matrix approximation technique which can used for time series linktion method. In the first-order exponential smoothing method, pre-prediction [7]. TSVD decomposes a matrix X of size M × N into

diction is performed on the observations y1 , · · · , yt such that more

three matrices. The best K rank approximation of the original ma-

weight is given to the most recent observation.

trix is given by
X ≈ UK ΣK VKT ,Pt+1 = α yt + (1 − α )Pt ,

(1)

(5)

where UK and VK are orthogonal matrices of size M × K and N × K

where α is a learning coefficient which takes values 0 <=α <= 1 andrespectively, and ΣK is a K × K diagonal matrix.

Pt+1 is the predicted value for one time period ahead. However,A three dimensional tensor can be reduced to a two dimensional

the first-order exponential smoothing prediction formula does not

matrix using the collapsed tensor (CT) and collapsed weighted ten-

capture change in trends such as increasing and decreasing trendssor (CWT) techniques proposed in [8]. CT method removes theeffectively.

time series information of a three dimensional tensor by taking thesum in the time direction as shown below.

such thatT

X = ∑ Zt ,In the second-order exponential smoothing method, prediction is
performed with the addition of the term bt which captures the trend

(2)

t=1

Pt+1 = α yt + (1 − α )Pt + bt and(6)bt = β (Pt − Pt−1 ) + (1 − β )bt−1 ,(7)where Zt is the matrix at time t and X is the sum of all matrices from time 1 to time T . CWT method assigns temporal weights
to the elements and reduces a three-dimensional matrix to a twodimensional matrix while maintaining the time series informationof the three-dimensional tensor as shown below.

Auto Regressive (AR) models are also widely used in data forecasting [21, 22]. The Vector Auto Regressive (VAR) model is usedT

X = ∑ (1 − Θ)T −t Zt ,where β is the trend coefficient.

(3)

t=1

where Zt is the matrix at time t, X is the final matrix after summingin prediction of future observations in multivariate time series data,and has the parameters; time lag set L, weight Al and Gaussian noise
1

Three-way outer product is defined as: X = a ◦ b ◦ c.

εt .

follow the Poisson distribution, we use the Kullback- Leibler (KL)yt = ∑ Al yt−l + εt .

(8)

l∈L

divergence. The update rules based on KL divergence are definedby the equations below as discussed by [16].

The AR model is used in various research fields, for example, in∑ xi j vk jsocial economics, Simultaneous Auto Regressive (SAR) model isj∑ uik vk jcommonly used in the analysis of spatial data since it incorporates,
uik ← uik k∑ vk jspatial auto-correlations into regression models [2]. VAR has also

(9)jbeen used in space-time prediction problems of sensor data by com∑ xik uik

bining AR model with tensor decomposition [21].

i∑ uik vk j

Recently, Recurrent Neural Networks especially Long-Shortvk j ← vk j

Term Memory (LSTM) have been studied for time series data fore-

k

∑ uik

.

(10)

i

casting [14]. This is because of their potential to capture long term

The obtained decomposed matrices can be regarded as matrices independencies in sequential data. The standard LSTM model takes awhich the features of the original matrix are reduced into a low di-sequence of vectors y1 , y2 , · · · , yt ∈ Rn as input and produces a single output vector ŷt+1 ∈ Rn , where ŷt+1 is the predicted value formension, that is, the features of each axis of the original matrix are

one time period ahead.

3 Preliminaries
As preliminaries of our method, we explain time series graph,

reduced into groups with K number of components.
3.3 Holt-Winters method
Holt-Winters is a forecasting method suitable for time series datathat has seasonality and trend [3, 10]. Seasonality means that thetime series data has trends that repeat every m cycles. There are two

NMF, and Holt-Winters method.

variations of the Holt-Winters technique; additive and multiplica-

3.1 Time series graphtive methods that differ in the seasonality component. The additive

Let G = (V1 ,V2 , E) be a bipartite graph, where V1 and V2 are setsmethod is suitable when the seasonal variations in the series areof nodes and E ∈ V1 ×V2 is a set of weighted links. The numbers ofroughly constant. On the other hand, the multiplicative method is

nodes in V1 and V2 are denoted by N and M, respectively. The bipar-

suitable when the seasonal variations change proportionally to thetite graph G is represented in a form of two-dimensional adjacency

level (average value) of the series. In this paper we choose the addi-

matrix X of size M × N. If there exists a link between nodes i andtive method because the change in our dataset is fairly constant. Thej, the (i, j) component of the matrix is assigned the weights of theadditive Holt-Winters method consists of three smoothing equations

links and zero if there is no link. We consider the graph structure

and a forecasting equation as shown below.

over a period of time. The graph structure changes with time (e.g.,

lt = α (yt − st−m ) + (1 − α )(lt−1 + bt−1 ).

(11)

bt = β (lt − lt−1 ) + (1 − β )bt−1 .

(12)

st = γ (yt − lt−1 − bt−1 ) + (1 − γ )st−m .

(13)

NMF decomposes a matrix of non-negative values to two matri-yt+h = lt + hbt + st−m+h .

(14)

ces of low dimensions such that they do not include negative val-

where y1 , y2 , · · · , yt are the observed values and m is the seasonal-ues [15]. This restriction enables NMF to provide different decom-

ity parameter which represents the length of the seasonal cycle, forposition results compared to other matrix decomposition methods

example, m = 3 for quarterly data, and m = 12 for monthly data.

such as singular value decomposition (SVD) or principal compo-

α , β , γ are smoothing parameters. lt is the smoothed estimate of thenent analysis (PCA). In addition, NMF results are easy to interpret,level at time t, bt is the smoothed estimate of the change in the trend

especially in tasks where the underlying features are interpreted as

at time t and st is the smoothed estimate of the seasonal component

non-negative.

at time t. The smoothing equations (11-13) minimize the squared

sales data in every month). We denote X (t) as an adjacency matrix
at each time t. We call the set of graphs over time range time seriesgraph.
3.2 NMF: Non-negative Matrix Factorization

In NMF, a non-negative value matrix X of size M × N is decomposed into two non-negative matrices, U and V of size M × K anderror and the forecast, (yt+h ) at h time periods ahead is calculated
as shown in (14).

K × N respectively, such that X ≈ U V . K is the base number of
NMF and is an arbitrary parameter. In this paper, we use the multi-

4 Proposed method

plicative update rules for NMF where matrices U and V are initial-

In our proposed method, we use NMF to extract the underly-

ized with random non-negative values. We use an iterative method

ing latent features and Holt-Winters method to capture the periodic

to update matrices U and V such that the divergence between thechange of the extracted latent features and predict future actions.

original matrix X and U V is minimized. In our experiments, weuse a supermarket dataset and, since sales events are expected toAlgorithm 1 Calculation of prediction matrix X (T +1)
Input Z = (X (1) , X (2) , ..., X (T ) ), XAve
Output X (T +1)
Uinit , Vinit ← XAve , Urandom , Vrandom(t)
2: for each X ∈ Z do:

1: NMF

3:

NMF:

U (t) , V (t) ← X (t) , Uinit , Vinit← U (t)

4:

Ulist

5:

Vlist ← V (t)6: end for
7: Holt − Winters:

U (T +1) ← Ulist

8: Holt − Winters:

V (T +1) ← Vlist

9: X (T +1) = U (T +1) V (T +1)

decomposition matrices Uinit , Vinit (line 1). NMF is then applied tothe matrices at each time t, with Uinit , Vinit as the initial matrices(lines 2 and 3). A list of decomposed matrices Ulist , Vlist is created
to which Holt-Winters is applied to forecast the values at time T + 1
Figure 1: Overview of proposed method: Step 1: NMF is applied to thematrices X (1) , · · · , X (T ) . Step 2: Holt-Winters method is applied to the sequence of U (t) and V (t) (1 <=t <= T ). Step 3: The predicted matrix at timeT + 1 is calculated as X (T +1) = U (T +1) V (T +1) .

(lines 7 and 8). The predicted matrix X (T +1) is obtained by the dot
product of U (T +1) , V (T +1) (line 9).
4.2 Matrix decomposition by time invariant features
In time series link prediction for bipartite graphs with two setsof nodes, we can assume that the latent features in one set are time

4.1 Outline of proposed method

invariant while the latent features in the other set change with time.

Figure1 shows the outline of the proposed method. Z is a three-

For example, since we use supermarket dataset for this paper, thisdimensional tensor which consists of the adjacency matrices ofdataset contains two sets of nodes which represent customers andthe bipartite graphs from time 1 to time T . We consider matrices

items. Then, we can assume that customers features (preferences)

X (1) , · · · , X (T ) of uniform size at each time t. NMF is applied toremain constant but products features change with time. This is be-

the matrices X (1) , · · · , X (T ) to obtain U (1) , V (1) , · · · , U (T ) , V (T ) .

cause some items such as fruits are seasonal and special products

Holt-Winters method is then applied to the sequence of U (t) and(T +1) and V (T +1) at
V (t) (1 <=t <= T ) to predict the values of U

are promoted during special events, and hence they have a greater
effect on sales. Therefore, it is important to consider a model wheretime T + 1, respectively. The predicted matrix at time T + 1 is cal-temporal change is represented in only one of the nodes of the bi-

culated as X (T +1) = U (T +1) V (T +1) .

partite graphs.

In NMF, the decomposition matrices are initialized with random

Given the adjacency matrices of the bipartite graphs from time 1

non-negative values and then the multiplicative update rules are ap-

to time T , we assume that the row features are time invariant but theplied. The final decomposition matrices depend on the initial de-column features change with time, and vice-versa. In the case ofcomposition matrices. Therefore, there is no guarantee that thetime invariant column features, NMF is applied to a matrix of size

reduced K features appear in the same order for each decomposition matrix in U (t) (V (t) ) (1 <=t <= T ). To solve this problem,

T M × N obtained by concatenating T matrices of size M × N. The

we first apply NMF to the average matrix defined as XAve (i, j) =1
Tmatrix is decomposed into two matrices U and V of sizes T M × K
and K × N respectively. The column features (matrix V ) will haveTthe same values over the entire time, and change is represented only

t=1

in the row features (matrix U ). Matrix U of size T M × K is then

∑ X (t) (i, j) to obtain Uinit , Vinit which we use as the initial ma-

trices at each time t. By using the same initial decomposition ma-

divided into T matrices, U 1 , U 2 , · · · , U T and Holt-Winters method

trix at each time t, we expect that the order of the features does

is used to forecast the values of U T +1 at time T + 1. The predicted

not change, that is, same features are likely to appear in the same

matrix at time T + 1 is calculated as X (T +1) = U (T +1) V . Figure 2

position in each decomposition matrix. This ensures that the latentfeatures are captured properly over the entire time.
Algorithm 1 shows the outline of the proposed method. The input Z is a list of same size matrices at each time t, and XAve is theaverage matrix. The output X (T +1) is the predicted matrix at timeT + 1. NMF is applied to the average matrix with random initial
non-negative value matrices Urandom , Vrandom to obtain the initialshows the proposed method with time invariant column features.
Algorithm 2 shows the algorithm for the proposed method withtime invariant column features. NMF is applied to matrix X of size
T M × N, which is then decomposed to matrices U and V of sizes
T M × K and K × N respectively (line 1). The T M × K matrix U isthen divided into T matrices to create a list of decomposition ma-

The Holt-Winters seasonality parameter, m, and the number offeatures, K, for NMF need to be selected manually. It is difficult
to search for optimum values for those parameters which give thebest performance. We employ an ensemble approach, let XKmdenote the matrix of scores calculated for K = 5, 10, · · · , 100 and
m = 1, 2, · · · , 12, the final matrix of ensemble scores is then calculated as,
X=

XKm
,
kX
Km kF
K∈(5,10,...,100) m∈(1,2,...,12)

∑

∑

(15)

where kXKm kF is the Frobenius norm for XKm .

5 Experiments
In this section, we describe experiments and their results by comparing our method proposed in Section 4 with the former methods
summarized in Section 2. We also investigate the results of our variation of 1) fixing time invariant features and 2) ensemble learning
described in Section 4.
In our experiments we use a point-of-sales dataset for a period of
24 months obtained from a supermarket. Supermarket data is periFigure 2: Proposed method by leveraging time invariant features V

odic because some products such as vegetables are seasonal and special products are promoted during special events such as Christmas

Algorithm 2 Calculation of prediction matrix X (T +1) by leverag-and Valentines. The dataset contains 25668 customers and 113688

ing the time invariant column features

items. From this dataset we extract the top 1000 frequent customers

Input Z = (X (1) , X (2) , ..., X (T ) )and top 500 items with highest number of sales. We transform the

Output X (T +1)

data to adjacency matrices of bipartite graphs, and use the data for

1: NMF:U, V ← X

2: Time Divide

U (1) , U (2) , ..., U (T ) ← U:

3: Ulist ← U (1) , U (2) , ..., U (T )4: Holt − Winters:

U (T +1) ← Ulist

5: X (T +1) = U (T +1) V

23 months as training dataset and the data for the last month as thetest dataset. We generate a three-dimensional tensor Z extended in
the time direction, that is, if customer i buys n items in month t, thenZ(i, j,t) = n. We eliminate the influence of large values in the data
by normalizing the data according to the equation below which was
proposed by [8].

trices. Holt-Winters is applied to the list of matrices to predict thevalues of U T +1 at time T + 1, (line 4). The predicted matrix at timeT + 1 is calculated as X T +1 = U T +1 V (line 5).
Assuming that the temporal change in one set of the nodes of a
1 + log(n) n > 0
Z(i, j,t) =0

n = 0.

(16)

bipartite graph is significantly larger than that in the other node, thisThe link prediction problem can be regarded as a matrix comple-

algorithm improves modelling of the latent features and prediction

tion problem with implicit feedback [1]. Given a M × N adjacency

accuracy.

matrix, we assume that the matrix is binary where the presence of

4.3 Holt-Winters and NMF methods with ensemble learning

links is represented by 1s and absence of links is represented by 0s.

Ensemble learning technique involves strategically combining

Since the objective of our study is to predict whether there is a linkseveral models so as to improve the stability and predictive perfor-or not between the customers and items, if customer i purchased

mance. Previous studies on ensemble learning show that combining

item j n times, the link information is represented as,1 n > 0
Y (i, j) =0 n = 0.

multiple individual models improves prediction accuracy compared
to using a single model [6]. The idea of combining multiple models
assumes that it is difficult for a single model to understand the underlying structure, but multiple models can capture different aspects(17)of data. Ensemble learning is helpful when it is difficult to choose

In addition, the Holt-Winters smoothing parameters α , γ , β and

optimum values of parameters, and when one wants to avoid large

seasonality component need to be decided before-hand. We use theerrors [4].

approach discussed in [3] to establish the values for the smoothing
parameters which minimize the sum of squared errors for one time

step ahead forecast. We set the seasonality component to m = 3

Table 1: Comparison with existing methods: all links prediction

months, by assuming seasonality appears every quarter year in sales

Prop.

Base-

CP+

default

line

HW

TN

393342389048

data. Moreover, Holt-Winters method can also be used with thetemporal information extracted by CP (described in Section 2) to

CP

HW

TSVD

376670373733409894369646
97121

CWT

FP

73425

77719

90097

93034

56873

FN

3924

2764

4623

4587

5203

4887

TP

29309

30469

28610

28646

28030

28346

described in Section 2; CP+HW, CP and TSVD CWT. For TSVD,

P.L

10273410818811870712168084903

125467CP and NMF decomposition, instead of using a fixed value of K,

C.P.L

29309

30469

28610

28646

28030

28346

we use an ensemble approach as shown below.

Prec.

28.53

28.16

24.10

23.54

33.01

22.59

Rec.

88.19

91.68

86.09

86.20

84.34

85.29

F-meas.

43.11

43.09

37.66

36.98

47.45

35.72

predict future observations, we refer to this as CP+HW. We compare the performance of our proposed methods with the methods

XK
X=
,
∑
kX
K kF
K∈(5,10,...,100)(18)

where XK is the matrix of scores calculated for K = 5, 10, · · · , 100,

Table 2: Comparison with existing methods: hidden links Prediction
Prop.

X is the final matrix of ensemble scores and kXK kF is the Frobenius norm for XK .

Baseline

default
P.L

CP+

CP

HW

TSVD

26731

0

36275

HW

25695

0

27343

CWT

In addition, we select a threshold value based on the Youden’s J

C.P.L

689

0

589

559

0

667

statistic 2 and use this threshold to determine the presence or ab-

Prec.

2.68

0

2.15

2.09

0

1.84

Rec.

54.77

0

46.82

44.44

0

53.02

F-meas.

5.11

0

4.12

3.99

0

3.55

sence of a link [5]. That is, predicted values above the threshold areconsidered as positive links and predicted values below the threshold indicate that there is no link between the corresponding customer and item. Further, we construct the receiver operating characteristic (ROC) curve and the area under the curve (AUC) measures the discrimination, that is, the ability to predict true positives
and true negatives correctly. We compare the performance of themethods in predicting all links and hidden links. In all links prediction, we compare how the different methods perform in predicting
the links in the test dataset. On the other hand, the prediction of
the hidden links addresses a difficult task, that is, how the different methods predict links that do not previously exist in the training
dataset. Since our objective is to predict the significant entries, we
treat all nonzero entries as ones (positive links) and the other entries
as zeros (no links). The total number of links in our test dataset was
33233 and the number of hidden links was 1258.
5.1 Comparison with existing methods
We compare the performance of the methods in terms of their precision, recall, and F-measure. We use Prop. default, TN, FP, FN, TP,
P.L, C.P.L, Prec., Rec. and F-meas. to denote proposed method, thenumber of true negatives, false positives, false negatives, true positives, predicted links, correctly predicted links, precision, recall,and F-measure, respectively. In addition, we use simple average as
the baseline method.
Table 1 shows the performance in the prediction of all links. HoltWinters method has the highest accuracy and F-measure. Excluding the Holt-Winters method, our method has the highest precision
and F-measure. The average prediction performs well because we
considered the most frequent users and items in our dataset and forecasts for only one time step ahead. However, for multiple-step ahead
forecasting and infrequent items, the performance reduces significantly because average method cannot capture data patterns.

Table 2 shows the results of the methods in predicting hidden
links. As expected, the F-measures for hidden links prediction arelower than the F-measures for all links prediction. The results showthat the Holt-Winters and average prediction methods cannot predict
hidden links at all and their F-measure is 0. Our proposed method
achieved the highest accuracy and F-measure values. This result
indicates that the proposed method is better in extracting hiddenfeatures by analyzing the underlying graph structure.
Figure 3 shows the performance of the methods in terms of thearea under the receiver operating characteristic curve (AUC) for all
links prediction. ROC curve shows the trade-off between sensitivity and specificity, and therefore curves close to the top-left cornerindicate a better performance [1]. The methods perform well because their ROC curves are close to the top-left corner. Holt-Winters
method is the best when the false positive rate is low. As the false
positive rate increases, the proposed method is the best and has thehighest AUC value.
5.2 Effectiveness of leveraging time invariant features
In Section 4 we described the method of matrix decomposition byleveraging time invariant features to improve accuracy of the link
prediction task. We use Prop. column and Prop. row to denotevariants of our proposed method by leveraging the time invariant
column features and the row features, respectively.
When predicting all links, table 3 shows that our method of leveraging time invariant row features achieved the best performance interms of precision and F-measure. Table 4 shows the performancein predicting hidden links. Our method of leveraging time invariant
row features achieved the highest precision and F-measure. These
results indicate that there is little change in row features with time
and there are large temporal changes in column features. In our

2

J = max(sensitivity + speci f icity − 1)

Table 3: Comparison with methods of time invariant features: all links pre-Table 4: Comparison with methods of time invariant features: hidden links

dictionprediction
Prop.

Prop.

Prop.

Prop.

Prop.

default

columnrow

default

columnrow

TN393342392518389586

25695

29120

23395

FP

73425

74249

75040

C.P.L

689

668

638

FN

3924

3705

4820

Precision

2.68

2.29

2.73

P.L

Prop.

TP

29309

29528

30554

Rec.

54.77

53.10

50.72

P.L

102734103777105594F-meas.

5.11

4.40

5.18

C.P.L

29309

29528

30554

Prec.

28.53

28.45

28.94

Rec.

88.19

88.85

86.37

F-meas.

43.11

43.10

43.35

dataset, rows represents customers while columns represent items.
The results show that the customers tastes remained constant while
the product features changed with time. This can be explained by
the fact that some products in supermarkets, especially food products such as vegetables change depending on seasons. However,
only a single dataset was used for our experiments, and the performance of each method may vary depending on datasets. Therefore,
by analyzing the characteristics of the dataset, it is possible to selectthe appropriate method and hence improve the accuracy.
Figure 4 shows the ROC curves for our method of leveraging timeinvariant features in all links prediction. Although the ROC curve of

Figure 4: ROC curves for our methods by leveraging time invariant features

the methods overlap largely, the AUC value of the method of matrix
decomposition by time invariant row features is slightly higher thancluding Holt-Winters method, the method of leveraging time invari-that of the proposed method and, therefore, it has the best perfor-

ant row features has the highest accuracy. Comparing the resultsmance in terms of AUC value.

of table 1, table 3 and table 5, we observe that ensemble learning

5.3 Comparison with method of ensemble learning

method led to a significant increase in the prediction accuracy of

In section 4 we described an extension to our proposed method

the proposed methods and Holt-Winters method. In hidden links

by applying ensemble learning. Here we vary the Holt-Winters sea-

prediction, table 6 shows that the proposed method achieved thesonality parameter, m and K for NMF as shown by equation (15)highest precision and F-measure. Although the F-measure is low,

described in section 4. Table 5 shows that the Holt-Winters method

it increased significantly from 5.11 to 5.86. Since we considered a

has the highest accuracy and F-measure in all links prediction. Ex-

single dataset in our experiments, the performance of the proposed
method and the method of matrix decomposition by time invariant
features might vary. Therefore, depending on the dataset one canselect the method with the best performance.

Table 5: Comparison with methods of ensemble learning: all links prediction

Figure 3: ROC curve comparison with existing methods

Prop.

Prop.

Prop.

Base-

CP+

default

column

row

line

HW

HW

TN

394418391455394137389048385875432022

FP

72349

75312

72630

77719

80892

34745

FN

3781

3215

3442

2764

4084

6832

TP

29452

30018

29791

30469

29149

26401

P.L

101801

105330102421108188110041

61146

C.P.L

29452

30018

29791

30469

29149

26401

Prec.

28.93

28.50

29.09

28.16

26.49

43.18

Rec.

88.62

90.33

89.64

91.68

87.71

79.44

F-meas.

43.62

43.33

43.92

43.09

40.69

55.95

Table 6: Comparison with methods of ensemble learning: hidden links pre-experiments, we confirmed that the proposed methods perform well,dictionespecially in predicting hidden links.
Prop.

Prop.

Prop.

Base-

CP+

HW

There are several directions which can be considered for future
work. First, in this paper we considered only a sales dataset, it isdefault

column

row

line

HW

P.L

19227

21897

21636

0

24210

0

C.P.L

600

555

636

0

545

0

Prec.

3.12

2.53

2.94

0

2.25

0

Rec.

47.69

44.12

50.56

0

43.32

0

information, therefore attribute-based prediction is important to fur-F-meas.

5.86

4.79

5.56

0

4.28

0

ther improve the accuracy of our methods.

important to evaluate the performance of our methods with different datasets. Second, our proposed method does not use attribute

References

Figure 5: ROC curves for methods of ensemble learning

Figure 5 shows the ROC curves for this method in all links prediction. As expected, the AUC values of the proposed methods increased compared to the values in Figure 4. All methods perform
well because the ROC curves are close to the top-left corner. The
proposed method has the highest AUC value.
In addition, we carried out additional experiments with LSTM
and VAR forecasting methods described in section 2. Although previous studies show that LSTM and VAR perform well in time series
forecasting, they did not outperform the proposed method on ourdataset. In addition, our dataset contains monthly sales data for aperiod of 24 months, the time period is short, making LSTM unsuitable. Further, LSTMs are non-linear models with different architectures, identifying a suitable architecture is challenging compared toour proposed method.

6 Conclusion
In this paper we addressed the time series link prediction problem. We proposed a method of extracting latent features from timeseries data and modelling the periodicity of the features by combining NMF and Holt-Winters methods. We also proposed extensions
to the proposed method, through matrix decomposition by time invariant features and applying ensemble learning to NMF and HoltWinters methods to improve the prediction accuracy. We compared
the performance of our methods with existing link prediction methods in predicting existing links and hidden links. As a result of the[1] Charu C Aggarwal et al. Recommender systems. Springer, 2016.
[2] Luc Anselin and Anil K Bera. Introduction to spatial econometrics.
Handbook of applied economic statistics, pages 237–290, 1998.
[3] Chris Chatfield and Mohammad Yar. Holt-winters forecasting: some
practical issues. The Statistician, pages 129–140, 1988.
[4] Robert T Clemen. Combining forecasts: A review and annotated bibliography. International journal of forecasting, 5(4):559–583, 1989.
[5] Predictivities Based ROC Curve. Medicalbiostatistics. com.
[6] Thomas G Dietterich. Ensemble methods in machine learning. In
International workshop on multiple classifier systems, pages 1–15.
Springer, 2000.
[7] Susan T Dumais, George W Furnas, Thomas K Landauer, Scott Deerwester, and Richard Harshman. Using latent semantic analysis toimprove access to textual information. In SIGCHI, pages 281–285,
1988.
[8] Daniel M Dunlavy, Tamara G Kolda, and Evrim Acar. Temporal link
prediction using matrix and tensor factorizations. ACM Transactions
on Knowledge Discovery from Data (TKDD), 5(2):10, 2011.
[9] Lise Getoor and Christopher P Diehl. Link mining: a survey. ACM
SIGKDD Explorations Newsletter, 7(2):3–12, 2005.
[10] Prajakta S Kalekar. Time series forecasting using holt-winters exponential smoothing. Kanwal Rekhi School of Information Technology,
4329008:1–13, 2004.
[11] Tamara G Kolda and Brett W Bader. Tensor decompositions andapplications. SIAM review, 51(3):455–500, 2009.
[12] Yehuda Koren. Collaborative filtering with temporal dynamics. In
SIGKDD, pages 447–456. ACM, 2009.
[13] Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization
techniques for recommender systems. Computer, (8):30–37, 2009.
[14] Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu.
Modeling long-and short-term temporal patterns with deep neuralnetworks. In ACM SIGIR, pages 95–104, 2018.
[15] Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects
by nonnegative matrix factorization. Nature, 401:788–791, 1999.
[16] Daniel D. Lee and H. Sebastian Seung. Algorithms for non-negative
matrix factorization. In NIPS, pages 556–562. MIT Press, 2001.
[17] David Liben-Nowell and Jon Kleinberg. The link-prediction problem
for social networks. JASIST, 58(7):1019–1031, 2007.
[18] Greg Linden, Brent Smith, and Jeremy York. Amazon. com recommendations: Item-to-item collaborative filtering. IEEE Internet
computing, (1):76–80, 2003.
[19] Xin Luo, Mengchu Zhou, Mingsheng Shang, Shuai Li, and Yunni
Xia. A novel approach to extracting non-negative latent factors from
non-negative big sparse matrices. IEEE access, 4(1), 2016.
[20] Spyros Makridakis and Michele Hibon. The M3-competition: results, conclusions and implications. International journal of forecasting, 16(4):451–476, 2000.
[21] Koh Takeuchi, Hisashi Kashima, and Naonori Ueda. Autoregressive
tensor factorization for spatio-temporal predictions. In ICDM, pages
1105–1110, 2017.
[22] Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon. Temporal regularized matrix factorization for high-dimensional time series prediction.
In NIPS, pages 847–855, 2016.