Spectrum Approximation Beyond Fast Matrix Multiplication: Algorithms and Hardness
Cameron Musco, Praneeth Netrapalli, Aaron Sidford, Shashanka Ubaru,, David P. Woodruff
TL;DR
This paper develops algorithms to efficiently approximate the singular value spectrum of matrices and related norms faster than traditional methods, revealing fundamental links to matrix multiplication complexity.
Contribution
It introduces the first algorithms for spectrum approximation and symmetric matrix norm estimation that run faster than matrix multiplication time, using novel techniques from stochastic trace estimation and polynomial approximation.
Findings
Efficient spectrum approximation algorithms for matrices in subcubic time.
Faster algorithms for approximating Schatten-$p$ norms, including nuclear norm.
Complexity of spectrum approximation is inherently tied to fast matrix multiplication.
Abstract
Understanding the singular value spectrum of a matrix is a fundamental task in countless applications. In matrix multiplication time, it is possible to perform a full SVD and directly compute the singular values . However, little is known about algorithms that break this runtime barrier. Using tools from stochastic trace estimation, polynomial approximation, and fast system solvers, we show how to efficiently isolate different ranges of 's spectrum and approximate the number of singular values in these ranges. We thus effectively compute a histogram of the spectrum, which can stand in for the true singular values in many applications. We use this primitive to give the first algorithms for approximating a wide class of symmetric matrix norms in faster than matrix multiplication time. For example, we give a …
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.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Spectrum Approximation Beyond Fast Matrix Multiplication:
Algorithms and Hardness
Cameron Musco
MIT
Praneeth Netrapalli
Microsoft Research
Aaron Sidford
Stanford University
Shashanka Ubaru
University of Minnesota
David P. Woodruff
Carnegie Mellon University
Understanding the singular value spectrum of a matrix is a fundamental task in countless numerical computation and data analysis applications. In matrix multiplication time, it is possible to perform a full SVD of and directly compute the singular values . However, little is known about algorithms that break this runtime barrier.
Using tools from stochastic trace estimation, polynomial approximation, and fast linear system solvers, we show how to efficiently isolate different ranges of ’s spectrum and approximate the number of singular values in these ranges. We thus effectively compute an approximate histogram of the spectrum, which can stand in for the true singular values in many applications.
We use our histogram primitive to give the first algorithms for approximating a wide class of symmetric matrix norms and spectral sums faster than the best known runtime for matrix multiplication. For example, we show how to obtain a approximation to the Schatten -norm (i.e. the nuclear or trace norm) in just time for with uniform row sparsity or time for dense matrices. The runtime scales smoothly for general Schatten -norms, notably becoming for any real .
At the same time, we show that the complexity of spectrum approximation is inherently tied to fast matrix multiplication in the small regime. We use fine-grained complexity to give conditional lower bounds for spectrum approximation, showing that achieving milder dependencies in our algorithms would imply triangle detection algorithms for general graphs running in faster than state of the art matrix multiplication time. This further implies, through a reduction of [WW10], that highly accurate spectrum approximation algorithms running in subcubic time can be used to give subcubic time matrix multiplication. As an application of our bounds, we show that precisely computing all effective resistances in a graph in less than matrix multiplication time is likely difficult, barring a major algorithmic breakthrough.
1 Introduction
Given , a central primitive in numerical computation and data analysis is to compute ’s spectrum: the singular values . These values can reveal matrix structure and low effective dimensionality, which can be exploited in a wide range of spectral data analysis methods [Jol02, US16]. The singular values are also used as tuning parameters in many numerical algorithms performed on [GVL12], and in general, to determine some of the most well-studied matrix functions [Hig08]. For example, for any , we can define the spectral sum:
[TABLE]
Spectral sums often serve as snapshots of ’s spectrum and are important in many applications. They encompass, for example, the log-determinant, the trace inverse, the Schatten -norms, including the important nuclear norm, and general Orlicz norms (see Section 1.2 for details).
While the problem of computing a few of the largest or smallest singular values of has been exhaustively studied [Par98, Saa11], much less is known about algorithms that approximate the full spectrum, and in particular, allow for the computation of summary statistics such as spectral sums. In time, it is possible to perform a full SVD and compute the singular values exactly.111 Note that an exact SVD is incomputable even with exact arithmetic [TB97]. Nevertheless, direct methods for the SVD obtain superlinear convergence rates and hence are often considered to be ‘exact’. Here, and throughout, denotes the current best exponent of fast matrix multiplication [Wil12]. However, even if one simply desires, for example, a constant factor approximation to the nuclear norm , no time algorithm is known. We study the question of spectrum approximation, asking whether obtaining an accurate picture of ’s spectrum is truly as hard as matrix multiplication, or if it is possible to break this barrier. We focus on spectral sums as a motivating application.
1.1 Our Contributions
Upper Bounds:
On the upper bound side, we show that significant information about ’s spectrum can be determined in time, for the current value of . We show how to compute a histogram of the spectrum, which gives approximate counts of the number of squared singular values in the ranges for some width parameter and for ranging from [math] to some maximum . Specifically our algorithm satisfies the following:
Theorem 1** (Histogram Approximation – Informal).**
Given , let be the number of squared singular values of on the range . Then given error parameter , with probability , Algorithm 1 outputs for all , satisfying:
[TABLE]
*For input parameter , let and . Let be the maximum number of nonzeros in a row of . The algorithm’s runtime can be bounded by: *
[TABLE]
for sparse or for dense , where is the time it takes to multiply a matrix by a matrix using fast matrix multiplication.
This primitive is useful on its own, giving an accurate summary of ’s spectrum which can be used in many downstream applications. Setting the parameter appropriately to balance costs (see overview in Section 1.3), we use it to give the first algorithms for computing relative error approximations to a broad class of spectral sums for functions , which are a) smooth and b) quickly growing, so that very small singular values cannot make a significant contribution to . This class includes for example the Schatten -norms for all , the SVD entropy, the Ky Fan norms, and many general Orlicz norms.
For a summary of our -norm results see Table 1. Focusing for simplicity on square matrices, with uniformly sparse rows, and assuming , are constants, our algorithms approximate in time for any real .222For any , denotes the number of nonzero entries in . For , we achieve runtime. In the important case of , this becomes . Note that , and for sparse enough , this bound is subquadratic. For dense , we use fast matrix multiplication, achieving time for all . For , this gives . Even without fast matrix multiplication, the runtime is , and so for .
Lower Bounds:
On the lower bound side, we show that obtaining time spectrum approximation algorithms with very high accuracy may be difficult. Our runtimes all depend polynomially on the error , and we show that improving this dependence, e.g., to , or even to a better polynomial, would give faster algorithms for the well studied Triangle Detection problem.
Specifically, for a broad class of spectral sums, including all Schatten -norms with , SVD entropy, , , and , we show that any approximation algorithm running in time yields an algorithm for triangle detection running in time. For and sufficiently small , such an algorithm would improve the state of the art in triangle detection, which currently requires time on dense graphs. Furthermore, through a reduction of [WW10], any subcubic time triangle detection algorithm yields a subcubic time algorithm for Boolean Matrix Multiplication (BMM). Thus, any spectral sum algorithm achieving subcubic runtime and accuracy for small enough constant , must (implicitly) implement fast matrix multiplication. This is in stark contrast to the fact that, for , for many spectral sums, including all Schatten -norms with , we obtain subcubic, and in fact for , runtimes without using fast matrix multiplication (see Table 1 for precise dependencies).
Our lower bounds hold even for well-conditioned matrices and structured matrices like symmetric diagonally dominant (SDD) systems, both of which admit nearly linear time algorithms for system solving [ST04]. This illustrates a dichotomy between linear algebraic primitives like applying to a vector and spectral summarization tasks like precisely computing , which in some sense require more global information about the matrix. Our analysis has ramifications regarding natural open problems in graph theory and numerical computation. For example, for graph Laplacians, we show that accurately computing all effective resistances yields an accurate algorithm for computing of certain matrices, which is enough to give triangle detection.
1.2 Related Work on Spectral Sums
The applications of approximate spectral sum computation are very broad. When is positive semidefinite (PSD) and , is the log-determinant, which is important in machine learning and inference applications [Ras04, DKJ*+*07, FHT08]. For , is the trace of the inverse, used in uncertainty quantification [BCF09] and quantum chromodynamics [SLO13].
When , where is the Schatten -norm of . Computation of the Schatten -norm, also known as the nuclear or trace norm, is required in a wide variety of applications. It is often used in place of the matrix rank in matrix completion algorithms and other convex relaxations of rank-constrained optimization problems [CR12, DTV11, JNS13, NNS*+*14]. It appears as the ‘graph energy’ in theoretical chemistry [Gut92, Gut01], the ‘singular value bound’ in differential privacy [HLM12, LM12], and in rank aggregation and collaborative ranking [LN15].
Similar to the nuclear norm, general Schatten -norms are used in convex relaxations for rank-constrained optimization [NHD12, NWC*+*12]. They also appear in image processing applications such as denoising and background subtraction [XGL*+*16], classification [LYCG14], restoration [XQT*+*16], and feature extraction [DHJZ15].
When (after is normalized by ), is the SVD entropy [ABB00], which is used in feature selection [VGLH06, BP14], financial data analysis [Car14, GXL15], and genomic data [ABB00] applications.
Despite their importance, prior to our work, few algorithms for fast computation of spectral sums existed. Only a few special cases of the Schatten -norms were known to be computable efficiently in time. These include the Frobenius norm () which is trivially computed in time, the spectral norm () which can be estimated via the Lanczos method in time [KW92], and the Schatten- norms for even integers , or general integers with PSD . These norms can be approximated in time via trace estimation [Woo14, BCKY16], since when is even or is PSD, is PSD and so its trace equals .
There are a number of works which consider estimating matrix norms in sublinear space and with a small number of passes over [LNW14, AKR15, LW16a, LW16b, BCKY16, LW17]. However, in these works, the main focus is on space complexity, and no non-trivial runtime bounds are given. We seem to be the first to tackle the arguably more fundamental problem of obtaining the best time complexity for simple norms like the Schatten- norms, irrespective of space and pass complexity.
Another interesting line of works tries to estimate the Schatten- norms of an underlying covariance matrix from a small number of samples from the distribution [KV16], or from entrywise sampling under various incoherence assumptions [KO17]. This model is different from ours, as we do not assume an underlying distribution or any incoherence properties. Moreover, even with such assumptions, these algorithms also only give non-trivial sample complexity when either is PSD and is an integer, or is a general matrix but is an even integer, which as mentioned above are easy to handle from the perspective of time complexity alone.
A number of works have focused on computing spectral sums when has bounded condition number, and relative error results exist for example for the log-determinant, , , and the Schatten -norms [BDKZ15, HMAS17, UCS17]. We are the first to give relative error results in time for general matrices, without the condition number dependency. Our histogram approach resembles spectral filtering and spectral density estimation techniques that have been considered in the numerical computation literature [ZWJ15, DNPS16, LSY16, US16, UCS17, USS17]. However, this literature typically requires assuming gaps between the singular values and existing work is not enough to give relative error spectral sum approximation for general matrices
1.3 Algorithmic Approach
Spectral Sums via Trace Estimation:
A common approach to spectral sum approximation is to reduce to a trace estimation problem involving the PSD matrix , using the fact that the trace of this matrix equals the sum of its singular values. In fact, this has largely been the only known technique, other than the full SVD, for obtaining aggregate information about ’s singular values [Hut90, SLO13, WLK*+*16, WWZ14, RKA15, FORF16, BDKZ15, HMAS17]. The idea is, letting , we have Writing the SVD , and defining the matrix function where , we have since, if is nonnegative, is PSD and its trace equals the sum of its singular values.
It is well known that this trace can be approximated up to accuracy by averaging samples of the form where is a random Gaussian or sign vector [Hut90, AT11]. While cannot be explicitly computed without a full SVD, a common approach is to approximate with a low-degree polynomial [BDKZ15, HMAS17]. If has degree , one can apply to any vector in time, and so estimate its trace in just time. Unfortunately, for many of the functions most important in applications, e.g., for odd , , , has a discontinuity at and cannot be approximated well by a low-degree polynomial near zero. While the approximation only needs to be good on the range , the required degree will still typically depend on where is the condition number, which can be unbounded for general matrices.
Singular Value Deflation for Improved Conditioning:
Our first observation is that, for many functions, it is not necessary to approximate on the full spectral range. For example, for (i.e., when ), setting , we have:
[TABLE]
Hence we can safely ‘ignore’ any and still obtain a relative error approximation to . The larger is, the larger we can set (corresponding to in Theorem 1) to be, since, after powering, the singular values below this threshold do not contribute significantly to . For , our ‘effective condition number’ for approximating becomes . Unfortunately, in the worst case, we may have and hence . Hiding dependences, this gives runtime when .
To improve the effective condition number, we can apply singular vector deflation. Our above bound on is only tight when the first singular value is very large and so dominates . We can remedy this by flattening ’s spectrum by deflating off the top singular vectors (corresponding to in Theorem 1), and including their singular values in the spectral sum directly.
Specifically, letting be the projection onto the top singular vectors of , we consider the deflated matrix , which has . Importantly, , and thus this singular value cannot dominate the -norm. As an example, considering and ignoring dependencies, our effective condition number after deflation is
[TABLE]
The runtime required to approximate via an iterative method (ignoring possible gains from fast matrix multiplication) is roughly . We then require time to approximate the polynomial trace of . The term comes from projecting off the top singular directions with each application of . Setting to balance the costs, we obtain runtime .
For a similar argument gives runtime . This is already a significant improvement over a full SVD. As grows larger, the runtime approaches reflecting the fact that for larger we can ignore a larger and larger portion of the small singular values in and correspondingly deflate off fewer and fewer top values.
Unfortunately, we get stuck here. Considering the important Schatten- norm, for a matrix with singular values each equal to and singular values each equal to , the tail of small singular values contributes a constant fraction of . However, there is no good polynomial approximation to on the range with degree (recall that we pick this function since ). So to accurately approximate , we either must deflate off all top singular values, requiring time, or apply a degree polynomial approximation, requiring the same amount of time.
Further Improvements with Stochastic Gradient Descent:
To push beyond this barrier, we look to stochastic gradient methods for linear systems. When using polynomial approximation, our bounds depend on the condition number of the interval over which we must approximate , after ignoring the smallest singular values and deflating off the largest. This is analogous to the condition number dependence of iterative linear system solvers like conjugate gradient or accelerated gradient descent, which approximate for using a polynomial of .
However, recent advances in convex optimization offer an alternative. Stochastic gradient methods [JZ13, SSZ14] sample one row, , of at a time, updating the current iterate by adding a multiple of . They trade a larger number of iterations for updates that take time, rather than time to multiply by a vector. At the same time, these methods give much finer dependencies on the singular value spectrum. Specifically, it is possible to approximately apply to a vector with the number of iterations dependent on the average condition number:
[TABLE]
is always at most the standard condition number, . It can be significantly smaller when has a quickly decaying spectrum, and hence . Further, the case of a quickly decaying spectrum with a few large and many small singular values is exactly the hard case for our earlier approach. If we can understand how to translate improvements on linear system solvers to spectral sum approximation, we can handle this hard case.
From Linear System Solvers to Histogram Approximation:
The key idea to translating the improved average condition number bounds for linear systems to our problem of approximating is to note that linear system solvers can be used to apply threshold functions to .
Specifically, given any vector , we can first compute . We can then apply a fast system solver to approximate . The matrix function has a number of important properties. All its singular values are between [math] and . Further, any singular value in with value is mapped to a singular value in which is . Correspondingly, any singular value is mapped to .
Thus, we can apply a low degree polynomial approximation to a step function at to to obtain , which approximates a step function at [FMMS16]. For some steepness parameter which affects the degree of the polynomial approximation, for we have and for , . On the intermediate range , falls somewhere between [math] and .
By composing these approximate threshold functions at different values of , it is possible to ‘split’ our spectrum into a number of small spectral windows. For example, is on the range and outside this range, with some ambiguity near and .
Splitting our spectrum into windows of the form for a width parameter , and applying trace estimation on each window lets us produce an approximate spectral histogram. Of course, this histogram is not exact and in particular, the ‘blurring’ of our windows at their boundaries can introduce significant error. However, by applying a random shifting technique and setting the steepness parameter small enough (i.e., ), we can ensure that most of the spectral weight falls outside these boundary regions with good probability, giving Theorem 1.
From Histogram Approximation to Spectral Sums:
If is small enough, and (the function in the spectral sum) and correspondingly (where ) are smooth enough, we can approximate by simply summing over each bucket in the histogram, approximating by its value at one end of the bucket.
This technique can be applied for any spectral sum. The number of windows required (controlled by ) and the histogram accuracy scale with the smoothness of and the desired accuracy in computing the sum, introducing polynomial runtime dependencies on these parameters.
However, the most important factor determining the final runtime is the smallest value (corresponding to in Theorem 1) which we must include in our histogram in order to approximate . The cost of computing the last bucket of the histogram is proportional to the cost of applying , and hence of approximately computing . Using stochastic gradient descent this depends on the average condition number of .
Again considering the Schatten -norm for illustration, we can ignore any singular values with . Hiding dependence, this means that in our histogram, we must include any singular values of with value . This gives us effective average condition number after deflating off the top singular values:
[TABLE]
where the last inequality follows from the observation that and . Comparing to (1), this bound is better by an factor.
Ignoring details and using a simplification of the runtimes in Theorem 1, we obtain an algorithm running in time to deflate singular vectors, plus time to approximate the spectral sum over the deflated matrix. Choosing to balance these costs, gives our final runtimes. For the nuclear norm, using the bound on from (2), we set which gives and runtime where is the maximum row sparsity. For dense this is , which is faster than state of the art matrix multiplication time. It can be further accelerated using fast matrix multiplication methods. See details in Section 7.
Returning to our original hard example for intuition, we have with singular values at and singular values at . Even without deflation, we have (again ignoring dependencies) Since and , this gives Thus, we can actually approximate in just time for this matrix.
With average condition number dependence, our performance is limited by a new hard case. Consider with singular values at and at . The average condition number without deflation is giving . Further, we can see that unless we deflate off nearly all top singular vectors, we do not improve this bound significantly.
1.4 Lower Bound Approach
We now shift focus to our lower bounds, which explore the fundamental limits of spectrum approximation using fine-grained complexity approaches. Fine-grained complexity has had much success for graph problems, string problems, and problems in other areas (see, e.g., [Wil15] for a survey), and is closely tied to understanding the complexity of matrix multiplication. However, to the best of our knowledge it has not been applied broadly to problems in linear algebra.
Existing hardness results for linear algebraic problems tend to apply to restricted computational models such as arithmetic circuits [BS83], bilinear circuits or circuits with bounded coefficients and number of divisions [Mor73, Raz03, Shp03, RS03], algorithms for dense linear systems that can only add multiples of rows to each other [KKS65, KS70], and algorithms with restrictions on the dimension of certain manifolds defined in terms of the input [Win70, Win87, Dem13]. In contrast, we obtain conditional lower bounds for arbitrary polynomial time algorithms by showing that faster algorithms for them imply faster algorithms for canonical hard problems in fine-grained complexity.
From Schatten -norm to Triangle Detection:
We start with the fact that the number of triangles in any unweighted graph is equal to , where is the adjacency matrix. Any algorithm for approximating to high enough accuracy therefore gives an algorithm for detecting if a graph has at least one triangle.
is not PSD, so is actually not a function of ’s singular values – it depends on the signs of ’s eigenvalues. However, the graph Laplacian given by where is the diagonal degree matrix, is PSD and we have:
[TABLE]
since has an all [math] diagonal. Further, it is not hard to see that . So this term and are easy to compute exactly. Thus, if we approximate up to additive error , we can determine if or and so detect if contains a triangle. for any unweighted graph on nodes, and hence computing this norm up to relative error for suffices to detect a triangle. If we have an time approximation algorithm for the Schatten -norm, we can thus perform triangle detection in time.
Our strongest algorithmic result for the Schatten -norm requires just time for dense matrices. Improving the dependence to for the current value of , would yield an algorithm for triangle detection running in time for general graphs, breaking a longstanding runtime barrier for this problem. Even a dependence would give a sub-cubic time triangle detection algorithm, and hence could be used to give a subcubic time matrix multiplication algorithm via the reduction of [WW10].
Generalizing to Other Spectral Sums
We can generalize the above approach to the Schatten -norm by adding self-loops to each node of , which corresponds to replacing with . We then consider . This is the sum over all vertices of the number of paths that start at and return to in four steps. All of these paths are either (1) legitimate four cycles, (2) triangles combined with self loops, or (3) combinations of self-loops and two-step paths from a vertex to one of its neighbors and back. The number of type (3) paths is exactly computable using the node degrees and number of self loops. Additionally, if the number of self loops is large enough, the number of type paths will dominate the number of type paths, even if there is just a single triangle in the graph. Hence, an accurate approximation to will give us the number of type paths, from which we can easily compute the number of triangles.
This argument extends to a very broad class of spectral sums by considering a power series expansion of and showing that for large enough , is dominated by along with some exactly computable terms. Thus, an accurate approximation to this spectral sum allows us to determine the number of triangles in . This approach works for any that can be represented as a power series, with reasonably well-behaved coefficients on some interval of , giving bounds for all with , the SVD entropy, , , and .
We further show that approximating for the used in our lower bound can be reduced to computing all effective resistances of a certain graph Laplacian up to error. Thus, we rule out highly accurate (with dependence for small ) approximation algorithms for all effective resistances, despite the existence of linear time system solvers (with error dependence) for Laplacians [ST04]. Effective resistances and leverage scores are quantities that have recently been crucial to achieving algorithmic improvements to fundamental problems like graph sparsification [SS08] and regression [LMP13, CLM*+*15]. While crude multiplicative approximations to the quantities suffice for these problems, more recently computing these quantities has been used to achieve breakthroughs in solving maximum flow and linear programming [LS14], cutting plane methods [LSW15], and sampling random spanning trees [KM09, MST15]. In each of these settings having more accurate estimates would be a natural route to either simplify or possibly improve existing results; we show that this is unlikely to be successful if the precision requirements are two high.
1.5 Paper Outline
Section 2: Preliminaries. We review notations that will be used throughout.
Section 3: Spectral Window Approximation. We show how to approximately restrict the spectrum of a matrix to a small window, which will be our main primitive for accessing the spectrum.
Section 4: Spectral Histogram Approximation. We show how our spectral window algorithms can be used to compute an approximate spectral histogram. We give applications to approximating general spectral sums, including the Schatten -norms, Orlicz norms, and Ky Fan norms.
Section 5: Lower Bounds. We prove lower bounds showing that highly accurate spectral sum algorithms can be used to give algorithms for triangle detection and matrix multiplication.
Section 6: Improved Algorithms via Polynomial Approximation. We demonstrate how to tighten dependencies in our runtimes using a more general polynomial approximation approach.
Section 7: Optimized Runtime Bounds. We instantiate the techniques of Section 6 give our best runtimes for the Schatten -norms and SVD entropy.
2 Preliminaries
Here we outline notation and conventions used throughout the paper.
Matrix Properties: For we assume without loss of generality that . We let denote the matrix’s singular values, denote the number of non-zero entries, and denote the maximum number of non-zero entries in any row. Note that .
Fast Matrix Multiplication: Let denote the current best exponent of fast matrix multiplication [Wil12, GU17]. Additionally, let denote the exponent such that it takes time to multiply a matrix by a matrix for any . for where and for [LG12, GU17]. For , .
Asymptotic Notation: We use notation to hide poly-logarithmic factors in the input parameters, including dimension, failure probability, and error . We use ‘with high probability’ or ‘w.h.p.’ to refer to events happening with probability at least for some constant , where is our smaller input dimension.
Other: We denote . For any and PSD , we denote
3 Approximate Spectral Windows via Ridge Regression
In this section, we give state-of-the-art results for approximating spectral windows over . As discussed, our algorithms will split ’s spectrum into small slices using these window functions, performing trace estimation to estimate the number of singular values on each window and producing an approximate spectral histogram.
In Section 3.1 we show how to efficiently apply smooth approximations to threshold functions of the spectrum provided we have access to an algorithm for solving regularized regression problems with the matrix. In Section 3.2 we then provide the fastest known algorithms for the regression problems in the given parameter regimes using both stochastic gradient methods and traditional solvers. Departing from our high level description in Section 1.3, we actually incorporate singular vector deflation directly into our system solvers to reduce condition number. This simplifies our final algorithms but has the same effect as the deflation techniques discussed in Section 1.3. Finally, in Section 3.3 we provide algorithms and runtime analysis for applying smooth approximations to window functions of the spectrum, which is the main export of this section.
3.1 Step Function Approximation
To compute a window over ’s spectrum, we will combine two threshold functions at the boundaries of the window. We begin by discussing how to compute these threshold functions.
Let be the threshold function at . for and [math] for . For some gap we define a soft step function by:
Definition 2** (Soft Step Function).**
* is a -soft step at if:*
[TABLE]
We use the strategy from [FMMS16], which, for with shows how to efficiently multiply a -soft step by any using ridge regression. The trick is to first approximately compute where . Then, note that . Additionally, the symmetric step function can be well approximated with a low degree polynomial. Specifically, there exists a polynomial of degree that is within additive of a true -soft step at and can be applied stably such that any error in computing remains bounded. The upshot, following from Theorem 7.4 of [AZL16] is:
Lemma 3** (Step Function via Ridge Regression).**
Let be an algorithm that on input , , returns such that with high probability. Then there is an algorithm which on input with , , , and , returns with
[TABLE]
where is a -soft step at (i.e. satisfies Definition 2). requires calls to along with additional runtime, where .
3.2 Ridge Regression
Given Lemma 3, to efficiently compute for satisfying Definition 2, it suffices to quickly approximate (i.e. to provide the algorithm used in the lemma). In this section we provide two theorems which give the state-of-the-art ridge regression running times achievable in our parameter regime, using sampling, acceleration, and singular value deflation.
Naively, computing using an iterative system solver involves a dependence on the condition number . In our theorems, this condition number is replaced by a deflated condition number depending on for some input parameter . We achieve this improved dependence following the techniques presented in [GOSS16]. We first approximate the top singular vectors of and then construct a preconditioner based on this approximation, which significantly flattens the spectrum of the matrix. By using this preconditioner in conjunction with a stochastic gradient based linear system solver, we further enjoy an average condition number dependence. The following theorem summarizes the results.
Theorem 4** (Ridge Regression – Accelerated Preconditioned SVRG).**
For any and , let . Let where is an input parameter. There is an algorithm that builds a preconditioner for using precomputation time for sparse or time for dense , and for any input , returns such that with high probability in
[TABLE]
time for sparse or time for dense .
Proof.
We give a proof in Appendix A. Note that the dependence in the runtime is and so is hidden by the notation. ∎
When is dense, the runtime of Theorem 4 is essentially the best known. Due to its average condition number dependence, the method always outperforms traditional iterative methods, like conjugate gradient, up to factors. However, in the sparse case, traditional approaches can give faster runtimes if the rows of are not uniformly sparse and is large. We have the following, also proved in Appendix A using the same deflation-based preconditioner as in Theorem 4:
Theorem 5** (Ridge Regression – Preconditioned Iterative Method).**
For any and , let and where is an input parameter. There is an algorithm that builds a preconditioner for using precomputation time , and for any input , returns such that with high probability in time.
3.3 Overall Runtimes For Spectral Windows
Combined with Lemma 3, the ridge regression routines above let us efficiently compute soft step functions of ’s spectrum. Composing step functions then gives our key computational primitive: the ability to approximate soft window functions that restrict ’s spectrum to a specified range. We first define our notion of soft window functions and then discuss runtimes. The corresponding Theorem 7 is our main tool for spectrum approximation in the remainder of the paper.
Definition 6** (Soft Window Function).**
Given with , and , is a -soft window for the range if:
[TABLE]
Theorem 7** (Spectral Windowing).**
For any with , , and , with , there is an algorithm that returns satisfying w.h.p.:
[TABLE]
where is a soft window function satisfying Def. 6. Let and where is an input parameter. The algorithm uses precomputation time for sparse or for dense after which given any it returns in time:
[TABLE]
for sparse or for dense .
Proof.
If then we can simply define for any satisfying Definition 2. Otherwise, given soft steps and satisfying Definition 2, we can define . Since we can verify that this will be a valid soft window function for (i.e. satisfy Definition 6). Further, we have for any :
[TABLE]
We can compute and each up to error via Lemma 3. This gives the error bound in the theorem, since we have both and so the computation in (4) does not amplify error. Our runtime follows from combining Theorems 4 and 5 with with Lemma 3. The errors in these theorems are measured with respect to . To obtain the error in as used by Lemma 3, we simply apply the theorems with which incurs an additional cost. Since the runtime is dominated by the computation of , which depends on the condition number when using SVRG (Theorem 4) or for a traditional iterative solver (Theorem 5). ∎
4 Approximating Spectral Sums via Spectral Windows
We now use the window functions discussed in Section 3 to compute an approximate spectral histogram of . We give our main histogram algorithm and approximation guarantee in Section 4.1. In Section 4.2 we show how this guarantee translates to accurate spectral sum approximation for any smooth and sufficiently quickly growing function . In Section 4.3 we apply this general result to approximating the Schatten -norms for all real , bounded Orlicz norms, and the Ky Fan norms.
4.1 Approximate Spectral Histogram
Our main histogram approximation method is given as Algorithm 1. The algorithm is reasonably simple. Assuming (this is w.l.o.g. as we can just scale the matrix), and given cutoff , below which we will not evaluate spectrum, we split the range into successive windows where . Here determines the width of our windows. In our final spectral approximation algorithms, we will set . is a random shift, which insures that, in expectation, the boundaries of our soft windows do not overlap too many singular values. This argument requires that most of the range is not covered by boundary regions. Thus, we set the steepness parameter where will control the error introduced by the boundaries. Finally, we iterate over each window, applying trace estimation to approximate the singular value count in each window.
In our final algorithms, the number of windows and samples required for trace estimation will be . The dominant runtime cost will come from computing the window for the lowest range , which will incur a dependence on the condition number of with .
Theorem 8** (Histogram Approximation).**
Let , be output by Algorithm 1. Let , for , and be the number of squared singular values of on the range . Then, for sufficiently small constants , with probability , for all , output by Algorithm 1 satisfies:
[TABLE]
That is, approximates the number of singular values of the range up to multiplicative error and additive error . Note that by setting , the error on each bucket is just multiplicative on its size plus the size of the two adjacent buckets, which contain singular values in nearby ranges. For simplicity we assume passed to the algorithm has . This is without loss of generality: we can estimate in time via the power or Lanczos methods [KW92, MM15], and scale down the matrix appropriately.
Ignoring logarithmic and dependencies, the runtime of Algorithm 1 is dominated by the calls to for the bucket corresponding to the smallest singular values, with . This runtime is given by Theorem 7. Since balancing the deflation parameter in that theorem with the minimum squared singular value considered can be complex, we wait to instantiate full runtimes until employing Algorithm 1 for specific spectral sum computations.
Proof.
We use the notation of Algorithm 1, where , is chosen uniformly in and . With this notation, we have . Let be a -soft window for (Definition 6) and let be the interval on which is nonzero. is an estimation of the trace of such a window applied to . We first show that, if these traces are computed exactly, they accurately estimate the singular value counts in each range . We have:
[TABLE]
We can bound the second term using the random shift . Since , each singular value falls within at most two extended ranges: and for some . Let be the set of indices whose singular values fall within two ranges. We have only if for some . Letting , this holds only if , which occurs with probability at most since is chosen uniformly in the range . Thus we have . Further, is distributed uniformly in the range , so we know for certain that if the constant on is set small enough:
[TABLE]
Let . Note that is fixed (i.e. not a random variable). Regardless of the random shift , by (4.1), we always have . We also have .
Let be if and [math] otherwise. We can upper bound the second term of (4.1) by:
[TABLE]
The first bound follows from the fact that for , we have and that . The second bound follows from (4.1), which shows that . Let be the number of squared singular values falling in . Note that like , is fixed (i.e., not a random variable.) Thus, by linearity of expectation, we have:
[TABLE]
Letting as in Algorithm 1, by a Markov bound, with probability :
[TABLE]
if is set small enough. By a union bound this holds for all simultaneously with probability . Plugging back into (4.1), we have, with probability , simultaneously for all :
[TABLE]
where the second bound follows from the fact that, regardless of the setting of the shift , so . Plugging (4.1) into (4.1) we have:
[TABLE]
We conclude by showing that, before the final rounding step of Algorithm 1, with high probability. In the final rounding step, if (which can only occur if ) and are sufficiently small, then we will have and so will round down to . Otherwise, the term will be absorbed into the relative error on .
Overall, combined with (9) we will have, with probability for all :
[TABLE]
which gives the theorem if we adjust by making sufficiently small.Thus we conclude by showing that in fact with high probability. Setting as in Algorithm 1, for chosen from , with high probability by a standard trace estimation result [AT11].
Further, let . By Theorem 7: which by Cauchy-Schwarz gives since . Thus, overall we have, before the rounding step in which is set to [math] if :
[TABLE]
∎
4.2 Application to General Spectral Sums
While Theorem 8 is useful in its own right, we now apply it to approximate a broad class of spectral sums. We need two assumptions on the sums that we approximate. First, for the histogram discretization to be relatively accurate, we need our function to be relatively smooth. Second, it is expensive to compute the histogram over very small singular values of (i.e. with very small in Algorithm 1) as this makes the condition number in Theorem 7 large. So it is important that small singular values cannot contribute significantly to our sum. We start with the following definition:
Definition 9** (Multiplicative Smoothness).**
* is -multiplicatively smooth if for some , for all , .*
We have the following claim, proven in Appendix D.
Claim 10**.**
Let be a -multiplicatively smooth function. For all and
[TABLE]
For the example of the Schatten -norm, for , and so is -multiplicatively smooth.
We now give our general approximation theorem, showing that any spectral sum depending on sufficiently smooth and rapidly growing can be computed using Algorithm 1:
Theorem 11** (Spectral Sums Via Approximate Histogram).**
Consider any and any function satisfying:
- •
Multiplicative Smoothness: For some , is -multiplicatively smooth (Defn. 9).
- •
Small Tail: For any there exists such that for ,
Given error parameter and spectral norm estimate , for sufficiently small constant , if we run Algorithm 1 on with input parameters , and then with probability , letting , be the outputs of the algorithm and :
[TABLE]
For parameter , letting and , the algorithm runs in
[TABLE]
time for sparse or for dense .
That is, we accurately approximate by discretizing over the histogram output by Algorithm 1. Note that we can boost our probability of success to by repeating the algorithm times and taking the median of the outputs.
Proof.
Let . Note that and so it is a valid matrix on which to apply Algorithm 1 and Theorem 8. Recall that we use the notation: , for , , and .
Since , by chain rule:
[TABLE]
So is multiplicatively smooth. By this smoothness, and Claim 10, for any with :
[TABLE]
if the constant on is set small enough. Small enough also ensures and . So, both the multiplicative and additive error terms in Theorem 8 do not hurt us significantly. For now, assume that . We will first prove the result with this assumption and then show that it can be relaxed to as given in the theorem statement. Let be the set of indices with for some . That is, the set of singular values covered by our histogram. Let be the set of indices with – that is, singular values which are not included in the histogram, but may be included in the additive error for the last bucket covering . Applying Theorem 8:
[TABLE]
if we set small enough. The second inequality arises because each contributing to appears at most three times as an additive error term for , , and . In this inequality we include bucket in the histogram, which only increases the right hand side.
On the lower bound side, we use our small tail assumption, that for with , . Using the notation of Algorithm 1, we have . So for any with , and thus falls in some bucket of our histogram so . We thus have:
[TABLE]
Applying Theorem 8 again, we have for sufficiently small :
[TABLE]
We conclude the theorem by noting that we can actually set instead of as used above. The additive error term on each bucket was bounded using a Markov bound, and to union bound over buckets, we lost a factor of . However, in expectation, the total contribution of additive error to our spectral sum estimation is just and so by a Markov bound is with probability if we set our constants small enough.
It just remains to discuss runtime, i.e. to calculate the runtime of Algorithm 1 with inputs , and . The number of outer loops of the algorithm is . The number of inner loops is . And finally, within each loop the most expensive operation is computing . Our final runtime follows from plugging this into Theorem 7 noting that and for all . Note that we perform the precomputation step to construct a preconditioner for just once, incurring cost or . ∎
4.3 Application to Schatten , Orlicz, and Ky Fan Norms
Theorem 11 is very general, allowing us to approximate any function satisfying a simple smoothness condition as long as the smaller singular values of cannot contribute significantly to . We now give some specific applications of the result.
Schatten -Norms
Theorem 11 already gives the fastest known algorithms for Schatten -norm estimation. We will not go into all runtime tradeoffs now as our best runtimes will be worked out in detail in Sections 6 and 7, but as an example:
Corollary 12** (Schatten -norms via Histogram Approximation).**
For any with uniformly sparse rows (i.e. ), given error parameter and , if we run Algorithm 1 on with , and for sufficiently small constant then with probability , letting be the outputs of the algorithm we have:
[TABLE]
Further the algorithm runs in time:
[TABLE]
For dense inputs this can be sped up to using fast matrix multiplication or without fast matrix multiplication.
For constant the first runtime is , and for the nuclear norm (), for constant the second runtime gives which is at worst . For dense matrices, the nuclear norm estimation time is using fast matrix multiplication. It is already , for the current , without using fast matrix multiplication.
Note that we can compute the spectral norm approximation used to scale via the Lanczos or power method in time. depends on which we are estimating. However, as we will discuss in the proof, we can use a rough estimate for which suffices. As a more general strategy, can be identified via binary search. We can start with and successively decrease running Algorithm 1 up to the stated runtime bounds. If it does not finish in the allotted time, we know that we have set too small. Thus, we can output the result with the smallest such that the algorithm completes within in the stated bounds.
Proof.
We invoke Theorem 11 with . We have so and our setting of suffices. Additionally, for any , we can set and so our setting of suffices. Thus the accuracy bound follows from Theorem 11.
We now consider runtime. For :
[TABLE]
We can bound and so have For we have:
[TABLE]
Using the fact that we have the tradeoff between and :
[TABLE]
As mentioned, depends on the value of . We can simply lower bound by , which we estimate up to multiplicative error when performing deflation. We can use this lower bound to set . Our estimated will only be smaller than the true value, giving a better approximation guarantee and the above condition number bound will still hold.
We now analyze our runtime. Recall that for , . Correspondingly, and so . Plugging into the first runtime of Theorem 11, using the uniform sparsity assumption and the fact that we have:
[TABLE]
For we just set and have runtime by our bound . For , not trying to optimize terms, we write the runtime as
[TABLE]
and balancing the first two coefficients on , set which gives by (11) and so . We then have . Finally, the is dominated by the term so we drop it.
Finally, for dense we apply the third runtime which gives
[TABLE]
We now balance the terms, again not optimizing dependence. Writing , for where and for [GU17]. Assuming we set:
[TABLE]
which gives for all (so our assumption that was valid.) This yields total runtime . Without using fast matrix multiplication, the first term in the runtime becomes and so we balance costs by setting: which gives and total runtime .
∎
Bounded Orlicz Norms
An Orlicz norm of a vector is given by , where is convex, nonnegative, and has . It is easy to see that applied to the vector of singular values, an Orlicz matrix norm is a special case of a spectral sum and can be approximated with Theorem 11 under sufficient conditions. A simple example shows that Orlicz norms for bounded by an envelope of and can be approximated in a similar time as in Corollary 12.
Corollary 13** (Bounded Orlicz norms via Histogram Approximation).**
For any convex function , let be the minimal real for which for all arguments to with , . Let be the maximum positive real for which for all we have .
Then for any with uniformly sparse rows (i.e. ), given error parameter and , if we run Algorithm 1 on with , and for sufficiently small constant then with probability , letting be the outputs of the algorithm we have:
[TABLE]
*Further the algorithm runs in time: *
[TABLE]
For dense inputs this can be sped up to via fast matrix multiplication.
Note that depends on which can be estimated for example using Corollary 12, or via binary search as described in the proof of Corollary 12.
Proof.
We invoke Theorem 11 with . While we do not show that exactly satisfies the multiplicative smoothness condition, the upper bound on directly implies the result of Claim 10 with . So our setting of suffices as the bound in this Claim is the only smoothness condition used to prove Theorem 11.
Additionally, by the lower bound on , for any , we can set . For any :
[TABLE]
which gives:
[TABLE]
where the last bound follows from the fact that by the setting of , Finally for any , . So overall we have:
[TABLE]
Thus our setting of suffices and the accuracy bound follows from Theorem 11.
It remains to discuss runtime. We have and . The runtimes follow from Theorem 11 via the same arguments used in Corollary 12. ∎
Ky Fan Norms
The Ky Fan -norm of a matrix is the sum of its top singular values: (note that these norms are typically called the ‘Ky Fan -norms’, however we use to avoid overloading notation on ). Such a norm is not strictly speaking a spectral sum. However it can still be approximated using our histogram method. We have the following corollary of Theorem 8:
Corollary 14** (Ky Fan norms via Histogram Approximation).**
For any , given rank , error parameter and , if we run Algorithm 1 on with , , and for sufficiently small constants then with probability , letting be the outputs of the algorithm we have:
[TABLE]
where is the smallest integer with or if no such integer exists (and hence the term of the sum is just [math]). Further the algorithm runs in time:
[TABLE]
Note that to compute the top singular values explicitly would require time using for example a block Krylov method [MM15]. Also note that depends on the norm we are attempting to compute. As discussed in Corollary 12, we can approximate via binary search, successively decreasing it until our call to Algorithm 1 exceeds the stated runtime.
Proof.
We can assume without loss of generality that and since rescaling will not effect our approximation factor. By Theorem 8, since we set , , and have :
[TABLE]
if we set small enough compared to . This give if is set small enough. Thus, since these scaled bucket sizes strictly overestimate the true bucket sizes we have , where is the smallest integer with . This gives, since and since by our setting of , at most an fraction of falls outside of the ranges :
[TABLE]
On the other side, let be a vector that, for each has entries each set to . Define the top- norm of to be the sum of its largest entries, denoted by . Our estimate of is equal to . Note that we can also add arbitrary zeros as padding entries to and not change this norm.
Similarly, let be a vector with values each set to for and entries set to .
By (12), . This fact combined with the fact that the entries in are scaled up by a factor ensure that the entries of dominate those of and so . Further, since we have scaled up each entry by at most factor and scaled up each bucket size by a factor. This bound gives our final multiplicative approximation after adjusting the constant on .
It remains to discuss runtime. We invoke Theorem 7, setting . This gives:
[TABLE]
where we us the fact that for , . With our settings of , and , Algorithm 1 performs calls to , each which requires regression calls due to our setting of and . Plugging into the second runtime of Theorem 7 gives the corollary. ∎
5 Lower Bounds
In this section we give hardness results for high accuracy spectrum approximation. Specifically, we show how to detect if an undirected unweighted graph contains a triangle using accurate approximation algorithms for various important spectral sums such as the Schatten -norms, log-determinant, the SVD entropy, and the trace inverse. Our spectral sum bounds further imply hardness for the important primitives of computing effective resistances in a graph or leverage scores in a matrix.
In the seminal work of [WW10] it was shown that any truly subcubic time algorithm for triangle detection yields a truly subcubic time algorithms for Boolean matrix multiplication (BMM). Consequently, these results show that computing any of these quantities too precisely is in a sense as difficult as BMM. Furthermore, as it is a longstanding open question whether or not there is any subcubic time combinatorial algorithm for BMM, i.e. an algorithm which avoids the powerful algebraic manipulations inherent in the fastest algorithms for BMM, these results can be viewed as showing that we do not expect simple iterative methods to yield precise algorithms for these problems without an immense breakthrough in linear algebra.
In Section 5.1 we give a general result on reducing spectral sums to triangle detection. Then in Section 5.2 we use this to show hardness for computing various well studied spectral sums. In Section 5.3 we conclude by showing hardness for computing graph and numerical linear algebra primitives, i.e. effective resistances and leverage scores.
5.1 Reductions From Triangle Detection
Here we provide our main technical tool for reducing spectral sum computation to triangle detection. As discussed in Section 1.4, our reduction leverages the well known fact that the number of triangles in any unweighted graph is equal to where is the adjacency matrix for . Consequently, given any function whose power series is reasonably behaved, we can show that for suitably small the quantity is dominated by the contribution of for . Therefore computing approximately lets us distinguish between whether or not or .
We formalize this in the following theorem. As it simplifies the result, we focus on the case where is defined on the interval , however, this suffices for our purposes and can be generalized via shifting and scaling of .
Theorem 15** (Many Spectral Sums are as Hard as Triangle Detection).**
Let be an arbitrary function, such that for we can express it as
[TABLE]
Then given the adjacency matrix of any simple graph that has no self-loops and spectral sum estimate for scaling and accuracy satisfying
[TABLE]
we can detect if has a triangle in time.
Consequently, given an algorithm which on input outputs in time we can produce an time triangle detection algorithm.
Proof.
Let , , , , and be as in the theorem statement and let . By Gershgorin’s circle theorem, and since , . Consequently is symmetric PSD, for all , and therefore applying (13) yields:
[TABLE]
is enough to insure that the first three terms of this power series dominate. Specifically:
[TABLE]
where the last inequality uses the fact that for all . Further, since , , and we have:
[TABLE]
Now, clearly in time we can compute in as well as:
[TABLE]
So we can detect if or if and hence whether or not has a triangle. ∎
Note that in this reduction, so long as is small (i.e. ) then is a very well conditioned matrix (its condition number is at most a constant). Consequently, our lower bounds apply even when approximately applying for example or to high precision is inexpensive. The theorem (and the results in Section 5.3) suggests that the difficulty in computing spectral sums arises more from the need to measure the contribution from multiple terms precisely, than from the difficulty in manipulating for the purposes of applying it to a single vector.
Also, note that the matrix in this reduction is symmetrically diagonally dominant (SDD). So, even for these highly structured matrices which admit near linear time application of [ST04] as well as approximate factorization [KS16], accurate spectral sums are difficult. We leverage this in Section 5.3.
5.2 Hardness for Computing Spectral Sums
Here we use Theorem 15 to show hardness for various spectral sum problems. To simplify our presentation, we focus on the case of dense matrices, showing bounds of the form . However, note that Theorem 15 and the approach we use also yields lower bounds on the running time for sparse matrices and can be stated in terms of .
Schatten -norm for all
For , using the Taylor Series about we can write
[TABLE]
This series converges since for and for , . Note that when is a non-negative integer, only the first terms of the expansion are nonzero. When is non-integral, the sum is infinite. We will apply Theorem 15 slightly differently for different values of . We first give our strongest result:
Corollary 16** (Schatten -Norm Hardness).**
Given an algorithm which on input returns in we can produce an algorithm that detects if an -node graph contains a triangle in time.
Proof.
For , for . So we apply Theorem 15 with and hence and . ∎
Note that for very close to a similar bound holds as . If Theorem 31 gives an algorithm running in time. Significant improvement to the dependence in this algorithm therefore either requires loss in the dependence or would lead to time triangle detection for the current value of . We next extend to all .
Corollary 17** (Schatten -Norm Hardness, ).**
For any , , given algorithm which for any returns in time, we can detect if an -node graph contains a triangle in time.
Proof.
We have for all as well as and similarly . We apply Theorem 15 with and , which gives the result. ∎
In the typical case when , the term above is negligible. The term is meaningful however. Our bound becomes weak as approaches (and meaningless when ). This is unsurprising, as for very close to , , which can be computed exactly in time. The bound also becomes weak for , which is natural as our reduction only uses PSD , for which which can be computed in time. However, we can remedy this issue by working with a (non-PSD) square root of which is easy to compute:
Corollary 18** (Schatten -Norm Hardness, ).**
For any , given an algorithm which for any returns in time, we can detect if an -node graph with edges contains a triangle in time.
Note that for , is just a constant. Again, the bound is naturally weak when as goes to [math]
Proof.
For as in Theorem 15. Let be the Laplacian of where is the diagonal degree matrix. We can write where is PSD since . Letting be the vertex edge incidence matrix of , and , we have . Thus, and so approximating this norm gives triangle detection by Corollary 17. Note that and further matrix can be formed in this amount of time, giving our final runtime claim. ∎
Note that for , since has maximum row sparsity , we obtain a runtime via Theorem 33 of for the current value of , even when , implying that significantly improving this dependence would either improve matrix multiplication or come at a cost in the polynomials of the other parameters.
SVD Entropy:
Corollary 19** (SVD Entropy Hardness).**
Given algorithm which for any returns for in time, we can detect if an -node graph contains a triangle in time.
Proof.
For , using the Taylor Series about we can write where , , and for . for all and while . Applying Theorem 15 with and gives the result. ∎
Log Determinant:
Corollary 20** (Log Determinant Hardness).**
Given algorithm which for any returns in time, we can detect if an -node graph contains a triangle in time.
Proof.
Using the Taylor Series about we can write where , for . Therefore for all , , and . Applying Theorem 15 with and gives the result. ∎
In Appendix B, Lemma 52 we show that a similar result holds for computing . In [BS83] it is shown that, given an arithmetic circuit for computing , one can generate a circuit of the same size (up to a constant) that computes . This also yields a circuit for matrix multiplication by a classic reduction.333Matrix multiplication reduces to inversion by the fact that
. See [Isa08]. Our results, combined with the reduction of [WW10] of Boolean matrix multiplication to triangle detection, show that a sub-cubic time algorithm for the approximating or up to sufficient accuracy, yields a sub-cubic time matrix multiplication algorithm, providing a reduction based connection between determinant and matrix multiplication analogous to the circuit based result of [BS83].
Trace of Exponential:
Corollary 21** (Trace of Exponential Hardness).**
Given algorithm which for any returns in time, we can detect if an -node graph contains a triangle in time.
Proof.
Using the Taylor Series about we can write . We have , , and for all , . Applying Theorem 15 with and gives the result. ∎
5.3 Leverage Score and Effective Resistance Hardness
Here we show hardness for precisely computing all effective resistances and leverage scores of a graph. Our main tool is the following result (which in turn is an easy corollary of Theorem 15) for an algorithm that precisely computes the trace inverse of a strictly symmetric diagonally dominant (SDD) , i.e. and .
Corollary 22** (Trace of Inverse Hardness).**
Given an algorithm which for any strictly SDD with non-positive off-diagonal entries returns in time for , we can produce an algorithm which detects if an -node graph contains a triangle in time.
Proof.
For we can write , and then apply Theorem 15 with and . Checking that the in Theorem 15 is strictly SDD with non-positive off-diagonal entries yields the result. ∎
Using this we prove hardness for precisely computing effective resistances in a graph. Recall that for a weighted undirected graph its Laplacian, is given by if there is an edge between and and [math] otherwise and or equivalently where is the diagonal degree matrix and is the weighted adjacency matrix associated with . Note that this describes a clear bijection between a Laplacian and its associated graph and we therefore use them fairly interchangeably in the remainder of this section.
The effective resistance between vertices and is given by where denotes the Moore-Penrose pseudoinverse. In the following lemma we prove that compute all the effective resistances between a vertex and its neighbors in the graph can be used to compute the trace of the inverse of any strictly SDD matrix with non-positive off-diagonals and therefore doing this precisely is as hard as triangle detection via Corollary 22. Our proof is based off a fairly standard reduction between solving strictly SDD matrices with negative off-diagonals and solving Laplacian systems.
Lemma 23** (Effective Resistance Yields Trace).**
Suppose we have an algorithm which given Laplacian with -non-zero entries, entry , and error computes a approximation to the total effective resistance between and the neighbors of in the graph associated with , that is
[TABLE]
in time . Then there is an algorithm that computes the trace of the inverse of strict SDD matrix with non-zero entries and non-positive off-diagonals in time.
Proof.
Let be an arbitrary strictly SDD matrix with non-positive off-diagonals, i.e. , , and for all . Let , , and
[TABLE]
Now, clearly by our assumptions on we have that entrywise and therefore . Therefore, the off-diagonal entries of are non-positive and by construction . Consequently, is a symmetric Laplacian matrix with non-zero entries.
Now, consider any and that satisfy the following for some
[TABLE]
Since is a symmetric Laplacian and the associated graph is connected by construction we know that and there there must exist such and . Furthermore, since is strictly SDD it is invertible and since we have that
[TABLE]
and consequently
[TABLE]
Consequently, if we used the algorithm to get a multiplicative approximation as stated then and the result follows. ∎
Using this, we also show that computing leverage scores of matrix, a powerful and prevalent notion in convex analysis and numerical linear algebra, is also difficult to compute. This follows from the well known fact that effective resistances in graphs and leverage scores of matrices are the same up to scaling by known quantities.
Corollary 24** (Leverage Score Hardness).**
Suppose we have an algorithm which given can compute that is a multiplicative approximation to the leverage scores of , i.e.
[TABLE]
in time . Then there is a time algorithm for detecting if an -node graph contains a triangle.
Proof.
Let be a symmetric Laplacian. Let , i.e. the set of edges in the graph associated with . Let and be the incidence matrix associated with , i.e. for all we have and for some canonical choice of ordering of and and let all other entries of . Clearly and we can form in time.
It is well known and easy to check that . Consequently, for all we have
[TABLE]
Now if we compute using the assumed algorithm in time, then since is non-negative in an additional time this yields a multiplicative approximation to the total effective resistance between any and all its neighbors in the graph associated with . Consequently, by Lemma 23 and Corollary 22 the result follows. ∎
6 Improved Error Dependence via Polynomial Approximation
For constant , Theorem 11 and the resulting Corollary 12 matches our fastest runtimes for the Schatten -norms. However, it is possible to significantly improve the dependence by generalizing our approach. Instead of splitting our spectrum into many small spectral windows, we split into windows of constant multiplicative width and approximate via a low degree polynomial over each window. The degree of this polynomial only varies with .
In Theorem 11, fixing to be constant for illustration, each window has width and so there are windows. Additionally, we must set the steepness parameter to be . This loses us a total factor in our runtime as compared to our improved algorithm which will set and so .
We begin by showing that if we well approximate on each window, then we well approximate overall. In the following theorems for simplicity we work with PSD matrices, as we will eventually apply these lemmas to .
Lemma 25**.**
Consider , parameters , and gap parameter . Set , pick uniformly at random from , and set for . Let for all . Let be a -soft window for (Definition 6) and let be the interval on which is nonzero. Furthermore, for , let be a ‘well-behaved’ approximation to on in the following sense:
- •
Multiplicative Approximation: for .
- •
Approximate Monotonicity: for all for some .
- •
Range Preserving: .
Then, for any PSD with , with probability :
[TABLE]
Proof.
Note that the restriction that along with the approximate monotonicity property implies that is nonnegative, so is a valid spectral sum and , so we will be able to estimate this sum via stochastic trace estimation.
Let denote the unique index such that . Since , each lies in the support of at most two overlapping soft windows (in at most two ranges and ). Let be the set of indices whose singular values fall in the support of exactly one soft window and be its complement. We first bound the error introduced on singular values with indices in .
[TABLE]
The last inequality follows from the multiplicative approximation assumption on that for . Let us now consider a particular and bound . Note that there are precisely two non-zero terms in this summation – one corresponding to and the other, to , which we denote by . Using the above hypothesis on again, we see that
[TABLE]
For the term involving , we have and by the approximate monotonicity requirement that have:
[TABLE]
Due to the random positioning of the windows, in expectation, the total contribution of the singular values lying at the intersection of two windows is small. Specifically,
[TABLE]
where the bound on follows from the fact that only if for some . This holds only if , which occurs with probability since is chosen uniformly in the range . By a Markov bound applied to (18), with probability we have . Combining with (15), (16), and(17), we obtain, with probability 9/10:
[TABLE]
∎
We now show that, as long as the contribution of the smaller singular values of to is not too large, we can truncate this sum and still obtain an accurate approximation. Specifically:
Corollary 26**.**
For any PSD with , let be a function such that, for any there exists such that for , . Given parameters , and gap parameter , for define , and, as in Lemma 25. Let be a well-behaved approximation to on as in Lemma 25. With probability :
[TABLE]
Proof.
This follows from Lemma 25 along with the small tail assumption. Specifically, by (17):
[TABLE]
∎
We now show that well behaved polynomial approximations exist for the function for general real , whose spectral sums give the Schatten -norms.
Lemma 27** (Polynomial Approximation of Power Function).**
For all and let
[TABLE]
Then for all and we have and
[TABLE]
Proof.
Induction shows that the -th derivative of at , is given by
[TABLE]
Furthermore, direct calculation reveals that for all and
[TABLE]
Consequently, is the degree Taylor expansion of about evaluated at . Furthermore, since , the integral form of the remainder of a Taylor series expansion shows
[TABLE]
In the last step we took the worst case of in the integral and used that since it is the case that and therefore . Consequently, and
[TABLE]
∎
Corollary 28**.**
For any , with , there is a polynomial of degree such that:
- •
Multiplicative Approximation: for .
- •
Monotonicity: for all .
- •
Range Preserving: .
Proof.
Set and for large enough constant . Instantiate Lemma 27 with and to obtain . Set . It is clear that . Further,
[TABLE]
For , . Further, and so if we choose large enough we have . We finally show monotonicity. We have . All are positive, so this difference is monotonically decreasing in . We thus have, for :
[TABLE]
since . This gives us monotonicity. ∎
We now combine the approximations of Corollary 28 with Corollary 26 to give our improved Schatten norm estimation algorithm.
Lemma 29** (Schatten Norms via Polynomial Approximation).**
For any with , and error parameter , for sufficiently small , let , , and . For define , , and as in Lemma 25. Let be as defined in Corollary 28 with , for sufficiently small and . With probability :
[TABLE]
Proof.
We apply Corollary 26 to , with , , and . satisfies the necessary conditions by Corollary 28. Further, we have . Plugging in , into (19) we have:
[TABLE]
which gives approximation if we set small enough. ∎
Theorem 30** (Schatten Norm Polynomial Approximation Algorithm).**
For any , and there is an algorithm returning that runs in time:
[TABLE]
time for sparse or
[TABLE]
for dense .
Proof.
We apply Lemma 29, first scaling so that . We apply a trace estimator to approximate up to multiplicative error plus very small additive error for each of our windows. The trace estimation requires applications of , since the degree of as given in Corollary 28 is .
The cost to apply to a vector is given by Theorem 7. Following the argument in Corollary 12, if we write for some ,
[TABLE]
Similarly, . We can further optimize our dependence by noting that it suffices to approximate up to multiplicative error since even if there are singular values below , they will contribute at most a fraction of . So our total cost of approximating the trace for each of our windows using the first runtime of Theorem 7 is:
[TABLE]
Factoring out the term, we have:
[TABLE]
Note that if , the terms in the geometric sum are increasing so the largest term is and so the whole thing is dominated by . If , then the largest term in the geometric sum is as so similarly the whole thing is .
A similar argument gives the runtimes for standard iterative and dense methods. Adding in the cost for deflation-based preconditioning gives the final runtimes. Note that while depends on which is unknown, it again suffices to lower bound the truncation point using . This lower bound gives only a better approximation and our bounds on and still hold. ∎
7 Optimized Runtime Results
The runtimes given in Theorem 30 are quite complex, with many tradeoffs depending on the properties of . In this section, we instantiate the theorem showing how to optimize and achieve our final runtime bounds. For simplicity, we consider the case when .
7.1 Schatten -Norms
We first tackle the relatively simply case when is dense.
Theorem 31** (Schatten -Norms for Dense Matrices).**
For any and there is an algorithm returning which runs in time for and
[TABLE]
time for . If we do not use fast matrix multiplication the runtime is .
Note that for the important case of our runtime is or if we do not use fast matrix multiplication. As approaches [math], our runtime approaches which is natural, as gives the matrix rank, which seems difficult to determine. As approaches the runtime smoothly approaches , which is then required for all .
Proof.
The bound for follows immediately from Theorem 30 with set to [math]. For we have runtime:
[TABLE]
We can partially optimize this term, ignoring and factors for simplicity and setting to equalize the coefficients on . Our optimization is identical to the argument in Corollary 12 (set when not using fast matrix multiplication), giving the stated runtimes. ∎
We now tackle the more complex case when is sparse. We first consider .
Theorem 32** (Schatten -norms, , for Sparse Matrices).**
For any , and there is an algorithm returning with high probability in time:
[TABLE]
Note that if our matrix is uniformly sparse, and so the first runtime becomes . The second runtime can be better when is not uniformly sparse.
Proof.
The first runtime follows by setting in Theorem 30. Note that . For the second, we use the second runtime of Theorem 30 which gives
[TABLE]
Setting to balance the coefficients on the terms gives which yields the result. Note that the term is dominated by the cost of the regressions, even if we do not use fast matrix multiplication this term will be . ∎
Finally, we consider the most complex cost, for sparse matrices. We have:
Theorem 33** (Schatten -norms, , for Sparse Matrices).**
For any , and there is an algorithm returning with high probability in time:
[TABLE]
where .
Note that in the special case of the first runtime gives . The second term here is at worst , and could be much smaller for sparse .
Proof.
For the second runtime, we use the second runtime of Theorem 30, balancing costs exactly as in Theorem 32 (setting ). For the first runtime, we consider the first runtime of Theorem 30 which gives:
[TABLE]
Ignoring , , and dependence, one way we can balance the costs is by setting:
[TABLE]
which gives and final runtime:
[TABLE]
Note that we drop the term as it is dominated by the last term with in it. This gives our final result by noting that , the average row sparsity. ∎
7.2 Constant Factor Approximation without Row-Sparsity Dependence
Theorem 34** (Removing the Uniform Sparsity Constraint).**
Let , , and . There is an algorithm returning with with high probability in time for , and in time for .
Proof.
For , let be the matrix whose rows agree with that of provided the number of non-zero entries in these rows are in the interval , while the remaining rows of are set to . Let . Since we can apply the triangle inequality,
[TABLE]
For any matrix obtained from a matrix by replacing some rows of with and preserving the remaining rows, . This follows from the fact that for all vectors , together with the min-max theorem for singular values. Hence,
[TABLE]
Combining (22) and (23), an -approximation to is an -approximation to .
We remove the zero rows from the , obtaining corresponding matrices . By definition of , the number of rows in is at most . For each , which is an matrix for , we can right-multiply it by an OSNAP matrix , with columns and non-zero entries per column, so that for all with probability , see [Coh16] (the fact that follows from the fact that is a subspace embedding of the row space of together with Lemma C.2 of [LNW14]). The time to compute is . Since has non-zero entries per row, each row of has a number of non-zero entries within an factor of the corresponding row of .
We compute -approximations to the for each , and then take the maximum over . To do so, for we run the algorithm of Theorem 32 on each , while for we run the algorithm of Theorem 33.
Let denote the parameter when run on . Then . We can assume is a square matrix by padding with zeros, where . Hence, . Also, .
The total time to apply Theorem 32 for across all is . The total time to apply Theorem 33 for across all is , using that the of that theorem is . ∎
7.3 SVD Entropy
In this section, we will show how to approximately estimate the SVD entropy of a matrix assuming its condition number to be bounded by . Recall that the SVD entropy of a matrix is given by , where is the matrix scaled inversely by its Schatten -norm. Using the results in Section 7.1, we obtain a approximation to i.e., . Define . This means that
[TABLE]
The following thoerem gives our main result for SVD entropy.
Theorem 35**.**
Given any PSD and error parameter . Then, there exists an algorithm that computes the approximate SVD entropy of and outputs , such that, with probability
[TABLE]
where . The runtime of the algorithm will be
[TABLE]
where , and .
The algorithm we consider to estimate the SVD entropy of follows the techniques of [HNO08]. The algorithm depends on the SVD Tsallis entropy of matrix , which is given by
[TABLE]
We can describe the multi-point interpolation method for SVD entropy estimation with additive and multiplicative error approximations using [HNO08, Algorithm 1], see Algorithm 2.
Here, given an error parameter , we compute an approximate SVD Tsallis entropy at different points , where . are defined as follows. Let and let be defined as:
[TABLE]
Now define an error parameter .
Algorithm 2 gives the stepwise algorithm. The runtime of the algorithm is dominated by the cost of computing , a -approximation for Schatten norm , particularly the smallest power Schatten norm. We get the smallest power when . We have such Schatten norms to be estimated. Thus, we obtain the runtime in Theorem 35 by using in the runtime for Schatten norm estimation given in Corollary 12, and observing that is between 1 and 2 (we set ).
Additive approximation
: Section 3.3.2 of [HNO08] tells us that if in Algorithm 2, the approximate is an additive approximation to for every , then we can use these to compute an additive- approximation to . Since , and since , obtaining multiplicative approximation to suffices. This approximation can be obtained using the results of Section 7.1.
Multiplicative approximation
: Article [HNO08] further extends the multi-point interpolation method to achieve multiplicative approximation for Shannon entropy (equivalently for SVD entropy) using the following modifications. We set the number of interpolation points . Then, section 3.4 in [HNO08] shows that if (in Algorithm 2) for every computed are to be -multiplicative approximation to (where is as defined above), then we can achieve multiplicative approximation for using the multi-point interpolation method with these . So, we need to obtain -relative error approximations to the at each .
If the matrix has reasonable (large) SVD entropy (meaning the singular values are uniform and there is no one singular value that is very large), then we can obtain approximation to the by simply computing approximation to using results from section 7.1. However, if the matrix has very small entropy, i.e., we have one singular value with very large magnitude and remaining singular values are very small, then achieving a multiplicative approximation will be difficult. This is because, we are approximating of , which is matrix scaled by a approximation of its Schatten -norm. The approximation obtained for might not be close to in this case.
This issue is equivalent to the ‘heavy element’ issue (one of the entries in the vector is very large) discussed in [HNO08]. Hence, to overcome the above issue, we can follow the workaround proposed in [HNO08]. Specifically, Lemma 5.5 in [HNO08] shows that a -approximation to together with a -approximation to , where are the singular values of with , suffices to get a -approximation to at . It follows that we just need to obtain -approximations for these latter two quantities for each .
It can be shown that when the Krylov method [MM15] (or power method) is used to deflate the top singular vector, we have that any unitarily invariant norm of the tail (remaining part of the spectrum) is preserved, see Appendix C for the proof. This means we can get -approximations to both and , where is matrix with the top singular vector exactly deflated off, by running the algorithm presented in this paper. That is, get the Schatten and Schatten norms for the matrix with the top singular vector deflated off. We can also approximate up to relative error using our algorithm. Then, the two quantities above ( and ) are exactly and , respectively. Since we have relative -approximations to the numerators and denominators of both these quantities, we obtain then up to -relative error the quantities and , as needed to achieve a multiplicative approximation to the SVD entropy. Note that we do not need to compute exactly here even when the matrix has very low SVD entropy. A similar argument can be seen in section 6 of [HNO08].
Acknowledgements
We thank Vladimir Braverman for pointing out an error in our original proof of Theorem 8, which has been corrected.
Appendix A Linear System Solvers
In this section we give runtimes for solving ridge regression using both traditional iterative methods and stochastic gradient descent equipped with deflation-based preconditioning. We start with a few standard notions from convex optimization, which are necessary for our stochastic method bounds.
Definition 36** (Strong convexity).**
A function is -strongly convex if, for all ,
[TABLE]
Equivalently, if .
Definition 37** (Smoothness).**
A function is -smooth if for all ,
[TABLE]
Equivalently, if .
The rate of convergence that iterative methods achieve on typically depends on the ratio , which corresponds to the condition number of a linear system. We next show how ridge regression can be recast as minimizing a convex function , and show that our error from the true ridge solution is proportional to error in minimizing .
Fact 38** (Ridge Regression as Convex Optimization).**
For any , , and , let and . is the minimizer of the convex function:
[TABLE]
which has gradient . is -strongly convex and -smooth.
Proof.
To check that is the minimizer, notice that:
[TABLE]
Since the function is quadratic, is the unique minimizer. , so by Definitions 36 and 37, is -strongly convex and -smooth. ∎
Fact 39** (Function Error Equals Norm Error).**
For any , letting , and be defined as in Fact 38,
[TABLE]
Proof.
Since :
[TABLE]
∎
We will focus on making multiplicative progress in which will lead to multiplicative progress in and a close approximation in iterations.
A.1 Unaccelerated SVRG
We first prove an unaccelerated and unpreconditioned runtime for the Stochastic Variance Reduced Gradient (SVRG) algorithm, introduced in [JZ13].
Theorem 40** (Standard SVRG Runtime).**
For any , , and , let and . There is an algorithm that returns with: in
[TABLE]
time.
SVRG proceeds in epochs. In each epoch we will make one full gradient computation – where is the value of our iterate at the beginning of the epoch. We will then make a number of stochastic gradient updates each ‘variance reduced’ using and show that we make constant factor progress on our function value in expectation. Stringing together epochs yields Theorem 40. We write our function given by (26) as a sum:
[TABLE]
where is the row of and is the rank- matrix which is its contribution to . We start with a well known lemma showing that it is possible to make constant progress on the value of in each epoch of SVRG:
Lemma 41** (SVRG Epoch).**
Consider a set of convex functions mapping . Let and . Suppose we have a probability distribution on and that starting from some initial point in each iteration we select with probability and set:
[TABLE]
for some step size . If is -strongly convex and if for all we have
[TABLE]
for some variance parameter then for all we have:
[TABLE]
Consequently, if we pick to be a small multiple of , then for we will decrease the error by a constant multiplicative factor in expectation.
Proof.
See for example Theorem 9 of [GHJ*+*16]. ∎
To apply the above Lemma we need the variance bound:
Lemma 42** (SVRG Variance Bound).**
If and then letting we have:
[TABLE]
Proof.
We have Write for simplicity of notation.
[TABLE]
where the last step uses Fact 39. ∎
We can now plug this variance bound into Lemma 41 to prove Theorem 40
Proof of Theorem 40.
Using Lemma 42 we can instantiate Lemma 41 with . If we set
[TABLE]
and (which we can compute explicitly) then after an step SVRG epoch:
[TABLE]
If we choose uniformly from this gives us . So stringing together of these epochs and letting in the first epoch gives the theorem. Each epoch requires stochastic gradient steps given by (29), taking time plus time to compute , giving us the stated runtime.
Note that, naively, to produce we need time, not time since our gradient step includes adding multiples of and both of which might be dense vectors in . However, in each epoch, we will just keep track of the coefficients of these two components in our iterate, allowing us to still compute (29) in time. ∎
A.2 Unaccelerated SVRG with Deflation-Based Preconditioning
If we assume that (i.e. our rows are uniformly sparse) we see that the runtime of Theorem 40 is dominated by . can be seen as an average condition number, which is always smaller than the condition number . This average condition number dependence means that SVRG can significantly outperform traditional iterative solvers that require a number of iterations depending on . However, this advantage is limited if is very concentrated in a few large singular values, and hence . We can perform better in such situations by deflating off these large singular values and preconditioning with our deflated matrix, significantly ‘flattening’ the spectrum of . This method was used in [GOSS16], and we follow their approach closely, giving our own proof for completeness and so that we can express runtimes in terms of the necessary parameters for our downstream results. In particular, we show that the preconditioning methods of [GOSS16] can be implemented efficiently for sparse systems.
Theorem 43** (Preconditioned SVRG Runtime).**
For any , , and , let . Let where is an input parameter. There is an algorithm that uses precomputation time for sparse matrices or for dense matrices after which, given any input , letting the algorithm returns with in
[TABLE]
time for sparse or
[TABLE]
time for dense .
Our ideal algorithm would compute the top singular vectors of and deflate these off our matrix to flatten the spectrum. However, for efficiency we will instead compute approximate singular vectors, using an iterative method, like simultaneous iteration or a block Krylov iteration. These algorithms give the following guarantee [MM15]:
Lemma 44**.**
There is an algorithm that, with high probability in time returns such that and, for all , letting , . For dense inputs the runtime can be sped up to by applying fast matrix multiplication at each iteration.444Each iteration requires time since we assume , which dominates the time required to orthogonalize the approximate singular directions.
[GOSS16] shows that we can build a good preconditioner from such a . Specifically:
Lemma 45** (Theorem 5 of [GOSS16]).**
For any , , and , given satisfying the guarantees of Lemma 44, let
[TABLE]
Then we have:
[TABLE]
Intuitively, after applying the preconditioner, all top singular values are capped at , giving a much flatter spectrum and better performance when optimizing with SVRG. To make use of the above bound, we first define a preconditioned ridge regression problem.
Definition 46** (Preconditioned Ridge Regression).**
For any , , and , let and . Letting be as described in Lemma 45, let and define the preconditioned regression objective function by:
[TABLE]
* is minimized at .*
Fact 47** (Preconditioned Solution).**
For any , if then where is the minimizer of .
Proof.
We have and so using the fact that can write:
[TABLE]
∎
By Fact 47, we can find a near optimal solution to our preconditioned system and by multiplying by obtain a near optimal solution to the original ridge regression problem. We now show that our preconditioned problem can be solved quickly. We first give a variance bound as we did in Lemma 42 in the non-preconditioned case.
Lemma 48** (Preconditioned Variance Bound).**
Let so that . For , let
[TABLE]
so that . Set then letting :
[TABLE]
Proof.
Write for simplicity of notation. Then we have:
[TABLE]
where the last step uses Fact 39. Finally we write:
[TABLE]
∎
We can now combine this variance bound with Lemmas 41 and 45 to prove our preconditioned SVRG runtime.
Proof of Theorem 43.
Using Lemma 48 we can instantiate Lemma 41 with . Letting be a lower bound on the strong convexity of , by Lemma 45 we have:
[TABLE]
If we set the step size which we will compute explicitly, then after an step SVRG epoch running on the functions each selected with probability , we make constant progress in expectation on . After iterations we have with and hence, by Fact 47. It remains to bound the cost of each SVRG epoch applied to .
For dense inputs, the argument is relatively simple. We can compute the approximate top singular vector space in time by Lemma 44. Using the factored structure of , we can explicitly form and also in time. This allows us to compute a full gradient in time (which we do once per epoch) and perform each stochastic gradient step in time (which we do times per epoch). This gives final runtime
[TABLE]
where .
For sparse we have to be much more careful to fully exploit sparsity when preconditioning. Computing takes time by Lemma 44. We will not explicitly form but will show how to apply it when needed. Recall that where we denote . First, in order to determine our step size we must compute .
[TABLE]
Applying Pythagorean theorem and submultiplicativity we have:
[TABLE]
which can be computed in time with just a single multiplication of by .
We can similarly compute our sampling probabilities quickly by noting that for , and for , . This again requires just time.
In each SVRG epoch we must compute one full gradient of the form:
[TABLE]
which takes time as can be applied to a vector in time. Naively, we then must make stochastic gradient steps, each requiring time. However we can do better. Instead of storing our iterate explicitly, we will store it as the sum:
[TABLE]
Note that , and is a scalar. At the beginning of each epoch we will set , , and . We can compute in time. We will maintain the invariant that is perpendicular to the span of while will give the component of within this span.
For ease of notation let denote where the standard basis vector. In this way we have for all . Each stochastic gradient step is of the form:
[TABLE]
We can precompute in time per epoch and then can compute the dot product in time. Using the fact that is always perpendicular to the span of , we can write the dot product as:
[TABLE]
We can precompute in time. We can also precompute and in time. With these values in hand, computing the above dot product takes just time.
Now, we can write . Since we have precomputed , can be computed and added to with the appropriate weight in time. We can then compute for the appropriate weight :
[TABLE]
which takes time . Finally, we set . Overall, our runtime per epoch is and so our total runtime is:
[TABLE]
where ∎
A.3 Accelerated and Preconditioned SVRG
We now combine the deflation-based preconditioning described above with accelerated gradient methods to give our strongest runtime bound for ridge regression using via stochastic solvers.
Theorem 49** (Accelerated Preconditioned SVRG Runtime).**
For any , , and , let and where is an input parameter. There is an algorithm that uses precomputation time for sparse matrices or for dense matrices after which, given any input , letting the algorithm returns with in
[TABLE]
time for sparse or
[TABLE]
time for dense .
The above runtime will follow from applying a blackbox technique for accelerating the runtime of convex optimization methods. While there are a number of improvements over this method [AZ17], we at most lose logarithmic factors and gain significant simplicity in our proofs.
Lemma 50** (Theorem 1.1 of [FGKS15]).**
Let be -strongly convex and let . For any and any , let . Let . Suppose that, for all , , , we can compute such that
[TABLE]
in time . Then we can compute such that in time
[TABLE]
Proof of Theorem 49.
We focus again on optimizing the preconditioned function in Definition 46 as by Fact 47 a near optimal minimizer for this function yields a near optimal solution for our original ridge regression problem. We split as where
[TABLE]
is as defined in Lemma 48 with .
We have . Letting , we can follow a similar calculation to Lemma 42 to show:
[TABLE]
Let and be an lower bound on the strong convexity of . Denote . The strong convexity of is lower bounded by which gives:
[TABLE]
So, by Theorem 43, for dense inputs, ignoring the precomputation cost to compute , and , which we only pay once, letting , we have . If then we already have runtime by Theorem 43. Otherwise, setting we can solve up to accuracy in time:
[TABLE]
For sparse inputs, we again ignore the one time precomputation cost of to compute and . We have:
[TABLE]
If then we already have runtime by Theorem 43. Otherwise, . So setting we have:
[TABLE]
∎
A.4 Preconditioned Iterative Methods
Finally, we describe how to combine deflation-based preconditioning with standard iterative methods, which can give runtime advantages over Theorem 49 in some parameter regimes.
Theorem 51** (Preconditioned Iterative Method).**
For any , , and , let and where is an input parameter. There is an algorithm that uses precomputation time after which, given any input , letting the algorithm returns such that with high probability in time
[TABLE]
Proof.
If we form the preconditioner as in Lemma 45, by Lemmas 2 and 4 of [GOSS16] we have the preconditioned condition number bound:
[TABLE]
Note that the bound in [GOSS16] is actually in terms of , however we write so that our theorem holds in the case that . This can be achieved with no effect on the asymptotic runtime by setting to .
We can now apply any accelerated linear system solver, such as Chebyshev iteration, Conjugate Gradient, or Accelerated Gradient Descent [Saa03, Nes13] to obtain with in where is the time required to multiply a single vector by this matrix. since can be applied in time. This gives the result combined with the precomputation time for from Lemma 44. ∎
Appendix B Additional Proofs: Lower Bounds
Lemma 52** (Determinant Hardness).**
Given algorithm which given returns in time, we can detect if an -node graph contains a triangle in time.
Proof.
Let be the adjacency matrix of an -node graph . Let denote its eigenvalues. Let for some which we will set later. We can write:
[TABLE]
The term in (31) is , and the next two are easy to compute. , and For we have:
[TABLE]
We will bound the terms by: since . However, in order to obtain a tighter result, we will use stronger bounds for . These bounds are very tedious but straightforward. Specifically:
[TABLE]
And similarly:
[TABLE]
Finally, if we set then we have:
[TABLE]
We then write:
[TABLE]
Since and if we compute for sufficiently small constant and subtract off , we will be able to determine if and hence detect if has a triangle. So any algorithm approximating to error in time yields a triangle detection algorithm running in time. ∎
Appendix C Krylov and Power Methods preserve Invariant Norms
In this section, we show that when the Krylov method [MM15] or the power method is used to deflate the top singular vectors of the matrix, any unitarily invariant norm of the tail (remaining part of the spectrum) is preserved. Let be the projector obtained for the top singular vectors of using the Krylov method [MM15] or the power method. Then, for any invariant norm of the tail to be preserved, we just need to show the following:
[TABLE]
Since is a approximation obtained in terms of the spectral norm , we have by the spectral low rank approximation guarantee of [MM15] that , where is the best rank approximation of . Consider first any for which . Then we have,
[TABLE]
For any with , i.e., with there is a large () relative gap between this singular value and . We have by the min-max characterization of singular values:
[TABLE]
If we just set , then we have:
[TABLE]
It follows that the bottom singular values are all [math], and so equal to the bottom singular values of .
Now, for , let denote the top singular vectors of and denote the bottom singular vectors. We set . First, we note that has columns and also . This is because, if we consider , the top left blocks is and the bottom right block is . The off-diagonal entries in the top right and bottom left blocks are all bounded by ; the proof of this latter statement is given in the latter part of this section, where we bound . By the Gershgorin circle theorem [GVL12], all eigenvalues of are in the range and so the matrix is full rank. So, we have:
[TABLE]
Next, we can write for some . By our argument above, every singular value of lies in . Then, we have . Splitting , where contains the first coordinates of the vector and contains the rest, we have:
[TABLE]
Then, we have
[TABLE]
where the inequality comes from bounding using the fact that all its entries are at most (see the latter part of this section), and bounding . Thus, finally we obtain:
[TABLE]
Plugging into (33) gives the proof of (32). Note that we can assume since if , we can just compute the SVD.
Bound for :
Let be an orthonormal basis for from the power method (or from block Krylov method). Then, for any singular vector with corresponding singular value with , we have .
To see this, we first write in its SVD form. Then, we have , where is the th diagonal entry of . On the other hand, we also have , where is the right singular vector corresponding to of . Since is i.i.d. Gaussian, we have with probability for a fixed (see proof of Lemma 2.12 in [Woo14]). So, we can union bound over all right singular vectors and the relation holds for all right singular vectors . We condition on this event.
This event implies that . It follows that, if we show , then.
Next, to show that , it suffices to show . Indeed, note that for (for the block Krylov method, we need to use higher number of iterations than required, which is [MM15]), by definition of , we have .
For this, we have since for any projection matrix and any matrix . Then, we have by the definition of the SVD, and that the columns of are orthonormal. So, we get , where is a matrix with i.i.d. Gaussian distribution due to the rotational invariance of the Gaussian distribution. Since the -th singular value is the smallest singular value of and , we then have with probability at least , for an arbitrary constant , using standard properties of minimum singular values of squared Gaussian matrices (see e.g., Fact 6 in section 4.3 of [Woo14]). Thus, we have , which completes the proof.
Appendix D Additional Proofs
Claim 10.
Let be a -multiplicatively smooth function on the range . For any and
[TABLE]
Proof.
Let denote the range . For we have:
[TABLE]
Similarly, letting we have:
[TABLE]
and hence . So plugging into (34) we have:
[TABLE]
which gives the result. ∎
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[ABB 00] Orly Alter, Patrick O Brown, and David Botstein. Singular value decomposition for genome-wide expression data processing and modeling. Proceedings of the National Academy of Sciences , 97(18):10101–10106, 2000.
- 2[AKR 15] Alexandr Andoni, Robert Krauthgamer, and Ilya P. Razenshteyn. Sketching and embedding are equivalent for norms. In Proceedings of the \nth 47 Annual ACM Symposium on Theory of Computing (STOC) , pages 479–488, 2015.
- 3[AT 11] Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. Journal of the ACM , 58(2):8, 2011.
- 4[AZ 17] Zeyuan Allen-Zhu. Katyusha: The first direct acceleration of stochastic gradient methods. Proceedings of the \nth 49 Annual ACM Symposium on Theory of Computing (STOC) , 2017.
- 5[AZL 16] Zeyuan Allen-Zhu and Yuanzhi Li. Faster principal component regression via optimal polynomial approximation to sgn(x). ar Xiv:1608.04773 , 2016.
- 6[BCF 09] Constantine Bekas, Alessandro Curioni, and I Fedulova. Low cost high performance uncertainty quantification. In Proceedings of the 2nd Workshop on High Performance Computational Finance , page 8. ACM, 2009.
- 7[BCKY 16] Vladimir Braverman, Stephen R Chestnut, Robert Krauthgamer, and Lin F Yang. Sketches for matrix norms: Faster, smaller and more general. ar Xiv:1609.05885 , 2016.
- 8[BDKZ 15] Christos Boutsidis, Petros Drineas, Prabhanjan Kambadur, and Anastasios Zouzias. A randomized algorithm for approximating the log determinant of a symmetric positive definite matrix. ar Xiv:1503.00374 , 2015.