Revisiting the Nystrom Method for Improved Large-Scale Machine Learning

Introduction

inline]Use a remark environment for numbering and better formatting of the remarks inline]Talwalkar: References and are often mixed up during citations, e.g., on page 2 (last paragraph), page 5 (above equation 5 and below equation 6). In particular, provides coherence-based bounds for the Nystrom method in the low-rank setting.

We reconsider randomized algorithms for the low-rank approximation of symmetric positive semi-definite (SPSD) matrices such as Laplacian and kernel matrices that arise in data analysis and machine learning applications. Our goal is to obtain an improved understanding, both empirically and theoretically, of the complementary strengths of sampling versus projection methods on realistic data. Our main results consist of an empirical evaluation of the performance quality and running time of sampling and projection methods on a diverse suite of dense and sparse SPSD matrices drawn both from machine learning as well as more general data analysis applications. These results are not intended to be comprehensive but instead to be illustrative of how randomized algorithms for the low-rank approximation of SPSD matrices behave in a broad range of realistic machine learning and data analysis applications.

In addition to being of interest in their own right, our empirical results point to several directions that are not explained well by existing theory. (For example, that the results are much better than existing worst-case theory would suggest, and that sampling with respect to the statistical leverage scores leads to results that are complementary to those achieved by projection-based methods.) Thus, we complement our empirical results with a suite of worst-case theoretical bounds for both random sampling and random projection methods. These bounds are qualitatively superior to existing bounds—e.g., improved additive-error bounds for spectral and Frobenius norm error and relative-error bounds for trace norm error. Importantly, by considering random sampling and random projection algorithms on an equal footing, we identify within our analysis deterministic structural properties of the input data and sampling/projection methods that are responsible for high-quality low-rank approximation.

In more detail, our main contributions are fourfold.

First, we provide an empirical illustration of the complementary strengths and weaknesses of data-independent random projection methods and data-dependent random sampling methods when applied to SPSD matrices. We do so for a diverse class of SPSD matrices drawn from machine learning and more general data analysis applications, and we consider reconstruction error with respect to the spectral, Frobenius, as well as trace norms. Depending on the parameter settings, the matrix norm of interest, the data set under consideration, etc., one or the other method might be preferable. In addition, we illustrate how these empirical properties can often be understood in terms of the structural nonuniformities of the input data that are of independent interest.

Second, we consider the running time of high-quality sampling and projection algorithms. For random sampling algorithms, the computational bottleneck is typically the exact or approximate computation of the importance sampling distribution with respect to which one samples; and for random projection methods, the computational bottleneck is often the implementation of the random projection. By exploiting and extending recent work on “fast” random projections and related recent work on “fast” approximation of the statistical leverage scores, we illustrate that high-quality leverage-based random sampling and high-quality random projection algorithms have comparable running times. Although both are slower than simple (and in general much lower-quality) uniform sampling, both can be implemented more quickly than a naïve computation of an orthogonal basis for the top part of the spectrum.

Third, our main technical contribution is a set of deterministic structural results that hold for any “sketching matrix” applied to an SPSD matrix. (A precise statement of these results is given in Theorems 1, 2, and 3 in Section 4.1.) We call these “deterministic structural results” since there is no randomness involved in their statement or analysis and since they depend on structural properties of the input data matrix and the way the sketching matrix interacts with the input data. In particular, they highlight the importance of the statistical leverage scores (and other related structural nonuniformities having to do with the subspace structure of the input matrix), which have proven important in other applications of random sampling and random projection algorithms.

Fourth, our main algorithmic contribution is to show that when the low-rank sketching matrix represents certain random projection or random sampling operations, then we obtain worst-case quality-of-approximation bounds that hold with high probability. (A precise statement of these results is given in Lemmas 2, 3, 4, and 5 in Section 4.2.) These bounds are qualitatively better than existing bounds (when nontrivial prior bounds even exist); they hold for reconstruction error of the input data with respect to the spectral norm and trace norm as well as the Frobenius norm; and they illustrate how high-quality random sampling algorithms and high-quality random projection algorithms can be treated from a unified perspective.

A novel aspect of our work is that we adopt a unified approach to these low-rank approximation questions—unified in the sense that we consider both sampling and projection algorithms on an equal footing, and that we illustrate how the structural nonuniformities responsible for high-quality low-rank approximation in worst-case analysis also have important empirical consequences in a diverse class of SPSD matrices. By identifying deterministic structural conditions responsible for high-quality low-rank approximation of SPSD matrices, we highlight complementary aspects of sampling and projection methods; and by illustrating the empirical consequences of structural nonuniformities, we provide theory that is a much closer guide to practice than has been provided by prior work. More generally, we should note that, although it is beyond the scope of this paper, our deterministic structural results could be used to check, in an a posteriori manner, the quality of a sketching method for which one cannot establish an a priori bound.

Our analysis is timely for several reasons. First, in spite of the empirical successes of Nyström-based and other randomized low-rank methods, existing theory for the Nyström method is quite modest. For example, existing worst-case bounds such as those of are very weak, especially compared with existing bounds for least-squares regression and general low-rank matrix approximation problems .This statement may at first surprise the reader, since an SPSD matrix is an example of a general matrix, and one might suppose that the existing theory for general matrices could be applied to SPSD matrices. While this is true, these existing methods for general matrices do not in general respect the symmetry or positive semi-definiteness of the input. Moreover, many other worst-case bounds make very strong assumptions about the coherence properties of the input data . Second, there have been conflicting views in the literature about the usefulness of uniform sampling versus nonuniform sampling based on the empirical statistical leverage scores of the data in realistic data analysis and machine learning applications. For example, some work has concluded that the statistical leverage scores of realistic data matrices are fairly uniform, meaning that the coherence is small and thus uniform sampling is appropriate ; while other work has demonstrated that leverage scores are often very nonuniform in ways that render uniform sampling inappropriate and that can be essential to highlight properties of downstream interest . inline]Talwalkar: In we don’t strictly claim that uniform sampling is the best. A more accurate statement is the quote at the end of Section 4: "The empirical results suggest a trade-off between time and space requirements, as noted by Scholkopf and Smola (2002)[Chapter 10.2]. Adaptive techniques spend more time to find a concise subset of informative columns, but as in the case of the K-means algorithm, can provide improved approximation accuracy." Third, in recent years several high-quality numerical implementations of randomized matrix algorithms for least-squares and low-rank approximation problems have been developed . These have been developed from a “scientific computing” perspective, where condition numbers, spectral norms, etc. are of greater interest , and where relatively strong homogeneity assumptions can be made about the input data. In many “data analytics” applications, the questions one asks are very different, and the input data are much less well-structured. Thus, we expect that some of our results will help guide the development of algorithms and implementations that are more appropriate for large-scale analytics applications.

In the next section, Section 2, we start by presenting some notation, preliminaries, and related prior work. Then, in Section 3 we present our main empirical results; and in Section 4 we present our main theoretical results. We conclude in Section 5 with a brief discussion of our results in a broader context.

Notation, Preliminaries, and Related Prior Work

In this section, we introduce the notation used throughout the paper, and we address several preliminary considerations, including reviewing related prior work.

Given $\mathbf{A}$ and a rank parameter $k$, the statistical leverage scores of $\mathbf{A}$ relative to the best rank-$k$ approximation to $\mathbf{A}$ equal the squared Euclidean norms of the rows of the $n\times k$ matrix $\mathbf{U}_{1}$:

denote the projection of $\mathbf{S}$ onto the top and bottom eigenspaces of $\mathbf{A}$, respectively.

inline]Address Ilse’s comments on the numerical analysts’ usual use of the terminology relative-error, and her concern about not knowing how our TCS idea of relative-error is calibrated (i.e. state that close to 1 is good, and why, and that we want to achieve close to 1 with as few col samples as possible). So our desidera are low relative-error achieved for few col samples

2. Preliminaries

In many machine learning and data analysis applications, one is interested in symmetric positive semi-definite (SPSD) matrices, e.g., kernel matrices and Laplacian matrices. One common column-sampling-based approach to low-rank approximation of SPSD matrices is the so-called Nyström method . The Nyström method— both randomized and deterministic variants—has proven useful in applications where the kernel matrices are reasonably well-approximated by low-rank matrices; and it has been applied to Gaussian process regression, spectral clustering and image segmentation, manifold learning, and a range of other common machine learning tasks . The simplest Nyström-based procedure selects columns from the original data set uniformly at random and then uses those columns to construct a low-rank SPSD approximation. Although this procedure can be effective in practice for certain input matrices, two extensions (both of which are more expensive) can substantially improve the performance, e.g., lead to lower reconstruction error for a fixed number of column samples, both in theory and in practice. The first extension is to sample columns with a judiciously-chosen nonuniform importance sampling distribution; and the second extension is to randomly mix (or combine linearly) columns before sampling them. For the random sampling algorithms, an important question is what importance sampling distribution should be used to construct the sample; while for the random projection algorithms, an important question is how to implement the random projections. In either case, appropriate consideration should be paid to questions such as whether the data are sparse or dense, how the eigenvalue spectrum decays, the nonuniformity properties of eigenvectors, e.g., as quantified by the statistical leverage scores, whether one is interested in reconstructing the matrix or performing a downstream machine learning task, and so on.

The following sketching model subsumes both of these classes of methods.

The choice of distribution for the sketching matrix $\mathbf{S}$ leads to different classes of low-rank approximations. For example, if $\mathbf{S}$ represents the process of column sampling, either uniformly or according to a nonuniform importance sampling distribution, then we refer to the resulting approximation as a Nyström extension; if $\mathbf{S}$ consists of random linear combinations of most or all of the columns of $\mathbf{A}$, then we refer to the resulting approximation as a projection-based SPSD approximation. In this paper, we focus on Nyström extensions and projection-based SPSD approximations that fit the above SPSD Sketching Model. In particular, we do not consider adaptive schemes, which iteratively select columns to progressively decrease the approximation error. While these methods often perform well in practice , rigorous analyses of them are hard to come by—interested readers are referred to the discussion in .

inline]Talwalkar: You should perhaps also also include column-projection approximations (defined in ) in your discussion / experiments. I would argue that a truly unified approach would include Nystrom, Column-Projection and Random Projection (and this indeed is what we do in our DFC work )

inline]Talwalkar: Considering parallel run times would be quite interesting. For instance, random projection can be trivially parallelized, and in a large-scale setting, reporting parallel runtime is appropriate.

3. The Power Method

One can obtain the optimal rank-$k$ approximation to $\mathbf{A}$ by forming an SPSD sketch where the sketching matrix $\mathbf{S}$ is an orthonormal basis for the range of $\mathbf{A}_{k},$ because with such a choice,

SPSD sketches produced using $q$ iterations of the power method have lower error than sketches produced without using the power method, but are roughly $q$ times more costly to produce. Thus, the power method is most applicable when $\mathbf{A}$ is such that one can compute the product $\mathbf{A}^{q}\mathbf{S}_{0}$ fast. We consider the empirical performance of sketches produced using the power method in Section 3, and we consider the theoretical performance in Section 4.

4. Related Prior Work

Motivated by large-scale data analysis and machine learning applications, recent theoretical and empirical work has focused on “sketching” methods such as random sampling and random projection algorithms. A large part of the recent body of this work on randomized matrix algorithms has been summarized in the recent monograph of Mahoney and the recent review article of Halko, Martinsson, and Tropp . Here, we note that, on the empirical side, both random projection methods (e.g., and ) and random sampling methods (e.g., ) have been used in applications for clustering and classification of general data matrices; and that some of this work has highlighted the importance of the statistical leverage scores that we use in this paper . In parallel, so-called Nyström-based methods have also been used in machine learning applications. Originally used by Williams and Seeger to solve regression and classification problems involving Gaussian processes when the SPSD matrix $\mathbf{A}$ is well-approximated by a low-rank matrix , the Nyström extension has been used in a large body of subsequent work. For example, applications of the Nyström method to large-scale machine learning problems include and , and applications in statistics and signal processing include .

Much of this work has focused on new proposals for selecting columns (e.g., ) and/or coupling the method with downstream applications (e.g., ). The most detailed results are provided by (as well as the conference papers on which it is based ). Interestingly, they observe that uniform sampling performs quite well, suggesting that in the data they considered the leverage scores are quite uniform, which also motivated the related work . This is in contrast with applications in genetics , term-document analysis , and astronomy , where the statistical leverage scores were seen to be very nonuniform in ways of interest to the downstream scientist; we return to this issue in Section 3.

On the theoretical side, much of the work has followed that of Drineas and Mahoney , who provided the first rigorous bounds for the Nyström extension of a general SPSD matrix. They show that when $\Omega(k\epsilon^{-4}\ln\delta^{-1})$ columns are sampled with an importance sampling distribution that is proportional to the square of the diagonal entries of $\mathbf{A}$, then

with probability exceeding $1-\delta$ and

We have described these prior theoretical bounds in detail to emphasize how strong, relative to the prior work, our new bounds are. For example, Equation (4) provides an additive-error approximation with a very large scale; the bounds of Kumar, Mohri, and Talwalkar require a sampling complexity that depends on the coherence of the input matrix , which means that unless the coherence is very low one needs to sample essentially all the rows and columns in order to reconstruct the matrix; Equation (5) provides a bound where the additive scale depends on $n$; and Equation (6) provides a spectral norm bound where the scale of the additional error is the (much larger) trace norm. Table 1 compares the bounds on the approximation errors of SPSD sketches derived in this work to those available in the literature. We note further that Wang and Zhang recently established lower-bounds on the worst-case relative spectral and trace norm errors of uniform Nyström extensions . Our Lemma 5 provides matching upper bounds, showing the optimality of these estimates.

A related stream of research concerns projection-based low-rank approximations of general (i.e., non-SPSD) matrices . Such approximations are formed by first constructing an approximate basis for the top left invariant subspace of $\mathbf{A},$ and then restricting $\mathbf{A}$ to this space. Algorithmically, one constructs $\mathbf{Y}=\mathbf{A}\mathbf{S},$ where $\mathbf{S}$ is a sketching matrix, then takes $\mathbf{Q}$ to be a basis obtained from the QR decomposition of $\mathbf{Y},$ and then forms the low-rank approximation $\mathbf{Q}\mathbf{Q}^{T}\mathbf{A}.$ The survey paper proposes two schemes for the approximation of SPSD matrices that fit within this paradigm: $\mathbf{Q}(\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})\mathbf{Q}^{T}$ and $(\mathbf{A}\mathbf{Q})(\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})^{\dagger}(\mathbf{Q}^{T}\mathbf{A}).$ The first scheme—for which provides quite sharp error bounds when $\mathbf{S}$ is a matrix of i.i.d. standard Gaussian random variables—has the salutary property of being numerically stable. On the other hand, although does not provide any theoretical guarantees for the second scheme, it points out that this latter scheme produces noticeably more accurate approximations in practice. In Section 3, we provide empirical evidence of the superior performance of the second scheme, and we show that it is actually an instantiation of the power method (as described in Section 2.3) with $q=2.$ Accordingly, the deterministic and stochastic error bounds provided in Section 4 are applicable to this SPSD sketch.

It is worth noting that in , the authors propose a modified Nyström method wherein the matrix $\mathbf{W}$ is replaced by $\mathbf{C}^{\dagger}\mathbf{A}(\mathbf{C}^{\dagger})^{T},$ so that the low rank approximation to $\mathbf{A}$ is given by $\mathbf{C}\mathbf{C}^{\dagger}\mathbf{A}(\mathbf{C}^{\dagger})^{T}\mathbf{C}^{T}.$ Note that $\mathbf{C}\mathbf{C}^{\dagger}$ is another expression for the orthoprojector $\mathbf{Q}\mathbf{Q}^{T}$ onto the range of $\mathbf{Y}=\mathbf{A}\mathbf{S},$ so this Nyström method is an instantiation of the projection-based low-rank approximations analyzed in . However, , unlike , considers the case where $\mathbf{C}$ is constructed by sampling from the columns of $\mathbf{A}$ adaptively. The low-rank approximation produced by the algorithm proposed in satisfies

5. An overview of our bounds

Our bounds in Table 1 (established as Lemmas 2–5 in Section 4.2) exhibit a common structure: for the spectral and Frobenius norms, we see that the additional error is on a larger scale than the optimal error, and the trace norm bounds all guarantee relative error approximations. This follows from the fact, as detailed in Section 4.1, that low-rank approximations that conform to the SPSD sketching model can be understood as forming column-sample/projection-based approximations to the square root of $\mathbf{A}$, and thus squaring this approximation yields the resulting approximation to $\mathbf{A}.$ The squaring process unavoidably results in potentially large additional errors in the case of the spectral and Frobenius norms— whether or not the additional errors are large in practice depends upon the properties of the matrix and the form of stochasticity used in the sampling process. For instance, from our bounds it is clear that Gaussian-based SPSD sketches are expected to have lower additional error in the spectral norm than any of the other sketches considered.

From Table 1, we also see, in the case of uniform Nyström extensions, a necessary dependence on the coherence of the input matrix since columns are sampled uniformly at random. However, we also see that the scales of the additional error of the Frobenius and trace norm bounds are substantially improved over those in prior results. The large additional error in the spectral norm error bound is necessary in the worse case . Lemmas 2, 3 and 4 in Section 4.2—which respectively address leverage-based, Fourier-based, and Gaussian-based SPSD sketches—show that spectral norm additive-error bounds with additional error on a substantially smaller scale can be obtained if one first mixes the columns before sampling from $\mathbf{A}$ or one samples from a judicious nonuniform distribution over the columns.

inline]Make sure the footnote associated with this table is on the right page after the document is finalized

Several trends can be identified; among them, we note that the bounds provided in this paper for Gaussian-based sketches come quite close to capturing the errors seen in practice, and the Frobenius and trace norm error guarantees of the leverage-based and Fourier-based sketches tend to more closely reflect the empirical behavior than the error guarantees provided in prior work for Nyström sketches. Overall, the trace norm error bounds are quite accurate. On the other hand, prior bounds are sometimes more informative in the case of the spectral norm (with the notable exception of the Gaussian sketches). Several important points can be gleaned from these observations. First, the accuracy of the Gaussian error bounds suggests that the main theoretical contribution of this work, the deterministic structural results given as Theorems 1 through 3, captures the underlying behavior of the SPSD sketching process. This supports our belief that this work provides a foundation for truly informative error bounds. Given that this is the case, it is clear that the analysis of the stochastic elements of the SPSD sketching process is much sharper in the Gaussian case than in the leverage-score, Fourier, and uniform Nyström cases. We expect that, at least in the case of leverage and Fourier-based sketches, the stochastic analysis can and will be sharpened to produce error guarantees almost as informative as the ones we have provided for Gaussian-based sketches.

Empirical Aspects of SPSD Low-rank Approximation

inline]Talwalkar: The empirical results are very detailed, which is great, but having a short summary of take-away messages (e.g., which methods perform best for which types of data) might be helpful.

In this section, we present our main empirical results, which consist of evaluating sampling and projection algorithms applied to a diverse set of SPSD matrices. In addition to understanding the relative merits, in terms of both running time and solution quality, of different sampling/projection schemes, we would like to understand the effects of various data preprocessing decisions. The bulk of our empirical evaluation considers two random projection procedures and two random sampling procedures for the sketching matrix $\mathbf{S}$: for random projections, we consider using SRFTs (Subsampled Randomized Fourier Transforms) as well as uniformly sampling from Gaussian mixtures of the columns; and for random sampling, we consider sampling columns uniformly at random as well as sampling columns according to a nonuniform importance sampling distribution that depends on the empirical statistical leverage scores. In the latter case of leverage score-based sampling, we also consider the use of both the (naïve and expensive) exact algorithm as well as a (recently-developed fast) approximation algorithm. Section 3.1 starts with a brief description of the data sets we consider; Section 3.2 describes the details of our SPSD sketching algorithms; and then Section 3.3 briefly describes the effect of various data preprocessing decisions. In Section 3.4, we present our main results on reconstruction quality for the random sampling and random projection methods; and, in Section 3.5, we discuss running time issues, and we present our main results for running time and reconstruction quality for both exact and approximate versions of leverage-based sampling.

We emphasize that we don’t intend these results to be “comprehensive” but instead to be “illustrative” case-studies—that are representative of a much wider range of applications than have been considered previously. In particular, we would like to illustrate the tradeoffs between these methods in different realistic applications in order, e.g., to provide directions for future work. For instance, prima facie, algorithms based on leverage-based column sampling might be expected to be more expensive than those based on uniform column sampling or random projections, but (based on previous work for general matrices ) they might also be expected to deliver lower approximation errors. Similarly, using approximate leverage scores to construct the importance sampling distribution might be expected to perform worse than using exact leverage scores, but this might be acceptable given its computational advantages. In addition to clarifying some of these issues, our empirical evaluation also illustrates ways in which existing theory is insufficient to explain the success of sampling and projection methods. This motivates our improvements to existing theory that we describe in Section 4.

Table 4 provides summary statistics for the data sets used in our empirical evaluation. In order to illustrate the complementary strengths and weaknesses of different sampling versus projection methods in a wide range of realistic applications, we consider four classes of matrices which are commonly encountered in machine learning and data analysis applications: normalized Laplacians of very sparse graphs drawn from “informatics graph” applications; dense matrices corresponding to Linear Kernels from machine learning applications; dense matrices constructed from a Gaussian Radial Basis Function Kernel (RBFK); and sparse RBFK matrices constructed using Gaussian radial basis functions, truncated to be nonzero only for nearest neighbors. Although not exhaustive, this collection of data sets represents a wide range of data sets with very different (sparsity, spectral, leverage score, etc.) properties that have been of interest recently not only in machine learning but in data analysis more generally.

inline] Talwalkar: it would be nice to remind the user why low-rank approximation is useful for each of the 4 categories of data that you look at. In particular, such motivation is important in cases where the matrices have a slowly decaying spectrum (e.g., Laplacian matrices).

To understand better the Laplacian data, recall that, given an undirected graph with weighted adjacency matrix $\mathbf{W}$, its normalized graph Laplacian is

where $\mathbf{D}$ is the diagonal matrix of weighted degrees of the nodes of the graph, i.e., $D_{ii}=\sum_{j\neq i}W_{ij}$. This Laplacian is an SPSD matrix, but note that not all SPSD matrices can be written as the Laplacian of a graph.

measures the similarity (correlation) of $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ in feature space .

When $\kappa$ is the usual Euclidean inner-product, so that

$\mathbf{A}$ is called a Linear Kernel matrix. Gaussian RBFK matrices, defined by

correspond to the similarity measure $\kappa(\mathbf{x},\mathbf{y})=\exp(-\|\mathbf{x}-\mathbf{y}\|^{2}_{2}/\sigma^{2}).$ Here $\sigma$, a nonnegative number, defines the scale of the kernel. Informally, $\sigma$ defines the “size scale” over which pairs of points $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ “see” each other. Typically $\sigma$ is determined by a global cross-validation criterion, as $\mathbf{A}^{\sigma}$ is generated for some specific machine learning task; and, thus, one may have no a priori knowledge of the behavior of the spectrum or leverage scores of $\mathbf{A}^{\sigma}$ as $\sigma$ is varied. Accordingly, we consider Gaussian RBFK matrices with different values of $\sigma$.

Finally, given the same data points, $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$, one can construct sparse Gaussian RBFK matrices

where $[x]^{+}=\max\{0,x\}.$ When $\nu$ is larger than $(d+1)/2,$ this kernel matrix is positive semidefinite . Increasing $\nu$ shrinks the magnitudes of the off-diagonal entries of the matrix toward zero. As the cutoff point $C$ decreases the matrix becomes more sparse; in particular, $C\rightarrow 0$ ensures that $\mathbf{A}^{(\sigma,\nu,C)}\rightarrow\mathbf{I}.$ On the other hand, $C\rightarrow\infty$ ensures that $\mathbf{A}^{(\sigma,\nu,C)}$ approaches the (dense) Gaussian RBFK matrix $\mathbf{A}^{\sigma}.$ For simplicity, in our empirical evaluations, we fix $\nu=\lceil(d+1)/2\rceil$ and $C=3\sigma$, and we vary $\sigma$. As with the effect of varying $\sigma$, the effect of varying the sparsity parameter $C$ is not obvious a priori— $C$ is typically chosen according to a global criterion to ensure good performance at a specific machine learning task, without consideration for its effect on the spectrum or leverage scores of $A_{ij}^{(\sigma,\nu,C)}$.

To illustrate the diverse range of properties exhibited by these four classes of data sets, consider Table 5. Several observations are particularly relevant to our discussion below.

Both the Linear Kernels and the Dense RBF Kernels are much denser and are much more well-approximated by moderately to very low-rank matrices. In addition, both the Linear Kernels and the Dense RBF Kernels have statistical leverage scores that are much more uniform—there are several ways to illustrate this, none of them perfect, and here, we illustrate this by considering the $k^{th}$ largest leverage score, scaled by the factor $n/k$ (if $\mathbf{A}$ were exactly rank $k$, this would be the coherence of $\mathbf{A}$). For the Linear Kernels and the Dense RBF Kernels, this quantity is typically one to two orders of magnitude smaller than for the Laplacian Kernels.

For the Dense RBF Kernels, we consider two values of the $\sigma$ parameter, again chosen (somewhat) arbitrarily. For both AbaloneD and WineD, we see that decreasing $\sigma$ from $1$ to $0.15$, i.e., letting data points “see” fewer nearby points, has two important effects: first, it results in matrices that are much less well-approximated by low-rank matrices; and second, it results in matrices that have much more heterogeneous leverage scores. For example, for AbaloneD, the fraction of the Frobenius norm that is captured decreases from $97.8$ to $42.1$ and the scaled $k^{th}$ largest leverage score increases from $2.44$ to $18.11$.

For the Sparse RBF Kernels, there are a range of sparsities, ranging from above the sparsity of the sparsest Linear Kernel, but all are denser than the Laplacian Kernels. Changing the $\sigma$ parameter has the same effect (although it is even more pronounced) for Sparse RBF Kernels as it has for Dense RBF Kernels. In addition, “sparsifying” a Dense RBF Kernel also has the effect of making the matrix less well approximated by a low-rank matrix and of making the leverage scores more nonuniform. For example, for AbaloneD with $\sigma=1$ (respectively, $\sigma=0.15$), the fraction of the Frobenius norm that is captured decreases from $97.8$ (respectively, $42.1$) to $90.6$ (respectively, $15.4$), and the scaled $k^{th}$ largest leverage score increases from $2.44$ (respectively, $18.11$) to $3.57$ (respectively, $48.4$).

As we see below, when we consider the RBF Kernels as the width parameter and sparsity are varied, we observe a range of intermediate cases between the extremes of the (“nice”) Linear Kernels and the (very “non-nice”) Laplacian Kernels.

2. SPSD Sketching Algorithms

The sketching matrix $\mathbf{S}$ may be selected in a variety of ways. We will provide empirical results for two sampling-based SPSD sketches and two projection-based SPSD sketches. In the former case, the sketching matrix $\mathbf{S}$ contains exactly one nonzero in each column, corresponding to a single sample from the columns of $\mathbf{A}.$ In the latter case, $\mathbf{S}$ is dense, and mixes the columns of $\mathbf{A}$ before sampling from the resulting matrix.

In the case of leverage-based sampling, $\mathbf{S}$ has a more complicated distribution. Recall that the leverage scores relative to the best rank-$k$ approximation to $\mathbf{A}$ are the squared Euclidean norms of the rows of the $n\times k$ matrix $\mathbf{U}_{1}:$

In the figures, we refer to sketches constructed by selecting columns uniformly at random with the label ‘unif’, leverage score-based sketches with ‘lev’, Gaussian sketches with ‘gaussian’, and Fourier sketches with ‘srft’.

3. Effects of Data Analysis Preprocessing Decisions

inline]Talwalkar: You may want to include k-means as a competitor in Section 3.3? It lacks a solid theoretical grounding, but in our work , k-means was clearly the best method.

Before proceeding with our main empirical results, we pause to describe the effects of various machine learning and data analysis “design decisions” on the behavior of SPSD sketching algorithms in general as well as on the behavior of the statistical leverage scores in particular. We should emphasize that, for “worst case” matrices, very little can be said in this regard. Thus, these observations are based on our experiences with a diverse range of data sets, including those from Section 3.1. While not completely general, these observations are likely to hold in modified form for many other realistic data, and they can potentially be useful as heuristic guides to practice. For example, if preprocessing does not significantly change the leverage score distribution, then one could compute the leverage scores on the raw data and use these to sample columns from the processed data or to certify that the data have low coherence. Likewise, the behavior of the leverage scores as the rank parameter $k$ is varied or as the $\sigma$ scale parameter of RBF kernels varies is of interest, as it is expensive to compute the leverage scores anew for each value of $k$ or $\sigma$ as part of a cross-validation computation. inline]Incorporate Ilse’s point that the leverage scores could potentially be cheaply updated— it’s a reasonable, but not clear how best to do so, so leave as research direction

Another preprocessing decision has to do with the choice of rank $k$ with which to describe the data. This is typically determined according to an exogeneously-specified “model selection” criterion that does not explicitly take into account the spectrum or leverage score structure of the input matrix. It enters our discussion since we consider sampling columns with probabilities proportional to their statistical leverage scores relative to a rank-$k$ space, and thus the leverage scores depend on $k$. In our experience, increasing $k$ tends to uniformize or homogeneize the leverage scores, often gradually, but sometimes quite substantially. (We should note, however, that there are exceptions to this, where one observes very strong localization on low-order eigenvectors of data matrices .)

Yet another preprocessing decision has to do with the choice of the $\sigma$ scale parameter in Gaussian RBFK matrices. As with the rank parameter, the scale parameter $\sigma$ in practice is determined according to an exogeneously-specified model selection criterion that does not explicitly take into account the spectrum or leverage score structure of the input matrix. In our experience, as $\sigma$ increases, the leverage scores become more and more uniform; and they become more heterogeneous as $\sigma$ decreases. Informally, as a data point “sees” more data points, any outlying effect is mitigated. Varying $\sigma$ also has an effect on the spectrum. As a general rule, letting $\sigma\rightarrow 0$ tends to make the spectrum of $\mathbf{A}^{\sigma}$ flatter, i.e., decay more slowly, and letting $\sigma\rightarrow\infty$ makes $\mathbf{A}^{\sigma}$ lower-rank. Recall that the diagonal entries of $\mathbf{A}^{\sigma}$ are identically one, and as $\sigma\rightarrow\infty,$ $\mathbf{A}^{\sigma}$ tends to the matrix of all ones. That is, increasing $\sigma$ corresponds to considering all the observations $\mathbf{x}_{i}$ as being equally dissimilar, so all columns are equally noninformative. On the other hand, as $\sigma\rightarrow 0,$ $\mathbf{A}^{\sigma}$ approaches the identity, and very dissimilar observations (in the sense that $\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|_{2}$ is large) are penalized more heavily than similar observations, and thus there is some nonuniformity in the columns of $\mathbf{A}^{\sigma}.$ In some cases, we observed that, as the scale $\sigma$ decreases, the leverage scores stabilize, identifying the same columns as being important or influential over a range of scales.

4. Reconstruction Accuracy of Sampling and Projection Algorithms

Here, we describe the performances of the SPSD sketches described in Section 3.2—column sampling uniformly at random without replacement, column sampling according to the nonuniform leverage score probabilities, and sampling using Gaussian and SRFT mixtures of the columns—in terms of reconstruction accuracy for the data sets described in Section 3.1. We describe general observations we have made about each class of matrices in turn, and then we summarize our observations. We consider only the use of exact leverage scores here, and we postpone until Section 3.5 a discussion of running time issues and similar reconstruction results when approximate leverage scores are used for the importance sampling distribution. In each case, we present results for both the “non-rank-restricted” case as well as the “rank-restricted” case. Recall that by non-rank-restricted, we mean that the error

is plotted; while by rank-restricted, we mean that the error

is plotted (viz., the matrix $\mathbf{W}$ in Eqn. (7) has been replaced with the low-rank approximation $\mathbf{W}_{k}$). Note that previous work has shown that relative-error guarantees can be obtained, e.g., with CUR matrix decompositions, not only when one projects onto the span of judiciously-chosen columns, analogously to Eqn. (7) and as our worst-case guarantees in this paper are formulated, but also when one restricts the rank of the low-rank approximation to be no greater than $k$ by projecting onto the best rank-$k$ approximation to the original matrix . We evaluate the “rank-restricted” case of the form of Eqn. (8), that depends on projecting onto the best rank-$k$ approximation of the subsample (and not the original matrix) since it is more algorithmically tractable; but we note that similar but “smoother” results (e.g., the error is much more monotonic as a function of the number of samples, when compared with the “rank-restricted” results we present below) are obtained empirically with this more expensive rank-restriction procedure. The data points plotted in each figure of this section represent the average errors observed over 30 trials.

Finally, we note that previous work has shown that the statistical leverage scores reflect an important nonuniformity structure in the columns of general data matrices ; that randomly sampling columns according to this distribution results in lower worst-case error (for problems such as least-squares approximation and low-rank approximation of general matrices) than sampling columns uniformly at random ; and that leverage scores have proven useful in a wide range of practical applications . In spite of this, ours is the first work to implement and evaluate leverage score sampling for low-rank approximation of SPSD matrices.

These and subsequent figures contain a lot of information, some of which is peculiar to the given data sets and some of which is more general. In light of subsequent discussion, several observations are worth making about the results presented in these two figures.

All of the SPSD sketches provide quite accurate approximations—relative to the best possible approximation factor for that norm, and relative to bounds provided by existing theory, as reviewed in Section 2.4—even with only $k$ column samples (or in the case of the Gaussian and SRFT mixtures, with only $k$ linear combinations of vectors). Upon examination, this is partly due to the extreme sparsity and extremely slow spectral decay of these data sets which means, as shown in Table 4, that only a small fraction of the (spectral or Frobenius or trace) mass is captured by the optimal rank $20$ or $60$ approximation. Thus, although an SPSD sketch constructed from $20$ or $60$ vectors also only captures a small portion of the mass of the matrix, the relative error is small, since the scale of the residual error is large.

The scale of the Y axes is different between different figures and subfigures. This is to highlight properties within a given plot, but it can hide several things. In particular, note that the scale for the spectral norm is generally larger than for the Frobenius norm, which is generally larger than for the trace norm, consistent with the size of those norms; and that the scale is larger for higher-rank approximations, e.g. compare GR $k=20$ with GR $k=60$, also consistent with the larger amount of mass captured by higher-rank approximations.

All in all, there seems to be quite complicated behavior for low-rank sketches for these Laplacian data sets. Several of these observations can also be made for subsequent figures; but in some other cases the (very sparse and not very low rank) structural properties of the data are primarily responsible.

4.2. Linear Kernels

Figure 3 shows the reconstruction error results for sampling and projection methods applied to several Linear Kernels. The data sets (Dexter, Protein, SNPs, and Gisette) are all quite low-rank and have fairly uniform leverage scores. Several observations are worth making about the results presented in this figure.

The scale of the Y axes is much larger than for the Laplacian data sets, mostly since the matrices are much more well-approximated by low-rank matrices, although the scale decreases as one goes from spectral to Frobenius to trace reconstruction error, as before.

These linear kernels (and also to some extent the dense RBF kernels below that have larger $\sigma$ parameter) are examples of relatively “nice” machine learning data sets that are similar to matrices where uniform sampling has been shown to perform well previously ; and for these matrices our empirical results agree with these prior works.

4.3. Dense and Sparse RBF Kernels

Figure 4 and Figure 5 present the reconstruction error results for sampling and projection methods applied to several dense RBF and sparse RBF kernels. Several observations are worth making about the results presented in these figures.

Recall from Table 5 that for smaller values of $\sigma$ and for sparser kernels, the SPSD matrices are less well-approximated by low-rank matrices, and they have more heterogeneous leverage scores. Thus, they are more similar to the Laplacian data than the Linear Kernel data; and this suggests (as we have observed) that leverage score sampling should perform relatively better than uniform column sampling and projection-based schemes when in these two cases. In particular, nowhere do we see that leverage score sampling performs much worse than other methods, as we saw with the rank-restricted Linear Kernel results.

4.4. Summary of Comparison of Sampling and Projection Algorithms

Before proceeding, there are several summary observations that we can make about sampling versus projection methods for the data sets we have considered.

Linear Kernels and to a lesser extent Dense RBF Kernels with larger $\sigma$ parameter have relatively low-rank and relatively uniform leverage scores, and in these cases uniform sampling does quite well. These data sets correspond most closely with those that have been studied previously in the machine learning literature, and for these data sets our results are in agreement with that prior work.

Sparsifying RBF Kernels and/or choosing a smaller $\sigma$ parameter tends to make these kernels less well-approximated by low-rank matrices and to have more heterogeneous leverage scores. In general, these two properties need not be directly related—the spectrum is a property of eigenvalues, while the leverage scores are determined by the eigenvectors—but for the data we examined they are related, in that matrices with more slowly decaying spectra also often have more heterogeneous leverage scores.

For Dense RBF Kernels with smaller $\sigma$ and Sparse RBF Kernels, leverage score sampling tends to do much better than other methods. Interestingly, the Sparse RBF Kernels have many properties of very sparse Laplacian Kernels corresponding to relatively-unstructured informatics graphs, an observation which should be of interest for researchers who construct sparse graphs from data using, e.g., “locally linear” methods, to try to reconstruct hypothesized low-dimensional manifolds.

In general, all of the sampling and projection methods we considered perform much better on the SPSD matrices we considered than previous worst-case bounds (e.g., ) would suggest. (That is, even the worst results correspond to single-digit approximation factors in relative scale.) This observation is intriguing, because the motivation of leverage score sampling (and, recall, that in this context random projections should be viewed as performing uniform random sampling in a randomly-rotated basis where the leverage scores have been approximately uniformized ) is very much tied to the Frobenius norm, and so there is no a priori reason to expect its good performance to extend to the spectral or trace norms. Motivated by this, we revisit the question of proving improved worst-case theoretical bounds in Section 4.

Before describing these improved theoretical results, however, we address in Section 3.5 running time questions. After all, a naïve implementation of sampling with exact leverage scores is slower than other methods (and much slower than uniform sampling). As shown below, by using the recently-developed approximation algorithm of , not only does this approximation algorithm run in time comparable with random projections (for certain parameter settings), but it leads to approximations that soften the strong bias that the exact leverage scores provide toward the best rank-$k$ approximation to the matrix, thereby leading to improved reconstruction results in many cases.

5. Reconstruction Accuracy of Leverage Score Approximation Algorithms

Algorithm 1 (which originally appeared as Algorithm 1 in ) takes as input an arbitrary $n\times d$ matrix $\mathbf{A}$, where $n\gg d$, and it returns as output a $1\pm\epsilon$ approximation to all of the statistical leverage scores of the input matrix. The original algorithm of uses a subsampled Hadamard transform and requires $r_{1}$ to be somewhat larger than what we state in Algorithm 1. That an SRFT with a smaller value of $r_{1}$ can be used instead is a consequence of the fact that Lemma 3 in is also satisfied by an SRFT matrix with the given $r_{1};$ this is established in .

Consider, next, Algorithm 2 (which originally appeared as Algorithm 4 in ), which takes as input an arbitrary $n\times d$ matrix $\mathbf{A}$ and a rank parameter $k$, and returns as output a $1\pm\epsilon$ approximation to all of the statistical leverage scores (relative to the best rank-$k$ approximation) of the input. An important technical point is that the problem of computing the leverage scores of a matrix relative to a low-dimensional space is ill-posed, essentially because the spectral gap between the $k^{th}$ and the $(k+1)^{st}$ eigenvalues can be small, and thus Algorithm 2 actually computes approximations to the leverage scores of a matrix that is near to $\mathbf{A}$ in the spectral norm (or the Frobenius norm if $q=0$). See for details. Basically, this algorithm uses Gaussian sampling to find a matrix close to $\mathbf{A}$ in the Frobenius norm or spectral norm, and then it approximates the leverage scores of this matrix by using Algorithm 1 on the smaller, very rectangular matrix $\mathbf{B}$. When $\mathbf{A}$ is square, as in our applications, Algorithm 2 is typically more costly than direct computation of the leverage scores, at least for dense matrices (but it does have the advantage that the number of iterations is bounded, independent of properties of the matrix, which is not true for typical iterative methods to compute low-rank approximations).

Finally, note that although choosing the number of iterations $q$ as we did in Algorithm 2 is convenient for worst-case analysis, as a practical implementational matter it is easier either to choose $q$ based on spectral gap information revealed during the running of the algorithm or to prespecify $q$ to be a small integer, e.g., $2$ or $3$, before the algorithm runs. Both of these have an interpretation of accelerating the rate of decay of the spectrum with a power iteration, but they behave somewhat differently due to the different stopping conditions. Below, we consider both variants.

5.2. Running Time Comparisons

Uniform sampling is always less expensive and typically much less expensive than the other methods, while (with one minor exception) sampling according to the exact leverage scores is always the most expensive method.

The “fast Fourier” methods underlying the SRFT can take advantage of the structure of the Linear Kernels to yield algorithms that are similar to Gaussian projections and much better than exact leverage score computation. Note that the reason that SRFT is worse than Gaussians here is that the matrices we are considering are not extremely large, and we are not considering very large values of the rank parameter. Extending in both those directions leads to Gaussian projections being slower than SRFT, as the trends in the figures clearly indicate.

Gaussian projections are not too much slower than uniform sampling for the extremely sparse Laplacian Kernels—this is due to the sparsity of the Laplacian Kernels, since Gaussian projections can take advantage of the fast matrix-vector multiply, while the SRFT-based scheme cannot—but this advantage is lost for the (denser) Sparse RBF Kernels, to the extent that there is little running time improvement relative to the Dense RBF Kernels. In addition, Gaussian projections are relatively slower, when compared to the SRFT and uniform sampling, for the Dense RBF Kernels than for the Linear Kernels, although both of those data sets are maximally dense.

These approximate leverage score-based algorithms can be orders of magnitude faster than exact leverage score computation; but, especially for “spec levscore” when $q$ is not prespecified to be $2$ or $3$, they can even be somewhat slower. Exactly which is the case depends upon the properties of the matrix and the parameters used in the approximation algorithm, including especially the number of power iterations.

The “spec levscore” and “power” approximations with $q>0$ are more expensive than the $q=0$ “frob lev” approximation, which is a result of the relatively-expensive matrix-matrix multiplication. For the Linear Kernels, both are much better than the exact leverage score computation, and for most other data at least “power” is somewhat less expensive than the exact leverage score computation. For example, this is particularly true for the Laplacian Kernels.

Recall that the cost associated with these SPSD sketches is two-fold: first, the cost to construct the sample—by sampling columns uniformly at random, by computing a nonuniform importance sampling distribution, or by performing a random projection to uniformize the leverage scores; and second, the cost to construct the low-rank approximation from the sample. For uniform sampling, the latter step dominates the cost, while for more sophisticated methods the former step typically dominates the cost. The approximate leverage score sampling methods are still sufficiently expensive that the cost of computing the sampling probabilities still dominates the cost to construct the low-rank approximation.

Most importantly, the running time of Algorithm 1 on these rectangular matrices is faster than performing a QR decomposition on $\mathbf{A}$ and is comparable to applying a SRFT to $\mathbf{A}$. This is expected, since the running time bottleneck for Algorithm 1 is the application of the SRFT.

In addition, the running time of Algorithm 1 is significantly faster than the other approximate leverage score algorithms. This too is expected, since these other algorithms are applied to $\mathbf{A}$ and ignore the rectangular structure of $\mathbf{X}$.

Figure 9 shows that these improved running time gains for Algorithm 1 can come at the cost of a slight loss in the reconstruction accuracy (relative to the exact computation of the leverage scores) of the low-rank approximations; the accuracy of the other approximate leverage score algorithms is discussed in the following subsection.

5.3. Reconstruction Accuracy Results

Here, we describe the performances of the various low-rank approximations that use approximate leverage scores in terms of reconstruction accuracy for the data sets described in Section 3.1. The results are presented in Figure 10 through Figure 14. The setup for these results parallels that for the low-rank approximation results described in Section 3.4, and these figures parallel Figure 1 through Figure 5. To provide a baseline for the comparison, we also plot the previous reconstruction errors for sampling with the exact leverage scores as well as the uniform column sampling sketch. Several observations are worth making about the results presented in these figures.

For both the dense and the sparse RBF data sets, for the non-rank-restricted case, the approximate leverage score algorithms tend to parallel the exact leverage score algorithm, and they are not substantially better. In particular, both “power” and “spec levscore” tend to saturate when the exact method saturates, but in those cases “frob levscore” tends not to saturate.

inline]This paragraph was not clear to Ilse; clarify Note that the difference between different approximate leverage score algorithms often corresponds to a difference in the spectral gaps of the corresponding matrices. From Table 5, if we fix $k$ and use the approximate leverage scores filtered through rank $k$ to form a Nyström approximation to $\mathbf{A}$, the accuracy of that approximation has a strong dependence on the spectral gap of $\mathbf{A}$ at rank $k,$ as measured by $\frac{\lambda_{k}}{\lambda_{k+1}}.$ In general, the larger the spectral gap, the more accurate the approximation. This phenomena can also be understood in terms of the convergence of the approximate leverage scores: the approximation algorithms (in Algorithm 2 and Algorithm 3) are essentially truncated versions of the subspace iteration method for computing the top $k$ eigenvectors of $\mathbf{A}.$ It is a classical result that the spectral gap determines the rate of convergence of the subspace iteration process to the desired eigenvectors: the larger it is, the fewer iterations of the process are required to get accurate approximations of the top eigenvectors. It follows immediately that the larger the spectral gap, the more accurate the approximate leverage scores generated by these approximation algorithms are. Our empirical results illustrate the complexities and subtle consequences of these properties in realistic machine learning applications of even modestly-large size.

5.4. Summary of Leverage Score Approximation Algorithms

Before proceeding, there are several summary observations that we can make about the running time and reconstruction quality of approximate leverage score sampling algorithms for the data sets we have considered.

The running time of computing the exact leverage scores is generally much worse than that of uniform sampling and both SRFT-based and Gaussian-based random projection methods.

The running time of computing approximations to the leverage scores can, with appropriate choice of parameters, be much faster than the exact computation of the leverage scores; and, especially for “frob levscore,” can be comparable to the running time of the random projection (SRFT or Gaussian) used in the leverage score approximation algorithm. For the methods that involve $q>0$ iterations to compute stronger approximations to the leverage scores, the running time can vary considerably depending on details of the stopping condition.

The approximate leverage scores computed from “power” and “spec levscore” approach those of the exact leverage scores, as $q$ is increased; and they obtain reconstruction accuracy that is no worse, and in many cases is better, than that obtained by the exact leverage scores. This suggests that, by not fitting exactly to the empirical statistical leverage scores, we are observing a form of implicit regularization.

The running time of Algorithm 1, when applied to “tall” matrices for which $n\gg d$, is faster than the running time of performing a QR decomposition of the matrix $\mathbf{A}$; and it is comparable to the running time of applying a random projection to $\mathbf{A}$ (which is the computational bottleneck of applying Algorithm 1). Thus, in particular, one could use this algorithm to compute approximations to the leverage scores to obtain a sketch that provides a relative-error approximation to a least-squares problem involving $\mathbf{A}$ ; or one could use the sketch thereby obtained as a preconditioner to an iterative method to solve the least-squares problem, in a manner analogous to how Blendenpik or LSRN do so with a random projection .

Previous work has showed that one can implement random projection algorithms to provide low-rank approximations with error comparable to that of the SVD in less time than state-of-the art Krylov solvers and other “exact” numerical methods . Our empirical results show that these random projection algorithms can be used in two complementary ways to approximate SPSD matrices of interest in machine learning: first, they can be used directly to compute a projection-based low-rank approximation; and second, they can be used to compute approximations to the leverage scores, which can be used to compute a sampling-based low-rank approximation. With the right choice of parameters, the two complementary approaches have roughly comparable running times, and neither one dominates the other in terms of reconstruction accuracy.

6. Projection-based Sketches

Finally, for completeness, we consider the performance of the two projection-based SPSD sketches proposed in , and we show how they perform when compared with the sketches we have considered. Recall that the idea of these sketches is to construct low-rank approximations by forming an approximate basis $\mathbf{Q}$ for the top eigenspace of $\mathbf{A}$ and then restricting $\mathbf{A}$ to that eigenspace. In more detail, given a sketching matrix $\mathbf{S},$ form the matrix $\mathbf{Y}=\mathbf{A}\mathbf{S}$ and take the QR decomposition of $\mathbf{Y}$ to obtain $\mathbf{Q},$ a matrix with orthonormal columns. The first sketch, which we eponymously refer to as the pinched sketch, is simply $\mathbf{A}$ pinched to the space spanned by $\mathbf{Q}:$

The second sketch, which we refer to as the prolonged sketch, is

It is clear that the prolonged sketch can be constructed using our SPSD Sketching Model by taking $\mathbf{Q}$ as the sketching matrix. In fact, a stronger statement can be made. As shown in , and as stated in Lemma 1 below, it is the case, for any sketching matrix $\mathbf{X},$ that when $\mathbf{C}=\mathbf{A}\mathbf{X}$ and $\mathbf{W}=\mathbf{X}^{T}\mathbf{A}\mathbf{X},$

By considering the two choices $\mathbf{X}=\mathbf{A}\mathbf{S}$ and $\mathbf{X}=\mathbf{Q}$, we see that in fact the prolonged sketch is exactly the sketch obtained by applying the power method with $q=2:$

It follows that the bounds we provide in Section 4 on the performance of sketches obtained using the power method pertain also to prolonged sketches.

In Figure 15, we compare the empirical performances of several of the SPSD sketches considered earlier with their pinched and prolonged variants. Specifically, we plot the errors of pinched and prolonged sketches for several choices of sketching matrices—corresponding to uniform column sampling, gaussian column mixtures, and SRFT-based column mixtures—along with the errors of non-pinched, non-prolonged sketches constructing using the same choices of $\mathbf{S}.$ In the interest of brevity, we provide results only for several of the datasets listed in Table 4, and we consider only the nonfixed-rank variants of the sketches.

In the spectral norm, the prolonged sketches are considerably more accurate than the pinched and standard sketches for all the datasets considered. Without exception, the prolonged Gaussian and SRFT column-mixture sketches are the most accurate in the spectral norm, of all the sketches considered. Only in the case of the Dexter Linear Kernel is the prolonged uniformly column-sampled sketch nearly as accurate in the spectral norm as the prolonged Gaussian and SRFT sketches. To a lesser extent, the prolonged sketches are also more accurate in the Frobenius and trace norms than the other sketches considered. The increased Frobenius and trace norm accuracy is particularly notable for the two RBF Kernel datasets; again, the prolonged Gaussian and SRFT sketches are considerably more accurate than the prolonged uniformly column-sampled sketches.

After the prolonged sketches, the pinched Gaussian and SRFT column-mixture sketches exhibit the least spectral, Frobenius, and trace norm errors. Again, however, we see that the pinched uniformly column-sampled sketches are considerably less accurate than the pinched Gaussian and SRFT column-mixture sketches. Particularly in the spectral and Frobenius norms, the pinched uniformly column-sampled sketches are not any more accurate than the basic uniformly column-sampled sketches.

From these considerations, it seems evident that the benefits of pinched and prolonged sketches are most prominent when the spectral norm is the error metric, or when the dataset is an RBF Kernel. In particular, pinched and prolonged sketches are not significantly more accurate (than the sketches considered in the previous subsections) in the Frobenius and trace norms for any of the datasets considered.

It is also evident from Figure 15 that the pinched sketches often have a much slighter increase in accuracy over the basic sketches than do the prolonged sketches. To understand why the pinched sketches are less accurate than the prolonged sketches, observe that the pinched sketches satisfy

while, as noted above, the prolonged sketches can be written in the form

Thus, pinched and prolonged sketches approximate the square root of $\mathbf{A}$ by projecting, respectively, onto the ranges of $\mathbf{A}\mathbf{S}$ and $\mathbf{A}^{3/2}\mathbf{S}.$ The spectral decay present in $\mathbf{A}$ is increased when $\mathbf{A}$ is raised to a power larger than one; consequently, the range of $\mathbf{A}^{3/2}\mathbf{S}$ is more biased towards the top $k$-dimensional invariant subspace of $\mathbf{A}$ than is the range of $\mathbf{A}\mathbf{S}.$ It follows that the approximate square root used to construct the prolonged sketches more accurately captures the top $k$-dimensional subspace of $\mathbf{A}$ than does that used to construct the pinched sketches.

Theoretical Aspects of SPSD Low-rank Approximation

In this section, we present our main theoretical results, which consist of a suite of bounds on the quality of low-rank approximation under several different sketching methods. As mentioned above, these were motivated by our empirical observation that all of the sampling and projection methods we considered perform much better on the SPSD matrices we considered than previous worst-case bounds (e.g., ) would suggest. We start in Section 4.1 with deterministic structural conditions for the spectral, Frobenius, and trace norms; and then in Section 4.2 we use these results to provide our bounds for several random sampling and random projection procedures.

In this section, we present three theorems that provide error bounds for the spectral, Frobenius, and trace norm approximation errors under the SPSD Sketching Model of Section 2.2. These are provided in Sections 4.1.1, 4.1.3, and 4.1.5, respectively, and they are followed by several more general remarks in Section 4.1.6. Note that these bounds hold for any, e.g., deterministic or randomized, sketching matrix $\mathbf{S}$. Thus, e.g., one could use them to check, in an a posteriori manner, the quality of a sketching method for which one cannot establish an a priori bound. Rather than doing this, we use these results (in Section 4.2 below) to derive a priori bounds for when the sketching operation consists of common random sampling and random projection algorithms.

Our results are based on the fact, established in , that approximations which satisfy our SPSD Sketching Model can be written in terms of a projection onto a subspace of the range of the square root of the matrix being approximated. The following fact appears in the proof of Proposition 1 in .

inline]Update this reference to unreleased version in anticipation of changing on arxiv

Let $\mathbf{A}$ be an SPSD matrix and $\mathbf{S}$ be a conformal sketching matrix. Then when $\mathbf{C}=\mathbf{A}\mathbf{S}$ and $\mathbf{W}=\mathbf{S}^{T}\mathbf{A}\mathbf{S},$ the corresponding low-rank SPSD approximation satisfies

We start with a bound on the spectral norm of the residual error. Although this result is trivial to prove, given prior work, it highlights several properties that we use in the analysis of our subsequent results.

assuming $\boldsymbol{\Omega}_{1}$ has full row rank.

Apply Lemma 1 with the sampling matrix $\mathbf{S}^{\prime}=\mathbf{A}^{q-1}\mathbf{S}$ (where, recall, $q\geq 1$) to see that

Next, recall that $\boldsymbol{\Omega}_{i}=\mathbf{U}_{i}^{T}\mathbf{S}$ and that $\mathbf{A}^{1/2}$ has eigenvalue decomposition $\mathbf{A}=\mathbf{U}\boldsymbol{\Sigma}^{1/2}\mathbf{U}^{T},$ where

It can be shown ([32, Theorems 9.1 and 9.2]) that, because $\boldsymbol{\Omega}_{1}$ has full row rank,

The latter inequality follows from the fact that the $2q-1$ radical function is subadditive when $q\geq 1$ and the identity $\left\|\boldsymbol{\Sigma}_{2}^{q-1/2}\right\|_{2}^{2}=\left\|\boldsymbol{\Sigma}_{2}\right\|_{2}^{2q-1}.$ This establishes the stated bound. ∎

Remark. The assumption that $\boldsymbol{\Omega}_{1}$ has full row rank is very non-trivial. It is, however, satisfied by our algorithms below. See Section 4.1.6 for more details on this point.

Remark. The proof of Theorem 1 proceeds in two steps. The first step relates low-rank approximation of an SPSD matrix $\mathbf{A}$ under the SPSD Sketching Model of Section 2.2 to column sketching (e.g., sampling or projecting) from the square-root of $\mathbf{A}$. A weaker relation of this type was used in , but the stronger form that we use here in Equation (11) was first proved in . The second step is to use a deterministic structural result that holds for sampling/projecting from an arbitrary matrix. The structural bound of the form of Equation (12) was originally proven for $q=1$ in , where it was applied to the Column Subset Selection Problem. The bound was subsequently improved in , where it was applied to a random projection algorithm and extended to apply when $q>1$. Although the analyses of our next two results are more complicated, they follow the same high-level two-step approach.

Before proceeding with the analogous Frobenius and trace norm bounds, we pause to describe a geometric interpretation of this result.

inline]Explain interpretation of $\boldsymbol{\Omega}_{2}\boldsymbol{\Omega}_{1}^{\dagger}$, how it shows up in convergence of power method, and how it being small implies $\mathbf{S}$ is an effective sketching matrix.

1.2. A geometric interpretation of the sketching interaction matrix

It is evident from the bound in Theorem 1 that the smaller the spectral norm of the sketching interaction matrix $\boldsymbol{\Omega}_{2}\boldsymbol{\Omega}_{1}^{\dagger},$ the more effective $\mathbf{S}$ is as a sketching matrix. If, additionally, the columns of $\mathbf{S}$ are orthonormal, we can give this norm a natural geometric interpretation as the tangent of the largest angle between the spaces spanned by $\mathbf{S}$ and $\mathbf{U}_{1}.$

To verify our claim, we first recall the definition of the sine between the range spaces of two matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}:$

Note that this quantity is not symmetric: it measures how well the range of $\mathbf{M}_{1}$ captures that of $\mathbf{M}_{2}$ [29, Chapter 12]. Since $\mathbf{U}_{1}$ and $\mathbf{S}$ (by assumption here) have orthonormal columns, we see that

The second to last equality holds because $\mathbf{U}_{1}^{T}\mathbf{S}$ has $k$ rows and we assumed it has full row rank. Accordingly, it follows that

The second to last equality holds because of the fact that, for any matrix $\mathbf{M},$

this identity can be established with a routine SVD argument.

Thus, when $\mathbf{S}$ has orthonormal columns and $\mathbf{U}_{1}^{T}\mathbf{S}$ has full row-rank, $\left\|\boldsymbol{\Omega}_{2}\boldsymbol{\Omega}_{1}\right\|_{2}$ is the tangent of the largest angle between the range of $\mathbf{S}$ and the eigenspace spanned by $\mathbf{U}_{1}.$ If $\mathbf{U}_{1}^{T}\mathbf{S}$ does not have full row-rank, then our derivation above shows that $\sin^{2}(\mathbf{S},\mathbf{U}_{1})=1,$ meaning that there is a vector in the eigenspace spanned by $\mathbf{U}_{1}$ which has no component in the space spanned by the sketching matrix $\mathbf{S}.$

We note that $\tan(\mathbf{S},\mathbf{U}_{1})$ also arises in the classical bounds on the convergence of the orthogonal iteration algorithm for approximating invariant subspaces of a matrix (see, e.g. [29, Theorem 8.2.2]).

1.3. Frobenius Norm Bounds

Next, we state and prove the following bound on the Frobenius norm of the residual error. The proof parallels that for the spectral norm bound, in that we divide it into two analogous parts, but the analysis is somewhat more complex.

The multiplicative eigengap $\gamma=\lambda_{k+1}(\mathbf{A})/\lambda_{k}(\mathbf{A})$ that appears in the statement of this theorem predicts the effect of using the power method when constructing sketches. Specifically, the additional errors of sketches constructed using $\mathbf{C}=\mathbf{A}^{q}\mathbf{S}$ are at least a factor of $\gamma^{q-1}$ times smaller than those constructed using $\mathbf{C}=\mathbf{A}\mathbf{S}.$

Then when $\mathbf{C}=\mathbf{A}^{q}\mathbf{S}$ and $\mathbf{W}=\mathbf{S}^{T}\mathbf{A}^{2q-1}\mathbf{S}$, the corresponding low-rank SPSD approximation satisfies

assuming $\boldsymbol{\Omega}_{1}$ has full row rank.

Apply Lemma 1 with the sampling matrix $\mathbf{S}^{\prime}=\mathbf{A}^{q-1}\mathbf{S}$ to see that

To bound this quantity, we first use the unitary invariance of the Frobenius norm and the fact that

By construction, $\mathbf{Z}$ has full column rank, thus $\mathbf{Z}(\mathbf{Z}^{T}\mathbf{Z})^{-1/2}$ is an orthonormal basis for the span of $\mathbf{Z}$, and

Next, we provide bounds for $T_{1}$, $T_{2}$, and $T_{3}$. Using the fact that $\mathbf{0}\preceq\mathbf{I}-\mathbf{F}(\mathbf{I}+\mathbf{F}^{T}\mathbf{F})^{-1}\mathbf{F}^{T}\preceq\mathbf{I},$ we can bound $T_{3}$ with

Likewise, the fact that $\mathbf{I}-(\mathbf{I}+\mathbf{F}^{T}\mathbf{F})^{-1}\preceq\mathbf{F}^{T}\mathbf{F}$ (easily seen with an SVD) implies that we can bound $T_{1}$ as

We proceed to bound $T_{2}$ by using the estimate

To develop the term involving a spectral norm, observe that for any SPSD matrix $\mathbf{M}$ with eigenvalue decomposition $\mathbf{M}=\mathbf{V}\mathbf{D}\mathbf{V}^{T},$

Using this estimate in Equation (16), we conclude that

Combining our estimates for $T_{1},$ $T_{2},$ and $T_{3}$ with Equation (15) gives

The claimed bound follows by identifying $\gamma$ and applying the subadditivity of the square-root function:

As before, we pause to describe a geometric interpretation of this result.

1.4. Another geometric interpretation of the sketching interaction matrix

Just as the spectral norm of the spectral interaction matrix is the tangent of the largest angle between the range of the sketching matrix and the dominant $k$-dimensional eigenspace of $\mathbf{A}$, the Frobenius norm of the spectral interaction matrix has a geometric interpretation. To see this, recall that the principal angles between the ranges of the matrices $\mathbf{S}$ and $\mathbf{U}_{1}$ are defined recursively by

The proof of this claim hinges on the fact that, for any matrix $\mathbf{M},$

as can be readily verified with an SVD argument. From this observation, we have that

so the theorem guarantees that the additional error is on the order of $\sqrt{\left\|\boldsymbol{\Sigma}_{2}\right\|_{2}\left\|\boldsymbol{\Sigma}_{2}\right\|_{\star}}.$ This is an upper bound on the optimal Frobenius error:

We see, in particular, that if the residual spectrum is flat, i.e. $\lambda_{k+1}(\mathbf{A})=\cdots=\lambda_{n}(\mathbf{A}),$ then equality holds and the additional error is on the scale of the optimal error.

inline]It’s worth asking if this bound can be replaced with the optimal error itself.

1.5. Trace Norm Bounds

Finally, we state and prove the following bound on the trace norm of the residual error. The proof method is analogous to that for the spectral and Frobenius norm bounds.

As in the case of the Frobenius norm error, we see that the multiplicative eigengap $\gamma=\lambda_{k+1}(\mathbf{A})/\lambda_{k}(\mathbf{A})$ predicts the effect of using the power method when constructing sketches: the additional errors of sketches constructed using $\mathbf{C}=\mathbf{A}^{q}\mathbf{S}$ are a factor of $\gamma^{2(q-1)}$ times smaller than the additional errrors of those constructed using $\mathbf{C}=\mathbf{A}\mathbf{S}.$

Then when $\mathbf{C}=\mathbf{A}^{q}\mathbf{S}$ and $\mathbf{W}=\mathbf{S}^{T}\mathbf{A}^{2q-1}\mathbf{S}$, the corresponding low-rank SPSD approximation satisfies

assuming $\boldsymbol{\Omega}_{1}$ has full row rank.

Since $\mathbf{A}-\mathbf{C}\mathbf{W}^{\dagger}\mathbf{C}^{T}=\mathbf{A}^{1/2}(\mathbf{I}-\mathbf{P}_{\mathbf{A}^{q-1/2}\mathbf{S}})\mathbf{A}^{1/2}\succeq\mathbf{0},$ its trace norm simplifies to its trace. Thus

where $\mathbf{Z}=\begin{pmatrix}\mathbf{I}\\ \mathbf{F}\end{pmatrix}$ is defined in Equation (13). The expression for $\mathbf{I}-\mathbf{P}_{\mathbf{Z}}$ given in Equation (14) implies that

Recall the estimate $\mathbf{I}-(\mathbf{I}+\mathbf{F}^{T}\mathbf{F})^{-1}\preceq\mathbf{F}^{T}\mathbf{F}$ and the basic estimate $\mathbf{I}-\mathbf{F}(\mathbf{I}+\mathbf{F}^{T}\mathbf{F})^{-1}\mathbf{F}^{T}\preceq\mathbf{I}.$ Together these imply that

The first equality follows from substituting the definition of $\mathbf{F}$ and identifying the squared Frobenius norm. The last equality follows from identifying $\gamma.$ We have established the claimed bound.

1.6. Additional Remarks on Our Deterministic Structural Results.

Before applying these deterministic structural results in particular randomized algorithmic settings, we pause to make several additional remarks about these three theorems.

First, for some randomized sampling schemes, it may be difficult to obtain a sharp bound on \big{\|}\boldsymbol{\Omega}_{2}\boldsymbol{\Omega}_{1}^{\dagger}\big{\|}_{\xi} for $\xi=2,F$. In these situations, the bounds on the excess error supplied by Theorems 1, 2, and 3 may be quite pessimistic. On the other hand, since $\mathbf{A}-\mathbf{C}\mathbf{W}^{\dagger}\mathbf{C}^{T}=\mathbf{A}^{1/2}(\mathbf{I}-\mathbf{P}_{\left(\mathbf{A}^{1/2}\right)^{2q-1}\mathbf{S}})\mathbf{A}^{1/2}$, it follows that $\mathbf{0}\preceq\mathbf{A}-\mathbf{C}\mathbf{W}^{\dagger}\mathbf{C}^{T}\preceq\mathbf{A}$. This implies that the errors of any approximation generated used the SPSD Sketching Model, deterministic or randomized, satisfy at least the crude bound $\left\|\mathbf{A}-\mathbf{C}\mathbf{W}^{\dagger}\mathbf{C}^{T}\right\|_{\xi}\leq\left\|\mathbf{A}\right\|_{\xi}$.

Second, we emphasize that these theorems are deterministic structural results that bound the additional error (beyond that of the optimal rank-$k$ approximation) of low-rank approximations which follow our SPSD sketching model. That is, there is no randomness in their statement or analysis. In particular, these bounds hold for deterministic as well as randomized sketching matrices $\mathbf{S}$. In the latter case, the randomness enters only through $\mathbf{S}$, and one needs to show that the condition that $\boldsymbol{\Omega}_{1}$ has full row rank is satisfied with high probability; conditioned on this, the quality of the bound is determined by terms that depend on how the sketching matrix interacts with the subspace structure of the matrix $\mathbf{A}$.

In particular, we remind the reader that (although it is beyond the scope of this paper to explore this point in detail) these deterministic structural results could be used to check, in an a posteriori manner, the quality of a sketching method for which one cannot establish an a priori bound.

Third, we also emphasize that the assumption that $\boldsymbol{\Omega}_{1}$ has full row rank (equivalently, that $\tan(\mathbf{S},\mathbf{U}_{1})<\infty$) is very non-trivial; and that it is false, in worst-case at least and for non-trivial parameter values, for common sketching methods such as uniform sampling. To see that some version of leverage-based sampling is needed to ensure this condition, recall that $\mathbf{U}_{1}^{T}\mathbf{U}_{1}=\mathbf{I}$ and thus that $\boldsymbol{\Omega}_{1}\boldsymbol{\Omega}_{1}^{T}=\mathbf{U}_{1}^{T}\mathbf{S}\mathbf{S}^{T}\mathbf{U}_{1}$ can be viewed as approximating $\mathbf{I}$ with a small number of rank-$1$ components of $\mathbf{U}_{1}^{T}\mathbf{U}_{1}$. The condition that $\boldsymbol{\Omega}_{1}$ has full row rank is equivalent to $\left\|\mathbf{U}_{1}^{T}\mathbf{U}_{1}-\mathbf{U}_{1}^{T}\mathbf{S}\mathbf{S}^{T}\mathbf{U}_{1}\right\|_{2}<1$. Work on approximating the product of matrices by random sampling shows that to obtain non-trivial bounds one must sample with respect to the norm of the rank-$1$ components , which here (since we are approximating the product of two orthogonal matrices) equal the statistical leverage scores. From this perspective, random projections satisfy this condition since (informally) they rotate to a random basis where the leverage scores of the rotated matrix are approximately uniform and thus where uniform sampling is appropriate .

Finally, as observed recently , methods that use knowledge of a matrix square root $\Phi$ (i.e., a $\Phi$ such that $\mathbf{A}=\Phi\Phi^{T}$) typically lead to $\Omega(n^{2})$ complexity. An important feature of our approach is that we only use the matrix square root implicitly—that is, inside the analysis, and not in the statement of the algorithm—and thus we do not incur any such cost.

2. Stochastic Error Bounds for Low-rank SPSD Approximation

In this section, we apply the three theorems from Section 4.1 to bound the reconstruction errors for several random sampling and random projection methods that conform to our SPSD Sketching Model. In particular, we consider two variants of random sampling and two variants of random projections: sampling columns according to an importance sampling distribution that depends on the statistical leverage scores (in Section 4.2.1); randomly projecting by using subsampled randomized Fourier transformations (in Section 4.2.2); randomly projecting by uniformly sampling from Gaussian mixtures of the columns (in Section 4.2.3); and, finally, sampling columns uniformly at random (in Section 4.2.4).

The results are presented for the general case of SPSD sketches constructed using the power method, i.e., sketches constructed using $\mathbf{C}=\mathbf{A}^{q}\mathbf{S}$ for a positive integer $q>1.$ The additive errors of these sketches decrease proportionally to the number of iterations $q,$ where the constant of proportionality is given by the multiplicative eigengap $\gamma=\lambda_{k+1}(\mathbf{A})/\lambda_{k}(\mathbf{A}).$ Accordingly, the bounds involve the terms $\gamma^{q-1}$ and $\gamma^{2(q-1)}.$ The bounds simplify considerably when $q=1$ (i.e., when there are no additional iterations) or $\gamma=1$ (i.e., when there is no eigengap). In either of these cases, the terms $\gamma^{q-1}$ and $\gamma^{2(q-1)}$ all become the constant 1.

In particular, for random sampling algorithms that use a leverage-based importance sampling distribution, as we use in Section 4.2.1, it is often said that the running time is no faster than that of computing $\mathbf{U}_{1}$. (This $O(n^{2}k)$ running time claim is simply the running time of the naïve algorithm that computes $\mathbf{U}_{1}$ “exactly,” e.g., with a variant of the QR decomposition, and then reads off the Euclidean norms of the rows.) However, the randomized algorithm of that computes relative-error approximations to all of the statistical leverage in a time that is qualitatively faster—in worst-case theory and, by using existing high-quality randomized numerical code , in practice—gets around this bottleneck, as was shown in Section 3. The computational bottleneck for the algorithms of is that of applying a random projection, and thus the running time for leverage-based Nyström extension is that of applying a (“fast” Fourier-based or “slow” Gaussian-based, as appropriate) random projection to $\mathbf{A}$ . See Section 3 or for additional details.

Here, the columns of $\mathbf{A}$ are sampled with replacement according to a nonuniform probability distribution determined by the (exact or approximate) statistical leverage scores of $\mathbf{A}$ relative to the best rank-$k$ approximation to $\mathbf{A}$, which in turn depend on nonuniformity properties of the top $k$-dimensional eigenspace of $\mathbf{A}$. To add flexibility (e.g., in case the scores are computed only approximately with the fast algorithm of ), we formulate the following lemma in terms of any probability distribution that is $\beta$-close to the leverage score distribution. In particular, consider any probability distribution satisfying

for some $\beta\in(0,1].$ Fix a failure probability $\delta\in(0,1]$ and approximation factor $\epsilon\in(0,1],$ and let

simultaneously with probability at least $1-6\delta-0.6.$

with probability at least $1-2\delta-0.2.$ Thus, the estimates

each hold, individually, with probability at least $1-2\delta-0.2.$ In particular, taking $q=1$, we see that

These three estimates used in Theorems 1, 2, and 3 yield the bounds given in the statement of the theorem. ∎

inline]update this comment: is it still true? Remark. The additive scale factors for the spectral and Frobenius norm bounds are much improved relative to the prior results of . At root, this is since the leverage score importance sampling probabilities highlight structural properties of the data (e.g., how to satisfy the condition in Theorems 1, 2, and 3 that $\boldsymbol{\Omega}_{1}$ has full row rank) in a more refined way than the importance sampling probabilities of .

Remark. These improvements come at additional computational expense, but we remind the reader that leverage-based sampling probabilities of the form used by Lemma 2 can be computed faster than the time needed to compute the basis $\mathbf{U}_{1}$ . The computational bottleneck of the algorithm of is the time required to perform a random projection on the input matrix.

Remark. Not surprisingly, constant factors such as $3200$ (as well as other similarly large factors below) and a failure probability bounded away from zero are artifacts of the analysis; the empirical behavior of this sampling method is much better. This has been observed previously .

2.2. Random Projections with Subsampled Randomized Fourier Transforms

simultaneously with probability at least $1-2\delta.$

each hold, individually, with probability at least $1-\delta.$ These estimates used in Theorems 1, and 3 yield the stated bounds for the spectral and trace norm errors.

The Frobenius norm bound follows from the same estimates and a simplification of the bound stated in Theorem 2:

We note that a direct application of Theorem 2 gives a potentially tighter, but more unwieldy, bound. ∎

This should be compared to the guarantee established in Lemma 4 below for Gaussian-based SPSD sketches constructed using the same number of measurements:

2.3. Random Projections with i.i.d. Gaussian Random Matrices

simultaneously with probability at least $1-2k^{-1}-4k^{-k/\epsilon^{2}}.$

with probability at least $1-k^{-1}-2k^{-k/\epsilon^{2}}$ and

These estimates used in Theorems 1 and 3 yield the stated spectral and trace norm bounds. To obtain the corresponding Frobenius norm bound, define the quantities

The following estimates hold for the coefficients in this inequality:

The Frobenius norm bound follows from using these estimates in Equation (22) and grouping terms appropriately:

Remark. The way we have parameterized these bounds for Gaussian-based projections makes explicit the dependence on various parameters, but hides the structural simplicity of these bounds. In particular, note that the Frobenius norm bound is upper bounded by a term that depends on the Frobenius norm of the error and a term that depends on the trace norm of the error; and that, similarly, the trace norm bound is upper bounded by a multiplicative factor that can be set to $1+\epsilon$ with an appropriate choice of parameters.

2.4. Sampling Columns Uniformly at Random

Here, the columns of $\mathbf{A}$ are sampled uniformly at random (with or without replacement). Such uniformly-at-random column sampling only makes sense when the leverage scores of the top $k$-dimensional invariant subspace of the matrix are sufficiently uniform that no column is significantly more informative than the others. For this case, we can prove the following.

simultaneously with probability at least $1-3\delta.$

with probability at least $1-\delta.$ Also,

where the summands $\mathbf{x}_{i}$ are distributed uniformly at random over the columns of $\boldsymbol{\Sigma}_{2}^{1/2}\mathbf{U}_{2}^{T}.$ Regardless of whether $\mathbf{S}$ selects the columns with replacement or without replacement, the summands all have the same expectation:

Now applying Markov’s inequality to (23), we see that

with probability at least $1-2\delta.$ Thus, we also know that

also with probability at least $1-2\delta.$ These estimates used in Theorems 1 and 3 yield the stated spectral and trace norm bounds.

To obtain the Frobenius norm bound, observe that Theorem 2 implies

Remark. As with previous bounds for uniform sampling, e.g., , these results for uniform sampling are much weaker than our bounds from the previous subsections, since the sampling complexity depends on the coherence of the input matrix. When the matrix has small coherence, however, these bounds are similar to the bounds derived from the leverage-based sampling probabilities. Recall that, by the algorithm of , the coherence of an arbitrary input matrix can be computed in roughly the time it takes to perform a random projection on the input matrix.

Discussion and Conclusion

We have presented a unified approach to a large class of low-rank approximations of Laplacian and kernel matrices that arise in machine learning and data analysis applications; and in doing so we have provided qualitatively-improved worst-case theory and clarified the performance of these algorithms in practical settings. Our theoretical and empirical results suggest several obvious directions for future work.

In addition, we should note that, in situations where one is concerned with the quality of approximation of the actual eigenspaces, one desires both a small spectral norm error (because by the Davis–Kahan sin$\Theta$ theorem and similar perturbation results, this would imply that the range space of the sketch effectively captures the top $k$-dimensional eigenspace of $\mathbf{A}$) as well as to use as few samples as possible (because one prefers to approximate the top $k$-dimension eigenspace of $\mathbf{A}$ with as close to a $k$-dimensional subspace as possible). Our results suggest that the leverage score probabilities supply the best sampling scheme for balancing these two competing objectives.

More generally, although our empirical evaluation consists of random sampling and random projection algorithms, our theoretical analysis clearly decouples the randomness in the algorithm from the structural heterogenities in the Euclidean vector space that are responsible for the poor performance of uniform sampling algorithms. Thus, if those structural conditions can be satisfied with a deterministic algorithm, an iterative algorithm, or any other method, then one can certify (after running the algorithm) that good approximation guarantees hold for particular input matrices in less time than is required for general matrices. Moreover, this structural decomposition suggests greedy heuristics—e.g., greedily keep some number of columns according to approximate statistical leverage scores and “residualize.” In our experience, a procedure of this form often performs quite well in practice, although theoretical guarantees tend to be much weaker; and thus we expect that, when coupled with our results, such procedures will perform quite well in practice in many medium-scale and large-scale machine learning applications.

Acknowledgments. AG would like to acknowledge the support, under the auspice of Joel Tropp, of ONR awards N00014-08-1-0883 and N00014-11-1-0025, AFOSR award FA9550-09-1-0643, and a Sloan Fellowship; and MM would like to acknowledge a grant from the Defense Advanced Research Projects Agency.

References