--- title: "Sampling for Approximate Maximum Search in Factorized Tensor" authors: - "Zhi Lu" - "Yang Hu" - "Bing Zeng" doi: "10.24963/ijcai.2017/334" url: "https://ijcai.org/proceedings/2017/334" pdf: "https://ijcai.org/proceedings/2017/0334.pdf" venue: "Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence" year: 2017 pages: "2400-2406" track: "Main track" keywords: - "Multidisciplinary Topics and Applications: Personalization and User Modeling" - "Machine Learning: Data Mining" - "Combinatorial & Heuristic Search: Combinatorial search/optimisation" --- # Sampling for Approximate Maximum Search in Factorized Tensor Zhi Lu, Yang Hu, Bing Zeng *Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence*, pages 2400-2406. DOI: [10.24963/ijcai.2017/334](https://doi.org/10.24963/ijcai.2017/334) ## Abstract Factorization models have been extensively used for recovering the missing entries of a matrix or tensor. However, directly computing all of the entries using the learned factorization models is prohibitive when the size of the matrix/tensor is large. On the other hand, in many applications, such as collaborative filtering, we are only interested in a few entries that are the largest among them. In this work, we propose a sampling-based approach for finding the top entries of a tensor which is decomposed by the CANDECOMP/PARAFAC model. We develop an algorithm to sample the entries with probabilities proportional to their values. We further extend it to make the sampling proportional to the $k$-th power of the values, amplifying the focus on the top ones. We provide theoretical analysis of the sampling algorithm and evaluate its performance on several real-world data sets. Experimental results indicate that the proposed approach is orders of magnitude faster than exhaustive computing. When applied to the special case of searching in a matrix, it also requires fewer samples than the other state-of-the-art method. ## Full text ## 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. Figure 1: Graph representation of factor matrices and sampled stars. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000000_3d7ad53f0cead3ce80d30d0c3f3172eace46c9859f5baa9e47fc73d8e1cb2798.png) 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. Figure 2: A sampled k-compound star with k = 2. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000001_5b9bf3feff58ef3086739aa15865d156ce82333f0e4fbfc89f52751d31d090a9.png) 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 ``` Table 1: Statistics for each data set. The last two columns show the largest and 1000-th largest entry value in the tensors. | Data | L1 | L2 | L3 | Top-1 | Top-103 | |-----------|--------|--------|-------|---------|-----------| | ML-S | 456 | 1,973 | 1,222 | 81.564 | 43.385 | | ML-L | 993 | 3,298 | 2,555 | 88.737 | 50.595 | | LastFM | 1,348 | 6,927 | 2,132 | 234.730 | 101.795 | | Delicious | 1,681 | 29,540 | 7,251 | 95.472 | 40.828 | | dblp | 43,928 | 2,572 | 17 | 4.190 | 2.452 | ## 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/ Figure 3: Recall of Corek for k ∈ {1, 2, 3}. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000002_eeaf298eded34bb9df34305164a47c27f3e02e696213a043714aae2b2c11d30c.png) ## 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. Figure 5: Running time in details for LastFM. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000003_9d5562632cf8759cc2a3b7d7776508da3d3646951013e31cb4f1f664212a4f96.png) 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 Figure 6: Performance of using prefiltering. The left bars in (b) and (c) are for t′ = s and the right are for t′ = s/10. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000004_ade3ef5512744c7a15d4c2e2100a2b70af6c2687752deab3943aafe230bfd265.png) Figure 7: Recall and running time with multiple users on LastFM. For better legibility, we only plot results for the first 100 users. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000005_91515faeb1258a3ec59030a2dddc903935f4f1ecdd1ea837bacd1ba4bcfe70d3.png) 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) Figure 8: Comparison with Diamond sampling on ML-L. ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000006_1893c4fea7d21b36236290deb7106c5de7237b4383a472123e8f6c6100689bb4.png) Figure 9: Plot of pairs (‖w‖1ˆxi/s, x3i) on ML-S. The solid line is the reference for equality. Two dashed lines are for x3i ± δx3i . ![Image](https://ijcai.org/proceedings/2017/fulltext/334/image_000007_da8971348a43c6617a24b0afc682a79b9947fa12a904cbfd46ae2a238a5ca8ce.png) 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. ## References 1. [Acar et al., 2009] Evrim Acar, Daniel M Dunlavy, and Tamara G. Kolda. Link prediction on evolving data using matrix and tensor factorizations. In the IEEE Internation- al Conference on Data Mining Workshops, pages 262-269, 2009. 2. [Adomavicius and Tuzhilin, 2015] Gediminas Adomavicius and Alexander Tuzhilin. Context-aware recommender sys- tems, pages 191-226. Springer, 2015. 3. [Andoni and Indyk, 2008] Alexandr Andoni and Piotr Indyk. Near-optimal hashing algorithms for approximate near- est neighbor in high dimensions. Communications of the ACM, pages 117-122, 2008. 4. [Bader and Kolda, 2007] Brett W. Bader and Tamara G. Kol- da. Efficient matlab computations with sparse and factored tensors. SIAM Journal on Scientific Computing, pages 205-231, 2007. 5. [Bader et al., 2012] Brett W. Bader, Tamara G. Kolda, et al. Matlab tensor toolbox version 2.5. Available online, 2012. 6. [Ballard et al., 2015] Grey Ballard, Tamara G. Kolda, Al- i Pinar, and C Seshadhri. Diamond sampling for approx- imate maximum all-pairs dot-product (MAD) search. In the IEEE International Conference on Data Mining, pages 11-20, 2015. 7. [Bertalmio et al., 2000] Marcelo Bertalmio, Guillermo S- apiro, Vincent Caselles, and Coloma Ballester. Image in- painting. In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, pages 417-424, 2000. 8. [Cantador et al., 2011] Ivan Cantador, Peter Brusilovsky, and Tsvi Kuflik. Second workshop on information het- erogeneity and fusion in recommender systems (HetRec 2011). In RecSys, pages 387-388, 2011. 9. [Coh, 1997] Approximating Matrix Multiplication for Pat- tern Recognition Tasks, 1997. 10. [Drineas et al., 2006] Petros Drineas, Ravi Kannan, and Michael W Mahoney. Fast monte carlo algorithms for ma- trices I: Approximating matrix multiplication. SIAM Jour- nal on Computing, 36(1):132-157, 2006. 11. [Dunlavy et al., 2011] Daniel M Dunlavy, Tamara G. Kolda, and Evrim Acar. Temporal link prediction using matrix and tensor factorizations. ACM Transactions on Knowl- edge Discovery from Data, pages 10:1-10:27, 2011. 12. [Harper and Konstan, 2016] F Maxwell Harper and Joseph A Konstan. The movielens datasets: History and context. ACM Transactions on Interactive Intelligent Systems, pages 19:1-19:19, 2016. 13. [Higham and Relton, 2016] Nicholas J Higham and Samuel D Relton. Estimating the largest elements of a matrix. SIAM Journal on Scientific Computing, pages C584-C601, 2016. 14. [Holodnak and Ipsen, 2015] John T Holodnak and Ilse CF Ipsen. Randomized approximation of the gram matrix: Exact computation and probabilistic bounds. SIAM Jour- nal on Matrix Analysis and Applications, pages 110-137, 2015. 15. [Hu et al., 2015] Yang Hu, Xi Yi, and Larry S Davis. Col- laborative fashion recommendation: A functional tensor factorization approach. In Proceedings of the 23rd ACM International Conference on Multimedia, pages 129-138, 2015. 16. [Jäschke et al., 2008] Robert Jäschke, Leandro Marinho, Andreas Hotho, Lars Schmidt-Thieme, and Gerd Stumme. Tag recommendations in social bookmarking systems. AI Communications, pages 231-247, 2008. 17. [Ke et al., 2004] Yan Ke, Rahul Sukthankar, and Larry Hus- ton. Efficient near-duplicate detection and sub-image re- trieval. In Proceedings of the 12th Annual ACM Interna- tional Conference on Multimedia, pages 869-876, 2004. 18. [Kolda and Bader, 2009] Tamara G. Kolda and Brett W. Bad- er. Tensor decompositions and applications. SIAM Review, pages 455-500, 2009. 19. [Koren et al., 2009] Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recom- mender systems. Computer, 42(8):30-37, 2009. 20. [Rendle et al., 2009a] Steffen Rendle, Leandro Balby Mar- inho, Alexandros Nanopoulos, and Lars Schmidt-Thieme. Learning optimal ranking with tensor factorization for tag recommendation. In Proceedings of the 15th Internation- al Conference on Knowledge Discovery and Data mining, pages 727-736. ACM, 2009. 21. [Rendle et al., 2009b] Steffen Rendle, Christoph Freuden- thaler, Zeno Gantner, and Lars Schmidt-Thieme. BPR: Bayesian personalized ranking from implicit feedback. In Proceedings of the 25th Conference on Uncertainty in Ar- tificial Intelligence, pages 452-461, 2009. 22. [Salton and Harman, 2003] Gerard Salton and Donna Har- man. Information retrieval. John Wiley and Sons Ltd., 2003. 23. [Shrivastava and Li, 2014] Anshumali Shrivastava and Ping Li. Asymmetric LSH (ALSH) for sublinear time maxi- mum inner product search (MIPS). In Advances in Neural Information Processing Systems, pages 2321-2329, 2014. 24. [Xiong et al., 2010] Liang Xiong, Xi Chen, Tzu-Kuo Huang, Jeff Schneider, and Jaime G Carbonell. Temporal collabo- rative filtering with bayesian probabilistic tensor factoriza- tion. In Proceedings of the SIAM International Conference on Data Mining, pages 211-222, 2010.