Sampling for Approximate Maximum Search in Factorized Tensor

1 Introduction

Matrix or tensor completion has received considerable attention in recent years. Many problems in application can be formulated as recovering the missing entries of a matrix or tensor. In some scenarios, such as image in-painting [Bertalmio et al., 2000], all lost entries are needed to be filled in. In some others, however, it is difficult and unnecessary to recover them all. Take recommender systems for example, the size of the matrix/tensor, determined by the numbers of users and items, is usually rather large. It is expensive computationally to compute all of the unknown values. On the other hand, for the recommendation purpose, we are only interested in a few entries that are the largest within a sub-array of the matrix/tensor. These entries are not only the central concerns for a personalized recommendation, but also meaningful in many other cases. In a similarity matrix, the top entries correspond to pairs of items that are most similar, which are of interest for applications like link prediction in graph [Dunlavy et al., 2011], duplicate detection [Ke et al., 2004] as well as information retrieval [Salton and Harman, 2003].

∗The corresponding author is Yang Hu. This work was supported by the National Natural Science Foundation of China (61602090) and the 111 Project (B17008).

In this work, we study the problem of efficiently identifying the top entries of a factorized tensor. Tensors, as multiway generalizations of matrices, hav been exploited more and more recently. Take recommender systems for example, while the traditional focus is on the user-item matrix [Koren et al., 2009], tensor is required for data representation in many emerging settings such as context-aware recommendation [Xiong et al., 2010; Adomavicius and Tuzhilin, 2015], and set-based recommendation [Hu et al., 2015] where the object to be recommended is a set of items such as a set of clothes that make an outfit.

Factorization models are essential tools for analyzing tensors. They assume that the tensor can be decomposed into some factors, which can be estimated from the observed entries by learning algorithms. Specifically, we focus on the CANDECOMP/PARAFAC [Kolda and Bader, 2009] decomposition model, a widely used technique for exploring and extracting the underlying structure of multi-way data. Given an N-order tensor X ∈ RL1×···×LN , the CP decomposition approximates it by N factor matrices:

where each factor matrix A(n) = [a(n) ∗1 a(n) ∗2 · · · a(n) ∗R ], n = 1, . . . ,N, is of size Ln ×R with a(n) ∗r ∈ RLn being the r-th column. The symbol "◦" represents the vector outer product. R is the tensor rank, indicating the number of latent factors. Element-wise, Eq.(1) is written as:

where i is short for the index tuple (i1, i2, . . . , iN). Given the factor matrices A(n) and a parameter t, we would like to find t index tuples {i1, . . . , it} that correspond to the top-t largest xi.

spond to nodes on the bottom, and the R columns of A(n) , indexed by r, correspond to nodes on the top. Edge (in, r) exists if a(n) inr is nonzero. Since all of the factor matrices have the same number of columns, they can share nodes for the columns. Therefore, the N factor matrices are represented by an (N+1)-partite graph with one partition, called the core partition, consisting of nodes for the common columns, and N peripheral partitions, each of which consists of nodes for the rows of a matrix.

To motivate our sampling procedure, we start with the case where all A(n) are binary matrices. If a core node r in core partition has a neighbor in, i.e. a(n) inr = 1, in each peripheral partition, we call subgraph (r, i1, . . . , iN) (shorted as (r, i)), a "star" (Figure 1 (b)). According to Eq.(2), xi is simply the number of distinct stars connecting the same group of neighbors in the peripheral partitions indexed by i. It can be shown that the probability of selecting a random star with endpoints i is proportional to xi [Ballard et al., 2015].

We further consider a "compound star", which is formed by two intersecting stars (r, i) and (r′ , i) (Figure 1 (c)). There are x2i compound stars with endpoints i. The probability of selecting a random compound star with endpoints i is proportional to x2i . More generally, if a k-compound star is formed by k intersecting stars, its selected probability will be proportional to xki .

Sampling random k-compound stars will accelerate identifying the top entries as compared to sampling basic stars but with higher complexity. Besides, the factor matrices A(n) are usually real-valued, which requires the sampling to have a nonuniform distribution. In the following, we present a sampling algorithm that handles these issues carefully.

3.1 Sampling A Star

To sample k-compound stars, we take node set that contains k nodes from the core partition as a k-compound core node.

matrices.

2 RelatedWork

Our problem subsumes many existing problems in the literature. When N = 2, it is exactly the MAD (Maximum Allpairs Dot-product Search) problem [Ballard et al., 2015] that finds the largest entries in the product of two matrices. Notice that MAD contains the MIPS (Maximum Inner Product Search) [Coh, 1997] problem as the special case with one matrix being a single column.

The most obvious solution to the problem is to compute the entries exhaustively. However, this becomes prohibitive as the sizes of factor matrices grow. There is some literature in approximate matrix multiplication. However, these methods are not suited even for MAD, since only a few entries among millions are of interest. The more efficient solution is to directly search the top ones, which has been extensively studied for the MIPS problem. Popular approaches include LSH (Locality Sensing Hashing) [Andoni and Indyk, 2008; Shrivastava and Li, 2014], space partition techniques like kd tree, and sampling-based approaches [Drineas et al., 2006; Holodnak and Ipsen, 2015]. Recently, Ballard et al. proposed a randomized approach called diamond sampling [Ballard et al., 2015] to the MAD problem. [Higham and Relton, 2016] derived an algorithms for estimating the largest elements of a matrix using a few matrix-vector products with it and its conjugate transpose. For tensor, however, there has not been any study conducted yet.

Inspired by the work of diamond sampling, we present an index sampling method for a factorized tensor, whose entries are computed by the CP decomposition model. We design a strategy to sample the index tuple i proportional to the k-th power of the entries, which amplifies the focus on the largest ones. We provide theoretical analysis for our sampling algorithm and derive concentration bound on its behavior. For applications in recommender systems, an algorithm that reuses the samples for multiple users is presented. We evaluate the proposed method on several real-world data sets. The results show that it is orders of magnitude faster than exact computing. When compared to previous approach for matrix sampling, our method requires much fewer samples.

3 Star Sampling

A graph is used to represent the factor matrices of CP decomposition. Each matrix A(n) ∈ RLn×R, n ∈ [N], is represented by a bipartite graph as shown in the dashed box in Figure 1 (a). The Ln rows of A(n) , indexed by in, corre-

Algorithm 1 The Corek Sampling Method

Given factor matrix A(n) ∈ RLn×R, n = 1, 2, . . . , N. Let s be the number of samples. 1: for r ∈ [R] × · · · × [R] ︸ ︷︷ ︸ k do 2: for n = 1, ..., N do 3: for in = 1, . . . , Ln do 4: e(n) inr ← a(n) inr1 · · · a(n) inrk 5: end for 6: ‖e(n) ∗r ‖1 ← ∑in |e(n) inr | 7: end for 8: wr ← ‖e(1) ∗r ‖1 · · · ‖e(N) ∗r ‖1 9: end for 10: for ℓ = 1, . . . , s do 11: Sample r with probability wr/‖w‖1 12: for n = 1, ..., N do 13: Sample in with probability |e(n) inr |/‖e(n) ∗r ‖1 14: end for 15: ξi,ℓ ← sign(e(1) i1r · · · e(N) iN r ) 16: if i = (i1, i2, . . . , iN) has not been sampled then 17: Create ˆxi ← ξi,ℓ 18: else 19: ˆxi ← ˆxi + ξi,ℓ 20: end if 21: end for

There are Rk k-compound nodes in total. Figure 2 shows a sampled star with two core nodes (r and r′ ) coalesced to a single compound core node (r, r′ ). We weigh the edge between a k-compound core node r = (r1, . . . , rk) and a node in in the n-th peripheral partition by

Let e(n) ∗r = [e(n) 1r , · · · , e(n) Lnr ]T be a vector containing weights from r to all nodes in the n-th peripheral partition, we assign weight for the k-compound node r as

where ‖ · ‖1 denotes 1-norm of a vector.

To find a k-compound star, we first select a k-compound node with probability wr/‖w‖1, where w is the vector containing all wr. Then, conditioned on the sampled compound node, we pick one node from each of the peripheral partitions. The node in the n-th peripheral partition is drawn with probability |e(n) inr |/‖e(n) ∗r ‖1. After an index tuple i = (i1, i2, . . . , iN) is obtained, a score ξi,ℓ is computed:

where ℓ denotes the ℓ-th sample. If the index tuple i has not been sampled before, a container is created with ˆxi = ξi,ℓ. Otherwise, we increase ˆxi by ξi,ℓ. The full procedure is shown in Algorithm 1.

3.2 Top-t Extraction

Let Ψp = {ip|p ∈ [P]} denotes the set of index tuples that have been sampled. P ≤ s since some i may be picked more than once. Two strategies for identifying the t largest entries from Ψp are presented. The first one is directly computing the exact entry value xi for each sample in Ψp using Eq.(2) and then extracting the t largest ones. Following the previous works [Ballard et al., 2015; Coh, 1997], in the second strategy, we utilize the scores ˆxi obtained during sampling to do prefiltering since the load for computing ˆxi is much lighter than that of computing xi. Moreover, as we show through theoretical analysis later, the expectation of ˆxi is proportional to the k-th power of xi. Therefore, ˆxi can be used to filter out the relatively small ones in Ψp. Specifically, we denote a subset by Ψt′ that contains t′ index tuples corresponding to the top-t′ largest ˆxi after the prefiltering. Then, the exact values are computed only for index tuples in Ψt′ , and the t index tuples with the largest xi in Ψt′ are taken as the output of our algorithm.

In the following part, we prove that in Algorithm 1, the probability that an index tuple i is sampled is proportional to xki and the expectation of ˆxi is also proportional to the k-th power of xi.

3.3 Theoretical Analysis

Let ϵi,r denote the event of choosing compound node r and peripheral nodes i in one iteration, while ϵi be that of choosing i. We first analyze p(ϵi), the probability that an entry is sampled, and then compute the expectation of ˆxi. We also provide error bound of the estimate.

Lemma 1. Suppose that all elements of the factor matrices A(n) , n ∈ [N], are nonnegative. We have p(ϵi) = xki /‖w‖1.

Proof. The probability p(ϵi) is marginal distribution of p(ϵi,r), which can be computed by:

Let ci,ℓ be a random variable that if i has been picked in the ℓ-th iteration, ci,ℓ = ξi,ℓ. Otherwise ci,ℓ = 0. The final score ˆxi can be written as ˆxi = ∑sℓ=1 ci,ℓ.

Lemma 2. The expectation of ˆxi equals to sxki /‖w‖1.

Proof. Since ci,ℓ are i.i.d for fixed i and varying ℓ, the expectation of ˆxi is:

Theorem 1. Fix δ > 0 and error probability σ ∈ (0, 1). Assume that all entries in the factor matrices are nonnegative. If the number of samples

Proof. Since ci,1, · · · , ci,s are independent random variables taking values in {0, 1}. The Chernoff bounds on the sum of Poisson trials shows Pr(|ˆxi - µ| ≥ δµ) ≤ 2 exp (-µδ2 /3), δ ∈ (0, 1), where µ = sxki /‖w‖1. By the choice of s, we have µ ≤ 3 log (2/σ)/δ2 . Then, Pr(|ˆxi - sxki /‖w‖1| ≥ δsxki /‖w‖1) ≤ σ. Multiplying ‖w‖1/s inside Pr[·] proofs the bound.

4 Sampling for Recommender Systems

In this section, we present an algorithm for efficiently identifying top entries in sub-arrays of a tensor, which is of particular interest for recommender systems. Without loss of generality, we assume that the first factor matrix A(1) corresponds to users, with each row characterizing a single user. For any user u = a(1) i1∗ , we are interested in top entries in the sub-tensor Xu, whose CP decomposition is

We find that for different users, only the probabilities for sampling k-compound nodes from the core partition are different. In the following step, they share the same distribution for sampling the peripheral nodes. Based on this observation, for each k-compound node in the core partition, we build a pool containing sub-index tuples i′ = (i2, . . . , iN) that have been picked given that node. For a new user, when a k-compound node is chosen, we can directly use the indices kept in its pool and only sample new i′ when necessary. By sharing among users, a huge number of samples can be saved. This procedure is illustrated in Algorithm 2.

Algorithm 2 Finding top-t entries for multiple users

Given factor matrix A(n) ∈ RLn×R, n = 1, 2, . . . , N. A(1) contains vectors characterizing users. Let s be the number of samples, and m = Rk .

1: Initialize m empty pools g1, g2, . . . , gm 2: Initialize fr = 0 for all compound nodes r. 3: for i1 = 1, 2, . . . ,L1 do 4: Let the current user be u = a(1) i1∗ 5: for all r = (r1, . . . , rk) do 6: wr ← |ur1 · · · urk |‖e(2) ∗r ‖1 · · · ‖e(N) ∗r ‖1 7: end for 8: for ℓ = 1, . . . , s do 9: Sample r with probability wr/‖w‖1. 10: Increase fr by 1. 11: end for 12: for all r do 13: if fr > |gr| then 14: Sample fr - |gr| index tuples i′ and append them to gr. 15: end if 16: Use fr index tuples i′ in gr and compute scores. 17: end for 18: Post-processing for finding top-t entries of u. 19: end for

5 Implementation Details

In this section, we discuss some details that could improve the performance. Follow the previous work [Coh, 1997], instead of first picking a node from core partition and then choosing N peripheral nodes in each iteration, we directly get s core nodes in one loop. Specifically, we compute the expected number of occurrence for each compound node. Let µr = swr/‖w‖1, the count for r is:

This requires O(Rk ) work while sampling r s times (line 10 of Algorithm 1) requires O(sk log R) work.

To conduct Corek sampling, we may extend the factor matrices to E(n) whose elements are from Eq.(3). Then, we can just conduct Core1 sampling with E(n) which requires O(LnRk ) space for each n ∈ [N]. This is costly in practice when k is large. On the other hand, given a core node r, the sampling of the peripheral nodes will only depend on one column in E(n) indexed by r. Therefore, instead of storing the whole E(n) , we only keep one column of it at a time. This requires O(L1 + · · · + LN) space in total. When k is small, we can directly store E(n) in exchange for speed.

6 Experiments

We evaluate our algorithms on five real-world data sets: Delicious Bookmarks 1 (Delicious), Last.FM 2 (LastFM), MovieLens 3 +IMDb 4 /Rotten Tomatoes 5 (ML-S), a larger MovieLens data set (ML-L) and dblp 6 . The first four are from the tag recommendation application [Cantador et al., 2011; Harper and Konstan, 2016]. Following previous practice [J¨aschke et al., 2008], users, items and tags that occur at least 5 times are used. The last is for link prediction. We follow [Acar et al., 2009] to extract only the proceedings from 1991 to 2007.

To learn the factor matrices of CP decomposition, for the tag data sets, we use the Bayesian Probabilistic Tensor Factorization (BPTF) method and optimize the pair-wise ranking function of Bayesian Personalized Ranking (BPR) [Rendle et al., 2009b; 2009a]. And for the dblp data set, Alternating Least Squares (ALS), provided in the Tensor Toolbox of [Bader et al., 2012; Bader and Kolda, 2007], is used. R is chosen following previous works, i.e. 64 for the first four data sets and 50 for the last one. We give the summary statistics in Table 1.

1 http://www.delicious.com

2 http://www.lastfm.com

3 http://www.grouplens.org

4 http://www.imdb.com

5 http://www.rottentomatoes.com

6 http://dblp.uni-trier.de/db/

6.1 Accuracy and Running Time

We experiment with k equals 1,2 and 3 respectively. In Figure 3, we plot recall, i.e. the percentage of true top-t entries identified, versus the number of samples s. The corresponding running time are illustrated in Figure 4. The time for exhaustive computing is shown as the baseline. Here we set t′ = s and focus on the effectiveness of the sampling stage firstly. For each s, we run 10 times and the average result is reported.

We can see that with the increase of k, much higher recall can be obtained with the same number of samples. For the tag recommendation data sets, starting from s = 105 , the time consumed by Core2 becomes similar with Core1 . For dblp, they cost almost the same time from s = 106 . Core3 consumes much more time than the other two when the number of samples are small due to its relatively heavy initialization. However, when the number of samples reaches 107 , its running time become quite similar with the other two. All three sampling methods are orders of magnitude faster than exhaustive computing.

To show more details on running time, we separate it into four stages: initialization, sampling, scoring and filtering. Initialization consists of computation of the weights for core nodes and some preprocessing. Sampling corresponds to sampling index tuples and may also include the computation of the distributions for peripheral partitions when extension matrices are not pre-saved. The scoring phase shows the time consumed on computing and saving scores for each sampled index tuple. The last is to identify the top-t entries through filtering. We plot them for LastFM in Figure 5.

The running time for initialization is dependent on k and will be less influential with the increase of s. When s is relatively small (Figure 5 (a)), the sampling phase in Core3 is dominated by its online computing of the probability vectors, which leads to considerable more running time when compared with Core1 and Core2 . And when s is large (Figure 5 (b)), all sampling phases are dominated by the sampling procedure, and are therefore similar in run time. When k = 3, the probability for top entries are amplified, leading to a smaller size of Ψp and a saving of computation during the filtering phase. We notice that with the increase of s, the running time of Core3 becomes quite similar with the other two.

6.2 Use Prefiltering

In this part, we evaluate the performance of using prefiltering during top-t extraction (Section 3.2). Due to space limit, we only show results on LastFM with t′ = s/10 in Figure 6. We notice that prefiltering leads to a slightly lower recall but save considerable computing time.

6.3 Sampling for Collaborative Filtering

In Figure 7, we show the performance of Algorithm 2 applied for collaborative filtering, where the samples are shared among different users. Here Core2 is used to find the top-100

largest entries on LastFM data set and s is set to 106 . We can see that only for the first few users we need some time to get the samples and build the sample pools. After that, the computation time for a new user is much reduced.

6.4 Comparison with Diamond Sampling

Since this is the first study of efficiently identifying the top entries in factorized tensors, there is no previous approach that we can compare with. On the other hand, when N = 2, the problem reduces to the MAD problem which finds the top entries in the product of two matrices. We therefore apply our algorithms to this problem and compare their performance with the diamond sampling method of [Ballard et al., 2015].

For each user, we can multiply its characterizing vector into the item matrix to get a new user-dependent item matrix. Note that the searchings for different users are conducted independently, i.e. we run Algorithm 1 instead of Algorithm 2. In Figure 8, we show the results on ML-L. We set t′ = s/10 and evaluate recall for t ∈ {1, 10, 100, 1000}. Since there are L1 pairs of matrices, the average results are drawn.

The performance of Core1 is slightly worse than diamond sampling since the score ˆxi in the latter is proportional to the square of entry magnitude. With the same number of samples, Core2 can achieve higher recall than diamond sampling with less computing time. With Core3 , the recall can be further improved with the cost of some more time. This illustrates the advantage of augmenting the sampling probability to higher power of the entry magnitude.

6.5 Accuracy of the Scores

To show the accuracy of the estimation ˆxi, we run Core3 on ML-S with s = 108 and plot the 104 pairs of (x3i , ‖w‖1ˆxi/s)

which have the largest x3i in Ψp after sampling in Figure 9. Given s = 108 , we set σ = 0.1 and compute the δ for each xki using the equality in Eq.(6) from Theorem 1. Then two dashed lines for xki ± δxki can be plotted As expected, the points concentrate around the diagonal and most pairs fall within the two dashed lines, which confirmed with the theoretical result.

7 Conclusion

In this work, we proposed a novel sampling-based algorithm for approximate maximum search in tensor which is factorized by CP decomposition. It is an extension of MAD and this is the first study of this important problem. Experimental results indicate that our method is orders of magnitude faster than exhaustive search. When applied to the MAD problem, it is also more efficient than other state-of-the-art method.