TL;DR
This paper introduces a fast randomized algorithm for matrix and tensor interpolative decomposition using CountSketch, offering significant speed improvements while maintaining accuracy, applicable to large-scale data.
Contribution
The paper presents a novel CountSketch-based randomized algorithm for matrix and tensor interpolative decomposition with theoretical guarantees and improved computational efficiency.
Findings
Achieves at least an order of magnitude speed-up on large matrices and tensors.
Maintains accuracy comparable to existing methods.
Provides theoretical performance guarantees for both matrix and tensor cases.
Abstract
We propose a new fast randomized algorithm for interpolative decomposition of matrices which utilizes CountSketch. We then extend this approach to the tensor interpolative decomposition problem introduced by Biagioni et al. (J. Comput. Phys. 281, pp. 116-134, 2015). Theoretical performance guarantees are provided for both the matrix and tensor settings. Numerical experiments on both synthetic and real data demonstrate that our algorithms maintain the accuracy of competing methods, while running in less time, achieving at least an order of magnitude speed-up on large matrices and tensors.
Click any figure to enlarge with its caption.
Figure 1
Figure 2| Algorithm for tensor ID | Complexity |
|---|---|
| Gram matrix | |
| Gaussian | |
| CountSketch (Proposal, Alg. 3) |
| Algorithm for matrix ID | Error | Run time (s) |
|---|---|---|
| Gaussian | 1.505e−15 | 20.38 |
| SRFT | 1.507e−15 | 18.40 |
| CountSketch (proposal) | 1.504e−15 | 0.59 |
| Method B | |||
|---|---|---|---|
| Method A | Our Proposal | SRFT | Gaussian |
| Our Proposal | - | 57 | 34 |
| SRFT | 3 | - | 2 |
| Gaussian | 26 | 58 | - |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
∎
11institutetext: O. A. Malik 22institutetext: Department of Applied Mathematics, University of Colorado Boulder
22email: [email protected]
33institutetext: S. Becker 44institutetext: Department of Applied Mathematics, University of Colorado Boulder
44email: [email protected]
Fast Randomized Matrix and Tensor Interpolative Decomposition Using CountSketch††thanks: This version of the article has been accepted for publication, after peer review (when applicable) and is subject to Springer Nature’s AM terms of use, but is not the Version of Record and does not reflect post-acceptance improvements, or any corrections. The Version of Record is available online at: http://dx.doi.org/10.1007/s10444-020-09816-9
Osman Asif Malik
Stephen Becker
(Received: date / Accepted: date)
Abstract
We propose a new fast randomized algorithm for interpolative decomposition of matrices which utilizes CountSketch. We then extend this approach to the tensor interpolative decomposition problem introduced by Biagioni et al. (J. Comput. Phys. 281, pp. 116–134, 2015). Theoretical performance guarantees are provided for both the matrix and tensor settings. Numerical experiments on both synthetic and real data demonstrate that our algorithms maintain the accuracy of competing methods, while running in less time, achieving at least an order of magnitude speed-up on large matrices and tensors.
Keywords:
Matrix Decomposition Tensor Decomposition Sketching
MSC:
15-02
††journal: Advances in Computational Mathematics
1 Introduction
Matrix decomposition is a fundamental tool used to compress and analyze data, and to improve the speed of computations. For data and computational problems involving more than two dimensions, analogous tools in the form of tensors and associated decompositions have been developed (Kolda and Bader, 2009). In many modern applications, matrices and tensors can be very large, which makes decomposing them especially challenging. One approach to dealing with this problems is to incorporate randomization in decomposition algorithms (Halko et al., 2011). In this paper, we consider the interpolative decomposition (ID) for matrices, as well as the tensor ID problem. By tensor ID, we mean the tensor rank reduction problem as introduced by Biagioni et al. (2015); we provide an exact definition in Section 1.2.2. We make the following contributions in this paper:
- •
We propose a new fast randomized algorithm for matrix ID and provide theoretical performance guarantees.
- •
We propose a new randomized algorithm for tensor ID. To the best of our knowledge, we provide the first performance guarantees for any randomized tensor ID algorithm.
- •
We validate our algorithms on both synthetic and real data.
- •
We propose a small modification to the standard CountSketch formulation which helps avoid certain rank deficiency issues and slightly strengthen our matrix ID results.
1.1 Tensors and the CP Decomposition
For a more complete introduction to tensors and their decompositions, see the review paper by Kolda and Bader (2009). A tensor is an -dimensional array of real numbers, also called an -way tensor. The number of elements in such a tensor is denoted by . Boldface Euler script letters, e.g. , denote tensors of dimension 3 or greater; bold capital letters, e.g. , denote matrices; bold lowercase letters, e.g. , denote vectors; and lowercase letters, e.g. , denote scalars. Uppercase letters, e.g. , are used to denote scalars indicating dimension size. A colon is used to denote all elements along a certain dimension. For example, and are the th row and th column of the matrix , respectively. If is a vector of column indices, then denotes the submatrix of consisting of the columns of whose indices are listed in . denotes the identity matrix. For a matrix , denotes its th singular value, and and denote the maximum and minimum singular values, respectively. The condition number of is defined as . The number of nonzero elements of is denoted by . For positive integers and , let and . The Hadamard product, or element-wise product, of matrices is denoted by . The Khatri–Rao product of matrices is denoted by . The tensor Frobenius norm is denoted by , where flattens the tensor into a column vector. A norm with no subscript will always denote the matrix spectral norm.
The singular value decomposition (SVD) decomposes matrices into a sum of rank-1 matrices (Golub and Van Loan, 2013). Similarly, the CP decomposition decomposes a tensor into a sum of rank-1 tensors:
[TABLE]
where denotes outer product, and each is a rank-1 tensor. Each is called an s-value, each is called a factor matrix, and all vectors have unit 2-norm. Usually, a tensor is said to be of rank- if is the smallest possible number of terms required in a representation of the form (1). We will use the term “rank” in a looser sense to mean the (not necessarily minimal) number of rank-1 terms in a representation of the form (1).
1.2 Interpolative Decomposition
1.2.1 Matrix Interpolative Decomposition
For a matrix , a rank- interpolative decomposition (ID) takes the form , where consists of a subset of columns from , and is a coefficient matrix which is well-conditioned in some sense. The fact that the decomposition is expressed in terms of the columns of means that inherits properties such as sparsity and non-negativity from . Moreover, expressing the decomposition in terms of columns of can increase interpretability. Algorithm 1 outlines one method to compute a matrix ID.
Fact 1
If the partial QR factorization on line 3 in Algorithm 1 is done using the strongly rank-revealing QR (SRRQR) decomposition developed by Gu and Eisenstat (1996), then Algorithm 1 has complexity (Cheng et al., 2005). Moreover, the decomposition it produces satisfies the following properties (Martinsson et al., 2011):
- (i)
Some subset of the columns of makes up the identity matrix, 2. (ii)
no entry of has an absolute value exceeding 2, 3. (iii)
, 4. (iv)
, 5. (v)
when or , and 6. (vi)
when .
In practice, using a variant of column pivoted QR instead of the SRRQR on line 3 of Algorithm 1 works just as well, and reduces the complexity of the algorithm to (Cheng et al., 2005).
There have been subsequent proposals for randomized versions of matrix ID (Liberty et al., 2007). Martinsson et al. (2011) propose a variant which incorporates Gaussian random sketching. It computes a sketch , where () is a matrix with iid standard normal entries, and then computes an ID . The same and then give an ID of . Woolfe et al. (2008) propose a similar fast randomized algorithm which uses a subsampled randomized fast Fourier transform (SRFT) instead of a Gaussian matrix. It computes a sketch , where is a diagonal matrix with each diagonal entry iid and equal to or with equal probability, is the fast Fourier transform (FFT), and is a subsampling operator that randomly samples rows.
1.2.2 Tensor Interpolative Decomposition
Biagioni et al. (2015) consider the problem of rank reduction of a CP tensor, which they call tensor ID. Suppose is an -way tensor with CP decomposition (1). Computing a rank-, , tensor ID of amounts to finding a representation
[TABLE]
where contains unique indices. Tensor ID has many applications. For example, in various algorithms, the rank of discretized separated representations of multivariate functions grows with each iteration, requiring repeated rank reduction of CP tensors (Beylkin and Mohlenkamp, 2002, 2006). Another example is the algorithm by Reynolds et al. (2017) for finding the element of maximum magnitude in a CP tensor which also requires repeated rank reduction.
Biagioni et al. (2015) approach the tensor ID problem by considering the matrix
[TABLE]
where is a diagonal matrix with entries . The tensor ID problem can now be reduced to identifying columns of using matrix ID. However, when the factor matrices have no special structure, has elements and is therefore typically infeasible to form. One way to tackle this problem is by forming the much smaller Gram matrix , which can be done using flops since
[TABLE]
compute its symmetric matrix ID, and use it to compute an ID of . This approach, however, can lead to accuracy issues since the Gram matrix can be ill-conditioned, since (Biagioni et al., 2015). Biagioni et al. (2015) therefore propose a randomized method which avoids the ill-conditioning issue and reduces the complexity. This is done by applying a kind of Gaussian sketch to , but instead of forming a full Gaussian matrix of size , a matrix of the form
[TABLE]
is used, where each is a matrix with elements that are iid standard normal random variables. The sketch can then be computed efficiently without ever forming or , since . Note that the elements of in (5) are not independent. This means that the theory for Gaussian matrix ID, which requires independence, cannot be used to provide guarantees for sketched matrix ID using .
1.3 Basics of CountSketch
Our proposed method uses a type of sketching called CountSketch (Charikar et al., 2004; Clarkson and Woodruff, 2017), which we now describe. Let be a random map such that each is iid and , let be a matrix with and all other entries equal to 0, and let be a diagonal matrix with each diagonal entry iid and equal to or with equal probability. The CountSketch operator is then defined as . Applying to does the following: The matrix changes the sign of each row of with probability , and the matrix then randomly adds each row of to one of target rows. Due to the special structure of , it can be applied implicitly with complexity (Clarkson and Woodruff, 2017).
Suppose has the special structure , where each . For such matrices, there is a variant of CountSketch which allows computing the sketch of without ever having to form the full matrix, which can be prohibitively large to store explicitly. This variant is called TensorSketch and is developed by Pagh (2013), Pham and Pagh (2013), Avron et al. (2014) and Diao et al. (2018). It works as follows:
- •
Define independent random maps such that each is iid and ; and
- •
define independent random sign functions such that .
Next, define as
[TABLE]
and as
[TABLE]
Notice that each row index of corresponds to a unique -tuple . and can therefore be considered functions on . With this in mind, let denote a diagonal matrix with the th diagonal entry equal to . If and are used instead of and in the definition of CountSketch above, we get TensorSketch, which we will denote by . The reason for choosing this formulation is that it can be computed efficiently using the following formula:
[TABLE]
where each is a CountSketch operator defined using and the diagonal matrix . The formula (8) follows from the discussion in Section A in the supplementary material of Diao et al. (2018). Other good sources for further details on TensorSketch are Pagh (2013), Pham and Pagh (2013) and Avron et al. (2014).
2 Other Related Work
We provided an overview of existing ID algorithms in Section 1.2. The matrix ID is related to the CX and CUR decompositions (Drineas and Kannan, 2003; Drineas et al., 2006, 2008; Mahoney and Drineas, 2009; Bien et al., 2010; Wang and Zhang, 2013; Boutsidis and Woodruff, 2017), also known as skeleton approximations (Goreinov et al., 1997a, b; Tyrtyshnikov, 2000), and the column subset selection problem (Frieze et al., 2004; Deshpande and Vempala, 2006; Deshpande et al., 2006; Boutsidis et al., 2009; Deshpande and Rademacher, 2010; Guruswami and Sinop, 2012; Boutsidis et al., 2014). Like ID, the CX decomposition takes the form , where contains a subset of the columns of . The crucial feature that distinguishes ID from a CX decomposition is the additional conditioning requirements on the coefficient matrix in ID; the matrix in a CX decomposition is not required to have the properties (i)–(iv) listed in Fact 1 (Drineas et al., 2008). A CUR decomposition takes the form , where and contain a subset of the columns and rows of , respectively. Consequently, setting would yield a CX decomposition. It is well-known that the matrix defined in this manner is typically ill-conditioned (Voronin and Martinsson, 2017). Since we require the coefficient matrix in our decomposition to be well-conditioned, the available algorithms for CX and CUR decomposition are not useful to us.
Various randomized algorithms have been utilized in the context of tensor decomposition before. Examples include the works of Wang et al. (2015), Battaglino et al. (2018), and Yang et al. (2018) for the CP decomposition; Drineas and Mahoney (2007), Tsourakakis (2010), da Costa et al. (2016) and Malik and Becker (2018) for the Tucker decomposition; and Zhang et al. (2018) and Tarzanagh and Michailidis (2018) for t-product based decompositions. Other notable works that use CUR-type algorithms or sampling are e.g. those by Mahoney et al. (2008), Caiafa and Cichocki (2010), Oseledets et al. (2008) and Friedland et al. (2011). The tensor ID which we consider is different from the various problems solved in these previous papers. The goal of tensor ID is not to compute a tensor decomposition from an arbitrary data tensor. Instead, the purpose of tensor ID is to compress a tensor which is already in CP format in an efficient and principled manner. To the best of our knowledge, the only work aside from that by Biagioni et al. (2015) which considers randomized tensor ID is the paper by Reynolds et al. (2016). They introduce a randomized alternating least-squares (ALS) algorithm, which is better conditioned but slower than the standard ALS algorithm for CP decomposition. Biagioni et al. (2015) conclude that standard ALS is much slower than their Gaussian sketching algorithm. We therefore do not compare our proposed tensor ID to the randomized ALS by Reynolds et al. (2016) since it is even slower.
To the best of our knowledge, TensorSketch was the first sketch with theoretical guarantees that could be applied particularly efficiently to matrices like in (3) with Kronecker structured columns. Recently, a number of works have appeared that provide guarantees for other methods designed for efficient sketching of Kronecker structured vectors. Sun et al. (2018) consider sketches of the form (5) where each has sub-Gaussian entries. They provide Johnson–Lindenstrauss (JL) style guarantees for the case when the sketch is a Khatri–Rao product of two smaller matrices, i.e., for in (5). Rakhshan and Rabusseau (2020) consider sketches which have tensor train or CP tensor structure, and with core tensors and factor matrices that have Gaussian entries. They provide JL style guarantees for these sketches for arbitrary orders of the random tensors. Their sketch with CP tensor structure includes the sketch in (5) as a special case. Another line of work considers the Kronecker fast JL transform, which is a structured variant of the fast JL transform of Ailon and Chazelle (2009). It was first proposed by Battaglino et al. (2018) with theoretical guarantees later provided by Jin et al. (2019) and Malik and Becker (2020). A variant of this transform is also considered by Iwen et al. (2019).
3 Fast Randomized Matrix ID Using CountSketch
Algorithm 2 explains our proposal for CountSketch matrix ID. Proposition 1 provides guarantees for the method. A proof is provided in Section 7.1, which also contains a more detailed version of the bound in (10).
Proposition 1 (CountSketch matrix ID)
Suppose , and are defined as in Algorithm 2. Let be a real number and a positive integer such that
[TABLE]
Suppose that the matrix ID on line 5 of Algorithm 2 utilizes SRRQR. Then, the output of Algorithm 2 satisfies properties (i)–(iv) in Fact 1. Moreover, the outputs satisfy
[TABLE]
with probability at least .
The condition in (9) is very similar to that for SRFT matrix ID by Woolfe et al. (2008). The only difference is that instead of a term of the form , their work only has a factor . In practice, the condition in (9) is very conservative. We find that a small oversampling factor, e.g. , works well in practice, producing errors of the same size as the other randomized ID methods.
Remark 1
The semi-coherent matrices defined by Avron et al. (2010) are adversarial to CountSketch, and therefore to our proposed method. For such matrices, using may result in a large error for our method. Some care is therefore necessary when applying our method together this choice of . In Section 6.1.2, we do extensive testing of our method on real-world matrices to demonstrate that using works well in practice. We also provide an example of a matrix with semi-coherent structure on which our method fails when choosing like this.
Remark 2
In cases when the target rank is quite large (e.g. ), an issue we encountered is that can be rank deficient due to rank deficiency of . This issue can be dealt with easily by defining slightly differently to ensure that each row contains at least one nonzero element. This is done by the following straightforward modification of the map in the definition of CountSketch in Section 1.3: Let each , where is a uniform random permutation of the elements of the vector , where each is iid uniformly random. With this modification, the guarantees of Proposition 1 still hold. In fact, the condition in (9) is slightly improved. We give a precise statement with proof in Section 7.2.
4 Extending the Results to Tensor ID
Let and be defined as in (1) and (2), respectively. Our approach to the tensor ID problem is similar to that of Biagioni et al. (2015): We sketch the matrix in (3) without forming it and compute a matrix ID of this sketch. The approximation is then constructed using the rank-1 components of corresponding to the columns of used in the ID of that matrix. The s-values used in the representation of are then computed as , for . The sketch we use is the efficient TensorSketch variant of CountSketch. Algorithm 3 outlines our proposed method for tensor ID. Proposition 2 provides guarantees for the method. A proof is provided in Section 7.3. To the best of our knowledge, there are no previous results like Proposition 2 for randomized tensor ID.
Proposition 2 (TensorSketch tensor ID)
Suppose , and are defined as in Algorithm 3. Let be a real number and a positive integer such that . Suppose that the matrix ID on line 6 of Algorithm 3 utilizes SRRQR. Then, the output of Algorithm 3 satisfies
[TABLE]
with probability at least .
As mentioned in Section 1.2.2, an issue with forming and then decomposing is that it can be ill-conditioned. Biagioni et al. (2015) point out that the sketched matrix typically is much better conditioned since and Gaussian matrices are well-conditioned. As Proposition 3 demonstrates, the matrix is also well-conditioned with high probability, when the sketch dimension is sufficiently large. A proof of Proposition 3 is provided in Section 7.4.
Proposition 3
Let be a real number, and let and be positive integers such that . Suppose is a TensorSketch matrix, and is an arbitrary matrix. Then with probability at least .
5 Complexity Analysis
In this section, we compare the complexity of our proposed methods with the other algorithms. We assume all QR factorizations are done using column pivoted QR instead of SRRQR, and ignore the cost of generating random variables. We also assume that where is a small positive integer (e.g. ) since this choice works well in practice. Since , and we assume for a small constant , we also make the assumption .
The costs of the different steps of Algorithm 2 are as follows:
- •
Computing the sketch : .
- •
Computing Matrix ID of , where : .
The total cost is therefore . The cost of standard matrix ID can be found in Remark 3 of Cheng et al. (2005), and the cost of SRFT matrix ID can be found in Remark 5.4 of Woolfe et al. (2008). The cost of Gaussian matrix ID is straightforward to compute similarly to our computation above. Table 1 summarize these matrix ID complexities.
For the tensor ID algorithms, we assume the input is an -way rank- CP tensor of size , and that each factor matrix has the same number of nonzeros, which we denote by . The costs of the different steps of Algorithm 3 are as follows:
- •
Computing the TensorSketched matrix : .
- •
Computing Matrix ID of , where : .
- •
Computing : .
The total cost is therefore . Although Biagioni et al. (2015) do not specify these, the complexities for the Gram matrix approach and Gaussian tensor ID can be computed from the descriptions in their paper. Table 2 summarize the complexities for the different tensor ID algorithms. The constant is the cost of computing one Gram matrix in (4), which we assume is the same for each , e.g. if the factor matrices were dense.
6 Numerical Experiments
The numerical experiments are done in Matlab R2018b and C. All results are averages over ten runs in an environment using four cores of an Intel Xeon E5-2680 v3 @2.50GHz CPU and 19 GB of RAM. All code used to generate our results can be found at https://github.com/OsmanMalik/countsketch-matrix-tensor-id, including implementations of our proposed methods. For all randomized methods, we use an oversampling parameter equal to 10 (i.e., in Algorithms 2 and 3).
6.1 Matrix ID Experiments
We compare the four methods in Table 1. For standard matrix ID, we use the implementation in RSVDPACK111Available at
https://github.com/sergeyvoronin/LowRankMatrixDecompositionCodes.. For the remaining methods, we use our own Matlab implementations which utilize Matlab’s column pivoted QR function. RSVDPACK only supports dense matrices. Moreover, since it is challenging to efficiently construct partial QR decompositions of sparse matrices, we did not attempt to write our own implementation of standard matrix ID for sparse matrices; see Section 11.1.8 of Golub and Van Loan (2013) for a discussion about the challenges of sparse QR. We therefore have to convert each sparse input matrix to dense format before applying standard matrix ID from RSVDPACK. Similarly, it is challenging to implement an efficient algorithm for FFT for sparse matrices, or for the accelerated FFT by Woolfe et al. (2008). We therefore also have to convert the input matrix to dense format before applying standard FFT in our implementation of SRFT matrix ID. However, by only sketching a subset of columns of the input matrix at a time, we can avoid having to convert all columns of the matrix to dense format at the same time. In the experiments, we use the modification described in Remark 2 when implementing our proposed CountSketch matrix ID.
Computing the spectral norm of the matrices we consider is not feasible due to their size. Therefore, when computing the error for each matrix decomposition, we utilize the randomized scheme for estimating the spectral norm suggested in Section 3.4 of Woolfe et al. (2008). Letting be the true error in spectral norm, our estimates satisfy the following properties: , and where . In other words, the estimate is smaller than the true spectral norm, but it is unlikely to be much smaller (with “much smaller” meaning more than two orders of magnitude smaller). This is good enough for our purposes, since we are primarily interested in comparing the performance of the different methods rather than establishing the exact errors. In the first experiment, , and in the second .
6.1.1 Experiment 1: Synthetic Matrices
We generate sparse matrices with , and density . We use different values of . The matrices have a true rank of , where . Similarly to experiments by Martinsson et al. (2011), we let , , decay exponentially to , and then remain constant at , .
The results for the first experiment are presented in Figure 1. Standard matrix ID encountered memory issues when . For the matrix sizes the standard method could handle, it was more accurate but much slower than the randomized methods. The accuracy of all randomized methods is comparable. Our proposed CountSketch matrix ID is the fastest, achieving a speed-up of about and when compared to Gaussian and SRFT ID, respectively.
6.1.2 Experiment 2: Real-World Matrices
We decompose a sparse matrix which comes from a computer vision problem and is part of the SuiteSparse Matrix Collection222The matrix can be downloaded from https://sparse.tamu.edu/Brogan/specular.. The matrix is of size 477,976 by 1,600, contains 7,647,040 nonzero elements, and has a rank of 1,442. We set the target rank to . Ideally, the methods should be able to produce decompositions with a very small error. We only attempt this with the three randomized methods, since the matrix is too large for standard matrix ID. Table 3 shows the result. All methods produce good approximations with a small error. Our proposed CountSketch matrix ID method is much faster than the other algorithms, achieving a speed-up of about and compared to Gaussian and SRFT matrix ID, respectively.
To further support our claim that works well in practice, we have done additional experiments. We consider 20 matrices from the SuiteSparse Matrix Collection333They are landmark, Franz7, ch7-8-b2, ch7-9-b2, ch8-8-b2, mk12-b2, shar_te2-b1, rel7, relat7b, relat7, abtaha2, abtaha1, specular, photogrammetry2, GL7d12, ch7-6-b2, ch7-7-b2, cis-n4c6-b3, mk11-b2, n4c6-b3. , and 3 different target ranks (10%, 50% and 90% of the number of columns). The matrices are of different sizes and come from different application areas. We compare the performance of Gaussian, SRFT and CountSketch (proposed method) matrix ID, all with , repeating each experiment 10 times and reporting averages. The results are in Table 4; on average, our method is the most accurate even compared to the Gaussian method which it outperforms in 34 of the 60 tests. Our method outperforms the SRFT method in 57 of the 60 tests. Out of the 26 cases when the Gaussian method outperforms our method, the difference is no more than 7% in 25 of those cases, and 31% in one case. Out of the 58 cases when the Gaussian method outperforms the SRFT method, the difference is no more than 14% in 56 of those cases, and 36%–45% in two cases.
As mentioned in Remark 1, semi-coherent matrices are adversarial to CountSketch and our proposed method. The matrix soc-sign-bitcoin-otc in the SuiteSparse Matrix Collection, which is the adjacency matrix of a graph, is a concrete example of when choosing results in a large error for our method. The semi-coherent structure of this matrix can be revealed by rearranging it so that rows and columns corresponding to nodes in the same strongly connected components of the graph are adjacent in the matrix. Some care is therefore necessary when applying our method together with the rule of thumb .
6.2 Tensor ID Experiments
We compare the three methods in Table 2. We have implemented all methods ourselves in Matlab and C.
6.2.1 Experiment 1: Synthetic Tensors
We generate sparse -way tensors using (1), where each factor matrix column is a random sparse vector with a density of 1%, and we use different values of . The number of rank-1 terms is , and we use a target rank of 1,000. The values of in (1) are defined as for , and for . The results for the experiment are presented in Figure 2. Gaussian and CountSketch tensor ID achieve similar accuracy. Although the Gram matrix approach has a better accuracy here, it can have issues reaching an error below the square root of machine precision due to poor conditioning; see the example in Section 5.1.1 of Biagioni et al. (2015). Our proposed method is much faster than both other methods for the larger tensors, achieving a speed-up of over the Gram matrix approach (for ) and over Gaussian tensor ID (for ).
6.2.2 Experiment 2: Real-World Tensor
The purpose of this experiment is to show how tensor ID can be useful in a data analysis task. We implement Algorithm 2 by Reynolds et al. (2017), which requires repeated rank reduction, and use it to find the maximum magnitude element in a CP tensor which comes from decomposing streamed data. The rank reduction step is done using tensor ID. The data we consider is a decomposed version of the Enron data set444The data set is available at http://frostt.io/tensors/enron. of size . The Enron data set keeps track of email correspondence between employees at Enron, and the four modes represent sender, receiver, keyword and date. The decomposition has rank 100, and was constructed using the streamed version of SPLATT555The streamed version of SPLATT is available at https://github.com/ShadenSmith/splatt-stream. (Smith et al., 2018), with the data streamed along the fourth mode (time). As suggested in the documentation of SPLATT-stream, we apply an additional Frobenius norm regularizer with regularization coefficient 1e−2 to the mode-4 factor matrix. We threshold the factor matrices outputted by SPLATT-stream by first normalizing them so that each column have unit 2-norm (the normalization constant is absorbed into the s-values) and then setting all elements with magnitude less than 1e−6 to zero. The relative error introduced by this thresholding is less than 2e−5.
Unlike the previous experiments, the matrices being sketched in this experiment have many rows containing only zeros. We therefore could speed up Gaussian tensor ID by only generating those columns of the Gaussian sketch matrices which are actually multiplied by nonzero elements. We used this improved version of Gaussian tensor ID in the experiment for a more fair comparison. The same modification does not yield a speed-up of Gaussian matrix or tensor ID in the previous experiments since there most rows of the matrices being sketched contain nonzero elements.
Finding the maximum magnitude element using a brute force approach would require computing every nonzero element in the tensor, which would be costly. Using the algorithm by Reynolds et al. (2017) together with our CountSketch tensor ID, we find the maximum in 11 seconds. The sketching portion of the algorithm takes more time if Gaussian tensor ID is used instead. We do not compare with the Gram matrix approach since it takes very long to run. With the results in the previous subsection in mind, we believe the speed-up would be more substantial for higher rank tensors. For all ten trials, and both when using CountSketch and Gaussian tensor ID for rank reduction, the same position for the maximum magnitude element is identified each time.
7 Proofs
7.1 Proof of Proposition 1
Our proof of Proposition 1 is an adaption of the proof for SRFT matrix ID provided by Woolfe et al. (2008). We show that their arguments hold when a CountSketch matrix is used for sketching instead of an SRFT matrix. Although much of our proof is identical to that provided by Woolfe et al. (2008), we choose to include it in detail. The reason for doing this is that the proofs of Propositions 2 and 4 rely on adapting the proof in the present section. Having a detailed proof here therefore makes those subsequent proofs easier to follow.
The following facts will be useful in the proof.
Fact 2 (Lemma 3.7 in Martinsson et al. (2011))
Let and be positive integers with . Suppose is a matrix such that is invertible. Then
[TABLE]
Fact 3 (Lemma 3.7 in Woolfe et al. (2008))
Let , , and be positive integers such that . Suppose , is a matrix whose columns constitute a subset of the columns of , , , and . Then
[TABLE]
Fact 4 (Lemma 3.9 in Martinsson et al. (2011))
Let , and be positive integers. Suppose , and . Then for all .
Fact 5 is a special case of a more general statement in Atkinson and Han (2009).
Fact 5 (Theorem 2.3.1 of Atkinson and Han (2009))
Let be a matrix and assume . Then is invertible and
[TABLE]
Lemma 1 is an adaption of Lemma 4.2 by Woolfe et al. (2008).
Lemma 1
Let , and be positive integers such that . Suppose is a CountSketch matrix, and is a matrix with orthonormal columns. Define as
[TABLE]
and define elementwise as
[TABLE]
Then .
Proof
For ,
[TABLE]
Since
[TABLE]
we can rewrite (17) as
[TABLE]
The second term on the last line in the equation above is just . Since
[TABLE]
and , the first term is just
[TABLE]
It follows that . ∎
Lemma 2 is an adaption of Lemma 4.3 by Woolfe et al. (2008).
Lemma 2
Let and be real numbers such that , and let , and be positive integers such that
[TABLE]
Suppose is a CountSketch matrix, is a matrix with orthonormal columns, and is the matrix defined in (16). Then
[TABLE]
with probability at least .
Proof
Using the definition in (16), we have
[TABLE]
Note that for each term in the sum above, and . This means that unless ( and ) or ( and ), we have
[TABLE]
since each is independent from all other random variables, and since for all . We can therefore rewrite (24) as
[TABLE]
The matrix has exactly one nonzero entry which is equal to 1 in each column. Consequently,
[TABLE]
The event happens with probability when . If follows that
[TABLE]
Using this fact, and the fact that each , (26) simplifies to
[TABLE]
Note that
[TABLE]
Moreover,
[TABLE]
Combining (29), (30) and (LABEL:eq:lemma-2-5) yields
[TABLE]
Since
[TABLE]
we have
[TABLE]
Using Markov’s inequality and the condition in (22), we have
[TABLE]
Consequently,
[TABLE]
∎
Lemma 3 is an adaption of Lemma 4.4 by Woolfe et al. (2008).
Lemma 3
Let , , , and satisfy the same properties as in Lemma 2. Furthermore, suppose , , and are defined as in Lemma 2, and let . If (23) is true, then the following hold:
[TABLE]
* is invertible, and*
[TABLE]
Proof
Using Lemma 1 and (23), we then have
[TABLE]
Since and , it follows from Fact 5 that is invertible and
[TABLE]
where the last inequality follows from (23). Consequently,
[TABLE]
∎
Lemma 4 is an adaption of Lemma 4.5 by Woolfe et al. (2008).
Lemma 4
Let and be positive integers with . Suppose is a CountSketch matrix. Then .
Proof
The matrix contains nonzero elements, all of magnitude 1. It follows that , and hence . ∎
Lemma 5 is an adaption of Lemma 4.6 by Woolfe et al. (2008).
Lemma 5
Let , , , and satisfy the same properties as in Lemma 2. Suppose is a CountSketch matrix, and is an arbitrary matrix. Then, with probability at least , there exists a matrix such that
[TABLE]
and
[TABLE]
Proof
Let be the SVD of , where and are unitary, and is diagonal with non-negative entries. Split into two matrices and so that . Let and . Then
[TABLE]
Define and let be the corresponding matrix defined in (16), but in terms of instead of . Then according to Lemma 1. For the remainder of the proof, we will assume that , which happens with probability at least according to Lemma 2. Then is invertible according to Lemma 3. Define and
[TABLE]
According to Fact 2 and Lemma 3, it follows that
[TABLE]
Combining (45) and (46), we have
[TABLE]
So (43) is satisfied. Next, let and be the matrices in the upper left and lower right corners of , respectively, so that
[TABLE]
It is easy to verify that
[TABLE]
Using (48), we can further rewrite
[TABLE]
Note that
[TABLE]
From (48), we know that
[TABLE]
Moreover, using (44), the fact that is unitary, and Lemma 4, we have
[TABLE]
Combining (49), (50), (51), (52), (53) and (46) we have
[TABLE]
which proves (42). ∎
We can now prove Proposition 1 in the main manuscript. The proof is an adaption of the discussion in Section 5.1 of Woolfe et al. (2008).
Proof (Proposition 1)
According to Fact 1, the outputs and computed on line 5 of Algorithm 2 satisfy the following: satisfies properties (i)–(iv) in Fact 1, including
[TABLE]
and
[TABLE]
since . Applying Fact 3, we have
[TABLE]
where is an arbitrary matrix. From Lemma 5, with probability at least , we can choose such that the bounds in (42) and (43) hold. Moreover, since , it follows that , and consequently,
[TABLE]
Combining (42), (43), (55), (56), (57), (58), and Fact 4 gives that
[TABLE]
with probability at least . Setting then yields the same bounds as in the statement in Proposition 1. ∎
7.2 Formal Statement and Proof of Claim in Remark 2
We express the statement in Remark 2 in slightly different terms here. Let be a hybrid deterministic/random function defined as
[TABLE]
where all are iid random variables that are uniformly distributed in . Furthermore, let be a uniform random permutation function. We then define as . Using instead of in the definition of CountSketch ensures that is of full rank. The guarantees of Proposition 1 still hold for this modified CountSketch, and in fact the bound in (9) is slightly improved.
Proposition 4
If defined in this way is used instead of when defining on line 3 in Algorithm 2, then Proposition 1 still holds, but with the condition in (9) improved to
[TABLE]
We have not seen anyone else consider this kind of modified CountSketch.
Proof (Proposition 4)
When using the modified CountSketch matrix proposed in Remark 2, the only thing that will change in the proof in Section 7.1 is Lemma 2. Notice that going from to only impacts and not , and since and remain independent, the argument that takes us from (24) to (26) remains valid for . Indeed, the key conditions that each is independent from all other random variables and remain true when we use instead of . However, the expectation in (28) will change, which impacts (29), due to the fact that the probability of the event when is not . Note that
[TABLE]
We can rewrite
[TABLE]
Notice that the first term on the right hand side of (63) is zero, since then will map and to distinct elements. The second and third term in (63) are equal. Considering the second term, we have
[TABLE]
since if , then if and only if . Furthermore,
[TABLE]
where the second equality is true since each , is iid. For the fourth term in the right hand side of (63), we have
[TABLE]
where the second equality again holds since each , is iid and . Combining (62), (63), (64), (65), (66), and using the fact that the second and third term in (63) are equal, we get
[TABLE]
Proceeding with the remainder of the proof of Lemma 2 as before, we now get a bound
[TABLE]
Using this new bound and the new condition
[TABLE]
together with Markov’s inequality, we get
[TABLE]
and consequently
[TABLE]
holds in this case too.
All the other lemmas will remain the same, with the only exception that Lemmas 3 and 5 now will use the new condition in (69) instead of the old one in (22). The proof of the proposition itself at the end of Section 7.1 will therefore remain identical. When using the modified CountSketch, the statements in Proposition 1 will therefore remain true with the new condition in (69). Setting then yields the desired bound. ∎
7.3 Proof of Proposition 2
The following fact will be useful.
Fact 6 (Theorem B.1 in the supplement of Diao et al. (2018))
Recall that . Let be a TensorSketch operator defined as in Section 1.3 in terms of CountSketch operators.
- (i)
Suppose and are matrices with rows. For , we have
[TABLE] 2. (ii)
Suppose is any matrix. If , then the following holds with probability at least :
[TABLE]
We break the proof into two parts. First, we prove Lemma 6 which is a variant of Proposition 1 for the case when a TensorSketch operator is used instead of a CountSketch operator. Then we prove the proposition itself.
Lemma 6
Let and be real numbers such that , and let , , and be positive integers such that and
[TABLE]
Suppose that the matrix ID on line 6 of Algorithm 3 utilizes SRRQR. Then the outputs and on that line will satisfy
[TABLE]
with probability at least .
Proof
Recall from Section 1.3 that TensorSketch is defined similarly to CountSketch, but using the hash function instead of , and using the diagonal matrix instead of . Letting be a matrix with for , and with all other entries equal to 0, we can write . This means that the proof in Section 7.1 largely can be repeated to prove the present lemma. Lemma 1 remains true in its present form when TensorSketch is used instead of CountSketch.
To see that Lemma 2 remains true with the new condition when is replaced by , let be a TensorSketch operator, and let be a matrix with orthonormal columns. Define , and let be defined as in (16), but in terms of the corresponding quantities from TensorSketch. Then , according to Lemma 1. Using Fact 6 (i), condition (74), and the fact that , we have
[TABLE]
All the other lemmas will remain the same when is used instead of , with the only exception that Lemmas 3 and 5 now will use the new condition in (74) instead of the old one in (22). Using exactly the same arguments as in the proof of Proposition 1 at the end of Section 7.1 will therefore give the bound in (75). ∎
Proof (Proposition 2)
Recall that and are defined as in (1) and (2), respectively, and the coefficients are defined as
[TABLE]
for . We then have
[TABLE]
Letting , we have
[TABLE]
where the inequality follows from Cauchy–Schwarz inequality. Combining (78) and (79) we get
[TABLE]
where the second inequality is a well-known relation (see e.g. equation (2.3.7) in Golub and Van Loan (2013)). Combining (80) and Lemma 6 gives that
[TABLE]
with probability at least . Setting then yields the same bounds as in the statement in Proposition 2. ∎
7.4 Proof of Proposition 3
Proof
Note that is of size , with . So is the smallest singular value of . Suppose
[TABLE]
To simplify notation, let . Using Theorem 8.6.1 in Golub and Van Loan (2013) and Fact 6 (ii), we have that with probability at least , the following hold:
[TABLE]
and
[TABLE]
We therefore have
[TABLE]
with probability at least . Setting gives us the bounds in Proposition 3. ∎
8 Conclusion
We have presented a new fast randomized algorithm for computing matrix ID, which utilizes CountSketch. We have then shown how this method can be extended to computing the tensor ID of CP tensors. For both the matrix and tensor settings, we provided performance guarantees. To the best of our knowledge, we provide the first performance guarantees for any randomized tensor ID algorithm. We conducted several numerical experiments on both synthetic and real data. These experiments showed that our algorithms maintain the same accuracy as other randomized methods, but with a much shorter run time, running at least an order of magnitude faster on the larger matrices and tensors.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Ailon and Chazelle (2009) Nir Ailon and Bernard Chazelle. The fast Johnson–Lindenstrauss transform and approximate nearest neighbors. SIAM Journal on Computing , 39(1):302–322, 2009.
- 2Atkinson and Han (2009) Kendall Atkinson and Weimin Han. Theoretical Numerical Analysis: A Functional Analysis Framework . Number 39 in Texts in Applied Mathematics. Springer-Verlag, New York, 3rd edition, 2009. ISBN 978-1-4419-0457-7.
- 3Avron et al. (2010) Haim Avron, Petar Maymounkov, and Sivan Toledo. Blendenpik: Supercharging LAPACK’s least-squares solver. SIAM Journal on Scientific Computing , 32(3):1217–1236, 2010.
- 4Avron et al. (2014) Haim Avron, Huy L. Nguyen, and David P. Woodruff. Subspace Embeddings for the Polynomial Kernel. In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2 , pages 2258–2266, Cambridge, MA, USA, 2014. MIT Press.
- 5Battaglino et al. (2018) Casey Battaglino, Grey Ballard, and Tamara G. Kolda. A practical randomized CP tensor decomposition. SIAM Journal on Matrix Analysis and Applications , 39(2):876–901, 2018.
- 6Beylkin and Mohlenkamp (2002) Gregory Beylkin and Martin J. Mohlenkamp. Numerical operator calculus in higher dimensions. Proceedings of the National Academy of Sciences , 99(16):10246–10251, 2002.
- 7Beylkin and Mohlenkamp (2006) Gregory Beylkin and Martin J. Mohlenkamp. Algorithms for Numerical Analysis in High Dimensions. SIAM Journal on Scientific Computing , 26(6):2133–2159, July 2006.
- 8Biagioni et al. (2015) David J. Biagioni, Daniel Beylkin, and Gregory Beylkin. Randomized interpolative decomposition of separated representations. Journal of Computational Physics , 281(C):116–134, January 2015. ISSN 0021-9991. doi: 10.1016/j.jcp.2014.10.009 .
