Tensor Decomposition with Missing Indices

1 Introduction

A tensor data is a set of values indexed by multiple indices. Some kinds of real world data are naturally represented as tensors, such as the number of times a user posts a word at a specific location in social media. Tensor decomposition, such as CP decomposition [Carroll and Chang, 1970; Harshman, 1970], is a fundamental tool to analyze the tensor data, which has various kinds of applications including recommendation [Rendle and Schmidt-Thieme, 2010], ranking semantic Web data [Franz et al., 2009], and many others [Kolda and Bader, 2009].

One of the major applications of tensor decomposition is the estimation of missing values. Suppose we have a movie rating dataset, which contains a list of quadruplet (user, movie, context, score) that indicates a user watched a movie at a context and evaluated its quality as score. As a tensor, (user, movie, context) represents its index and score represents its value. In general, there are many missing values, i.e., the number of observed samples (quadruplets) is much smaller than the number of all the possible combinations of users, movies, and contexts. By estimating the parameters from the observed quadruplets, tensor decomposition can predict the missing values.

In real applications, not only values but also a part of indices can be missing. Let us take a Twitter data as an example, which consists of a set of tweets and each tweet contains a lot of attributes such as a user ID, a message text, a timestamp, and a location. This data can be represented as a higher-order tensor and analyzing it by tensor decomposition looks promising. However, in Twitter, some attributes are optional and they are not always observed, i.e., some indices are missing. Indeed, only 0.42% of tweets contain the location attribute [Cheng et al., 2010]. We call a sample (a pair of indices and a value) whose indices are partially missing an incomplete sample, and call a sample whose indices are completely observed a complete sample.

Existing tensor decomposition methods are not directly applicable to the tensors with missing indices. A straightforward way to handle the tensor data containing both complete and incomplete samples is just ignoring the incomplete samples. That is, if we observe that user i posts word j at some unknown location, we simply discard this information. However, when the majority of samples are incomplete, this approach is problematic-it may cause severe degradation of solution due to the decrease of the sample size. Another way of handling missing indices is to discard the index through all samples. That is, we discard the location attributes of all tweets if there is at least one tweet that does not have its location attribute. This approach is again not desirable-we cannot use the location attributes at all even though some tweets have those information.

Present work. In this paper, we study the tensor decomposition problem where indices are partially missing. For this problem, we propose a probabilistic generative model for tensor decomposition that handles missing indices. Our main idea is to deal with the missing indices as latent variables. By inferring the latent variables, we can leverage incomplete samples as well as complete samples, to learn tensor decomposition. One of the advantages of the proposed model is its flexibility-it can be applied to the higher-order (>3) tensors, and to any tensor decomposition models (e.g., CP decomposition). For parameter inference, we develop an algorithm based on the variational MAP-EM algorithm [Gupta et al., 2011]. We perform experiments on both synthetic and real datasets. The experimental results show that utilizing both complete and incomplete samples leads to improving the accuracy of tensor decomposition, comparing the proposed model and baselines that do not appropriately deal with incomplete samples.

Contributions. The contributions of this paper are summarized as follows:

1. New Problem: we define a tensor decomposition problem where the indices are partially missing.
2. Model: we propose a probabilistic generative model for tensor decomposition, which handles the missing indices as latent variables. The proposed model is flexible, i.e., it can be applied to higher-order (>3) tensors, and to any tensor decomposition models.
3. Algorithm: we develop a parameter inference algorithm for the proposed model based on the variational MAPEM algorithm.

2 RelatedWork

In the history of the tensor decomposition research [Kolda and Bader, 2009], many decomposition methods have been proposed, such as CP decomposition [Carroll and Chang, 1970; Harshman, 1970], Tucker decomposition [Tucker, 1966], pairwise-interaction decomposition [Rendle and Schmidt-Thieme, 2010], and RESCAL [Nickel et al., 2011], to name a few. Since tensors are often very large, one of the main streams of the tensor decomposition research is about efficient computation, including sampling-based methods [Papalexakis et al., 2015], online/streaming computations [Sun et al., 2006; Maehara et al., 2016], memory efficient algorithms [Kolda and Sun, 2008; Oh et al., 2017], and parallel computations [Jeon et al., 2016].

Tensor data in real world tends to have missing values, which has triggered a lot of studies that tackle with them. [Acar et al., 2011a] developed an extension of CP decomposition that can handle missing values, which is solved as a weighted least square problem. [Jain and Oh, 2014] theoretically studied how many samples are needed to exactly reconstruct a data tensor. Similar to our study, some papers [Chu and Ghahramani, 2009; Rai et al., 2014; Zhao et al., 2016] proposed probabilistic models to address missing values. However, the problem these studies address is completely different from our problem, because their focus is on the missing values, whereas our focus is on the missing indices.

Although, to the best of our knowledge, there is no study to address the missing indices, the most similar one to ours would be coupled matrix and tensor factorizations (CMTF) [Acar et al., 2011b]. CMTF decomposes a set of tensors and matrices that share some of their modes. Although CMTF is designed for utilizing side information matrices in addition to the target data tensor, it is practically applicable to our problemwhere some indices are missing. That is, it could utilize incomplete samples by representing them as matrices, and decompose them together with the tensor that represents the set of complete samples. However, in this way, CMTF

simply ignores the missing indices. In contrast, our proposed model infers the missing indices and utilize them for improving tensor decomposition. We compare these methods in the experiments.

3 Setup

In this section, we introduce the notation and define the problem. To make the notation and the presentation easy to follow, we only discuss the case where we are given a 3rd-order tensor. However, we emphasize that our model can be naturally generalized to higher-order tensors.

3.1 Conventional Tensor Decomposition

For n ∈ N, we write by [n] the set {1, 2, . . . , n}. Let X be an I × J × K tensor that can contain missing values, where I, J,K ∈ N, and let D ⊆ [I] × [J] × [K] be a set of indices of observed values. We denote by ∅ the missing values, i.e., Xijk = ∅ if (i, j, k) /∈ D, or Xijk ∈ R otherwise. Using the squared loss, tensor decomposition of X is defined as follows:

Definition 1 (Squared-loss tensor decomposition). Let P be a space of parameters and ˆ Xijk : P → R be a predictive model. Given X ∈ {R ∪ {∅}}I×J×K , tensor decomposition with ˆ X is to compute ˆΘ = arg min Θ∈P L(X ; Θ) where

There are several choices for the predictive model ˆ X. Here we define CP decomposition model that is used as the example in this paper:

Definition 2 (Rank-R CP decomposition). Given R ∈ N, let Θ = {U, V,W} where U ∈ RI×R, V ∈ RJ×R, and W ∈ RK×R. The predictive model is defined as

Hereafter we write ˆ Xijk instead of ˆ X(Θ)ijk if it is not confusing.

3.2 Tensor Decomposition withMissing Indices

Next we consider the case where a part of indices are missing. In such a case it is difficult to represent data as a tensor. Instead, we consider the data as a set of indices and values. Suppose we have N ∈ N observations (i.e., samples) and denote by X = {xn}Nn=1 the values and by ˆ Z = {(ˆin, ˆjn, ˆ kn)}Nn=1 the corresponding indices that can be missing. As well as missing values, we denote the missing index by ∅.

The tensor decomposition problem with missing indices is then defined as follows:

Problem 1 (Tensor decomposition with missing indices). Suppose we have a set of values X ∈ RN , a set of indices ˆ Z ∈ {[I] ∪ {∅}} × {[J] ∪ {∅}} × {[K] ∪ {∅}}N , a tensor decomposition model (P, ˆ X), and an oracle O that returns the true indices for all missing indices. We want to find

ˆΘ = arg min Θ∈P L(X, Z; Θ), where Z = O( ˆ Z) is a set of true

indices fulfilled by the oracle and

4 ProposedModel

Since the oracle in Problem 1 is in reality not in hand, we solve the relaxed problem by estimating set of true indices Z to obtain the minimizer ˆΘ. To handle missing indices, we consider them as latent variables. Let in ∈ [I] be the latent variable denoting the true index which is not given, and ˆin ∈ {[I] ∪ ∅} be the given index, which can be observed or missing. We hypothesize the latent variable in is generated from the following distribution:

where δ(a, b) is the delta function taking 1 if a = b, or 0 otherwise, and Unif (1, I) is the categorical uniform distribution from 1 to I. pJ and pK are also defined in the same way. The above definition means that, in prior, we have no information about true index in if given index ˆin is missing, whereas true index in is always the same as the given index if it is observed. In our model, we infer those latent variables (i.e., missing indices) and learn tensor decomposition simultaneously, which improves the accuracy of tensor decomposition as will be shown in the experiments.

The generative process of our proposed model is written as follows:

where N(·|µ, σ2 ) is Gaussian distribution with mean µ and variance σ2 , and pΘ(·|λ) is the prior distribution of Θ with some parameter λ depending on the tensor decomposition model. Fig. 1 shows the graphical model of the proposed model. Note that the proposed model allows the duplicated indices, that is , there exist samples n1 and n2 such that ˆin1 = ˆin2 , ˆjn1 = ˆjn2 , and ˆ kn1 = ˆ kn2 .

The following proposition states that, if we could estimate the true indices correctly for all n, the MLE of Θ in this model is identical to the solution of Problem 1:

Proposition 3. If the estimted indices ¯ Z are the same as the true indices Z, the following holds:

Proof. We can write the log-likelihood as follows:

where C is some constant, which concludes the proof.

With regard to the above proposition, we can say that we can solve Problem 1 at least approximately by estimating the missing indices and Θ simultaneously.

To prevent from the overfitting, instead of MLE, we aim to obtain the marginal MAP as follows:

Since p(Θ|X, ˆ Z, λ) ∝ p(X, ˆ Z,Θ|λ), the marginal MAP can be obtained maximizing the joint distribution of the model written as follows:

where

4.1 Parameter Inference

In this section, we develop a parameter inference algorithm based on the variational MAP-EM algorithm. The joint distribution defined in the previous section can be written in the variational form as follows:

where

and q is any distribution over Z. Since KL(q||p) > 0, the joint distribution is always greater than L[q,Θ] that is called variational lower bound (VLB). In the variational MAP-EM algorithm, instead of directly maximizing the joint distribution, the VLB is maximized to obtain the marginal MAP by alternatively performing E-step and M-step as described in the following paragraphs.

In the E-step, L[q,Θ] is maximized with regard to q(Z) with Θ fixed. Assuming the mean field approximation q(Z) = ∏n q(in)q(jn)q(kn), and differentiating L[q,Θ] by q(in = i), we get the following updating equation:

where an,i = - 1 2 Eq(jn,kn) [(xn - ˆ Xijnkn )2 ] . We can also update q(jn) and q(kn) in the same way. Since q(in), q(jn), and q(kn) depend on each other, we repeatedly calculate them until they converge. Note that if ˆin = ˆi (i.e., ˆin is not missing), q(in) is instantly obtained as δ(in,ˆi) because p(in|ˆin) = δ(in,ˆi), which means that there is no need to calculate q for observed indices.

In the M-step, L[q, Θ] is maximized with regard to Θ with q(Z) obtained in the previous E-step fixed. Differentiating L[q, Θ] by θ ∈ Θ, we get:

where we denote zn = (in, jn, kn).

The concrete computations of the E-step and the M-step depend on the tensor decomposition model, in other words, depend on how to define ˆ Xijk and p(Θ|λ). In the next subsection, we describe the parameter inference algorithm for CP decomposition as an example. However, we emphasize again that any tensor decomposition models can be plugged-in to our proposed model.

Parameter inference for CP decomposition

Recall that the parameters of CP decomposition are Θ = {U, V,W}. Hereafter we write E [·] for Eq [·] if it is not confusing. In the E-step, we can calculate an,i as follows:

where ⊙ denotes the Hadamard product, and we omit the constant terms that are not needed to compute q(in).

In the M-step, pΘ(Θ|λ) needs to be defined. For CP decomposition, we define pΘ(Θ|λ) as follows:

where I denotes the identity matrix with the appropriate number of dimensions. Keeping in mind the above definition of pΘ(Θ|λ), we differentiate L[q,Θ] by Ui, and get the following updating equation:

V and W can also be updated in the same way. Since U, V , and W depend on each other, we alternatively update them until they converge. The sufficient statistics of the algorithm are E [Uin ], E [Vjn ], E [Wkn ], E [Uin U T in ], E [Vjn V T jn ], and E [Wkn W T kn ] for all n.

To summarize, the parameter inference algorithm for the proposed model with CP decomposition is shown in Algorithms 1-3. Every time q and Θ are updated, the related sufficient statistics have to be also updated.

Complexity. We show the time complexity of the proposed algorithm for higher-order tensors. Suppose we have Mth-order tensor where the number of dimensions of mth mode is Im. The time complexity of the proposed algorithm with CP decomposition model for each iteration is O(R3 ∑m Im +R2 (N+ ∑m N-mIm)), where N - m is the number of missing indices for mth mode. This complexity becomes the same as the conventional CP decomposition when there is no incomplete samples (i.e., N - m = 0 for all m).

5 Experiments

In this section, we report the results of the experiments conducted on both synthetic and real datasets. We compare the three variants of the proposed model, and three baselines as follows:

• MAP-EM (Proposed): Algorithm 1 with q inferred by the data.
• Uniform (Proposed): Algorithm 1 with q fixed as the categorical uniform distribution.
• Prior (Proposed): Algorithm 1 with q fixed as the distribution that is estimated from the data as a histogram. It is obtained as q(in = i) = N + i /N + I for all n where N + i is the number of samples where ˆin = i, and N + I is the number of samples where ˆin ̸= ∅.
• Minimal: CP decomposition using only complete samples.
• Complete: CP decomposition using only modes without missing indices. For example, when the 3rd index of the 3rd-order tensors can be missing, this method performs matrix decomposition using only 1st and 2nd modes.
• CMTF: coupled matrix and tensor factorization proposed in [Acar et al., 2011b]. This method can utilize a set of incomplete samples by representing it as a matrix.

The first three are our proposed models with q set in different way.

Reproducibility. Our code to reproduce the experiments is available at https://goo.gl/W4K86I. The real dataset used in the experiments in Section 5.2 is publicly available.

5.1 Synthetic Dataset

Settings. To see on what condition the parameters of the proposed model is correctly recovered, we generate two data tensors with a large/small number of samples. The size of the

Algorithm 1 Variational MAP-EM algorithm for the proposed model with CP decomposition (3rd-order)

Input: R, λ. Output: U, V , W. 1: Initialize q, U, V , andW. 2: Calculate all sufficient statistics. 3: while convergence criterion is met do 4: # E-step 5: for n = 1 to N do 6: while convergence criterion is met do 7: if ˆin = ∅ then 8: Update-q(in) (Alg. 2) 9: Update E [Uin ] and E [Uin U T in ] 10: end if 11: if ˆjn = ∅ then 12: Update-q(jn) (Alg. 2) 13: Update E [Vjn ] and E [Vjn V T jn ] 14: end if 15: if ˆ kn = ∅ then 16: Update-q(kn) (Alg. 2) 17: Update E [Wkn ] and E [Wkn W T kn ] 18: end if 19: end while 20: end for 21: # M-step 22: while convergence criterion is met do 23: Update-U (Alg. 3) 24: for each n where ˆin = ∅ do 25: Update E [Uin ] and E [Uin U T in ] 26: end for 27: Update-V (Alg. 3) 28: for each n where ˆjn = ∅ do 29: Update E [Vjn ] and E [Vjn V T jn ] 30: end for 31: Update-W (Alg. 3) 32: for each n where ˆ kn = ∅ do 33: Update E [Wkn ] and E [Wkn W T kn ] 34: end for 35: end while 36: end while

Algorithm 2 Update-q(in)

Output: q(in) 1: z ← 0 2: for i = 1 to I do 3: Calculate an,i (Eqn. 3). 4: z ← z + exp(an,i). 5: end for 6: for i = 1 to I do 7: q(in = i) ← exp(an,i)/z. 8: end for

tensors is 10 × 10 × 10. First, we generate each element of U true ∈ R10×R, V true ∈ R10×R, and W true ∈ R10×R from Gaussian distribution N(·|0, λ-1 ) with λ = 1.0 and R = 3. Then, we randomly divide the set of all indices D = [10] × [10] × [10] into training set Dtrain (90%) and test set Dtest (10%). For each index in Dtrain, we generate S samples from Gaussian distribution N(·|Xijk, 1), where Xijk = ∑r U true ir V true jr W true kr . S = 1 for small dataset,

Algorithm 3 Update-U

Output: U 1: for i = 1 to I do 2: Gi ←O # zero matrix of size R×R 3: Hi ← 0 # zero vector of size R 4: end for 5: for n = 1 to N do 6: z ←E [Vjn V T jn ] ⊙ E [Wkn W T kn ] 7: zz ←E [Vjn V T jn ] ⊙ E [Wkn W T kn ] 8: if ˆin = ∅ then 9: for i = 1 to I do 10: Hi ← Hi + q(in = i)xnz 11: Gi ← Gi + q(in = i)zz 12: end for 13: else 14: Hˆin ←Hˆin +xnz 15: Gˆin ←Gˆin +zz 16: end if 17: end for 18: for i = 1 to I do 19: Ui ← (Gi +λI)-1 Hi 20: end for

whereas S = 10 for large dataset. For each index in Dtest, we use Xijk itself without Gaussian noise.

We compare the cases where the different ratio of indices are missing. After generating the data tensors, we delete the 3rd index by probability α. We vary α from 0 to 0.9 by the step size 0.1. For each α, we perform 5 trials and report the mean and the standard deviation of the results. Note that since Complete discards the 3rd mode, it predicts the test samples as matrix decomposition, i.e., ˆ Xijk = ∑r UirVjr for all k. We use the root mean square error (RMSE) defined as follows: RMSE =

Results. Figs. 2(a) and 2(b) show the results. We do not compare Prior because the prior is the categorical uniform distribution in this experiment, which is the same as Uniform. On the large dataset (Fig. 2(a)), MAP-EM achieves the smallest RMSE among all compared methods, which means that the parameters are correctly recovered. On the other hand, on the small dataset (Fig. 2(b)), Uniform is better than MAP-EM when α ≥ 0.6, which indicates that the distribution q is not correctly inferred when the number of samples is not enough and the missing ratio is too large. CMTF and Minimal are almost always worse than MAP-EM and Uniform, meaning the proposed model works reasonably well for our problem where some indices are missing. Complete does not work at all when the missing ratio is not so large, because it discards all the information about the 3rd mode when at least one sample is incomplete. However, it shows comparable (or even better) results to other methods when the number of samples is too few and the missing ratio is too large, because, in such a case, the provided data is almost a matrix (i.e., most of 3rd indices are missing).

5.2 Real Dataset

Settings. We use the Twitter dataset that is publicly available. 1 This dataset consists of about 13M geotagged tweets. We construct a user-hashtag-location (3rd-order) tensor as follows. First, we discard all the tweets posted from outside the US. The number of remained tweets posted from inside the US is about 2.85M. Then we extract user ids, hashtags, and locations from those tweets. We use 49 states that appear in the dataset as locations. We filter user ids that occur less than 5 times, and hashtags that occur less than 30 times. As a result, we construct a 25297 × 2616 × 49 tensor X where Xijk = 1 if user i posts hashtag j at location k, or Xijk = 0 otherwise. The number of non-zero elements is 125, 017. Since this data is very sparse, we randamly sample the same number of zero elements as non-zero elements, that is, the number of observed (zero and non-zero) samples is N = 250, 034. We randomly divide the samples into training set (70%) and test set (30%).

As the same as the previous experiment, we delete the 3rd index by probability α that vary from 0.1 to 0.9 by the step size 0.2. To delete the indices, we randomly sample user ids by probability α, and then delete the 3rd index (i.e., location) of the samples. For each α, we perform 5 trials and report the mean and the standard deviation of the results. Parameters are set to λ = 1.0 and R = 10.

We use the area under ROC curve (AUC) [Fawcett, 2006] as the evaluation metric. The AUC is often used to evaluate binary classification. In this experiment, our objective is to predict the value of the target data tensor, which takes 0 or 1. Hence this task can be evaluated as binary classification. To calculate the AUC, we first sort the test samples in descending order in terms of its predicted values. Then we calculate the number of true positives and false positives at each position in the ranking, and calculate the AUC using those numbers.

Results. Fig. 2(c) shows the results. First of all, our three proposed methods (i.e., MAP-EM, Prior, and Uniform) outperform CMTF and Minimal, which suggests that the proposed model works well even on the real dataset. However, when there are too many missing indices (α ≥ 0.7), Complete shows the highest AUCs, which agrees with the previous experiment on the small synthetic data.

Comparing the three proposed methods, MAP-EM and Prior show almost the same AUCs, which are larger than that of Uniform especially when α ≤ 0.5. We observed that the q inferred by MAP-EM is very close to the prior q in this experiment. This observation suggests that, even though MAP-EM can infer the global trend, it cannot infer the correct q for each sample n when the missing ratio is large.

1 http://noisy-text.github.io/2016/geo-shared-task.html

Although we do not report the actual computational time, it almost follows the theoretical time complexity.

6 Conclusion

In this paper, we define and address a tensor decomposition problem where the indices can be missing. We propose a probabilistic generative model that handles the missing indices as a latent variables, and develop the parameter inference algorithm based on the variational MAP-EM algorithm. The experimental results show that our proposed model outperforms the baselines except when there are only too few samples and the ratio of missing indices is too large. There are options how to deal with the latent variables. Performing the full variational MAP-EM works the best when the enough samples are available, whereas when we are given only scarce samples, simply using the prior q achieves reasonable results and reduces the computational cost.

There are several open problems as follows. First, it is important to theoretically guarantee how many samples are needed to recover the missing indices correctly. Second, since any tensor decomposition models can be plugged-in to our proposed model, it is interesting to develop the concrete parameter inference algorithm for those models. Finally, we can also consider different models that are not assuming Gaussian. For example, for 0-1 or count data such as Twitter data, assuming Poisson or Bernoulli distribution for p(xn|in, jn, kn, Θ) would be promissing.

Acknowledgements

This research was partially supported by the program "Research and Development on Real World Big Data Integration and Analysis" of RIKEN, Japan, and by a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO).