A Bootstrap Method for Error Estimation in Randomized Matrix Multiplication

Introduction

The development of randomized numerical linear algebra (RNLA or RandNLA) has led to a variety of efficient methods for solving large-scale matrix problems, such as matrix multiplication, least-squares approximation, and low-rank matrix factorization, among others (Halko et al. 2011; Mahoney 2011; Woodruff 2014; Drineas and Mahoney 2016). A general feature of these methods is that they apply some form of randomized dimension reduction to an input matrix, which reduces the cost of subsequent computations. In exchange for the reduced cost, the randomization leads to some error in the resulting solution, and consequently, there is a tradeoff between cost and accuracy.

For many canonical matrix problems, the relationship between cost and accuracy has been the focus of a growing body of theoretical work, and the literature provides many performance guarantees for RNLA methods. In general, these guarantees offer a good qualitative description of how the accuracy depends on factors such as problem size, number of iterations, condition numbers, and so on. Yet, it is also the case that such guarantees tend to be overly pessimistic for any particular problem instance — often because the guarantees are formulated to hold in the worst case among a large class of possible inputs. Likewise, it is often impractical to use such guarantees to determine precisely how accurate a given solution is, or precisely how much computation is needed to achieve a desired level of accuracy.

In light of this situation, it is of interest to develop efficient methods for estimating the exact relationship between the cost and accuracy of RNLA methods on a problem-specific basis. Since the literature has been somewhat quiet on this general question, the aim of this paper is to analyze randomized matrix multiplication as a prototype setting, and propose an approach that may be pursued more broadly. (Extensions are discussed at the end of the paper in Section 6.)

2 Problem formulation

For example, the quantity $q_{0.99}(t)$ is the tightest upper bound on $\varepsilon_{t}$ that holds with probability at least 0.99. Hence, for any fixed $\alpha$, the function $q_{1-\alpha}(t)$ represents a precise tradeoff curve for relating cost and accuracy. Moreover, the function $q_{1-\alpha}(t)$ is specific to the input matrices ${\bf A}$ and ${\bf B}$.

To clarify the interpretation of $q_{1-\alpha}(t)$, it is helpful to plot the fluctuations of $\varepsilon_{t}$. In the left panel of Figure 1, we illustrate a simulation where randomly generated rows are incrementally added to a sketching matrix ${\bf S}$, with ${\bf A}$ and ${\bf B}$ held fixed. Each time a row is added to ${\bf S}$, the sketch size $t$ increases by 1, and we plot the corresponding value of $\varepsilon_{t}$ as $t$ ranges from 100 to 1,700. (Note that the user is typically unable to observe such a curve in practice.) In the right panel, we display 1,000 repetitions of the simulation, with each colored curve corresponding to one repetition. (The variation is due only to the different draws of ${\bf S}$.) In particular, the function $q_{0.99}(t)$ is represented by the thick black curve, delineating the top 1% of the colored curves at each value of $t$.

In essence, the right panel of Figure 1 shows that if the user had knowledge of the (unknown) function $q_{1-\alpha}(t)$, then two important purposes could be served. First, for any fixed value $t$, the user would have a sharp problem-specific bound on $\varepsilon_{t}$. Second, for any fixed error tolerance $\epsilon$, the user could select $t$ so that that “just enough” computation is spent in order to achieve $\varepsilon_{t}\leq\epsilon$ with probability at least $1-\alpha$.

The challenge we face is that a naive computation of $q_{1-\alpha}(t)$ by generating samples of $\varepsilon_{t}$ would defeat the purpose of sketching. Indeed, generating samples of $\varepsilon_{t}$ by brute force would require running the sketching method many times, and it would also require computing the entire product ${\bf A}^{T}{\bf B}$. Consequently, the technical problem of interest is to develop an efficient way to estimate $q_{1-\alpha}(t)$, without adding much cost to a single run of the sketching method.

3 Contributions

From a conceptual standpoint, the main novelty of our work is that it bridges two sets of ideas that are ordinarily studied in distinct communities. Namely, we apply the statistical technique of bootstrapping to enhance algorithms for numerical linear algebra. To some extent, this pairing of ideas might seem counterintuitive, since bootstrap methods are sometimes labeled as “computationally intensive”, but it will turn out that the cost of bootstrapping can be managed in our context. Another reason our approach is novel is that we use the bootstrap to quantify error in the output of a randomized algorithm, rather than for the usual purpose of quantifying uncertainty arising from data. In this way, our approach harnesses the versatility of bootstrap methods, and we hope that our results in the “use case” of matrix multiplication will encourage broader applications of bootstrap methods in randomized computations. (See also Section 6, and note that in concurrent work, we have pursued similar approaches in the contexts of randomized least-squares and classification algorithms (Lopes et al. 2018b; Lopes 2019).)

From a technical standpoint, our main contributions are a method for estimating the function $q_{1-\alpha}(t)$, as well as theoretical performance guarantees. Computationally, the proposed method is efficient in the sense that its cost is comparable to a single run of standard sketching methods (see Section 2). This efficiency is made possible by an “extrapolation” technique, which allows us to bootstrap small “initial” sketches with $t_{0}$ rows, and inexpensively estimate $q_{1-\alpha}(t)$ at larger values $t\gg t_{0}$. The empirical performance of the extrapolation technique is also quite encouraging, as discussed in Section 5. Lastly, with regard to theoretical analysis, our proofs circumvent some technical restrictions occurring in the analysis of related bootstrap methods in the statistics literature.

4 Related work

Our approach differs from the “test-vector approach” in some essential ways. One difference arises because the bounds on $\hat{\eta}_{t}$ are generally constructed from the vectors $\{{\bf w}_{i}\}$ using conservative inequalities. By contrast, our approach avoids this conservativeness by directly estimating $q_{1-\alpha}(t)$, which is an optimal bound on $\varepsilon_{t}$ in the sense of equation (4).

At a more technical level, the ability to avoid restrictions on ${\bf A}$ and ${\bf B}$ comes from our use of the Lévy-Prohorov metric for distributional approximations — which differs from the Kolmogorov metric that has been predominantly used in previous works on multiplier bootstrap methods. More specifically, analyses based on the Kolmogorov metric typically rely on “anti-concentration inequalities” (Chernozhukov et al. 2013; Chernozhukov et al. 2015), which ultimately lead to the mentioned variance assumptions. On the other hand, our approach based on the Lévy-Prohorov metric does not require the use of anti-concentration inequalities. Finally it should be mentioned that the techniques used to control the LP metric are related to those that have been developed for bootstrap approximations via coupling inequalities as in Chernozhukov et al. 2016.

This paper is organized as follows. Section 2 introduces some technical background. Section 3 describes the proposed bootstrap algorithm. Section 4 establishes the main theoretical results, and then numerical performance is illustrated in Section 5. Lastly, conclusions and extensions of the method are presented in Section 6, and all proofs are given in the appendices.

Preliminaries

We will use $c$ to denote a positive absolute constant that may change from line to line. The matrices ${\bf A}$, ${\bf B}$, and ${\bf S}$ are viewed as lying in a sequence of matrices indexed by the tuple $(d,d^{\prime},t,n)$. For a pair of generic functions $f$ and $g$, we write $f(d,d^{\prime},t,n)\lesssim g(d,d^{\prime},t,n)$ when there is a positive absolute constant $c$ so that $f(d,d^{\prime},t,n)\leq c\,g(d,d^{\prime},t,n)$ holds for all large values of $d,d^{\prime},t,$ and $n$. Furthermore, if $a$ and $b$ are two quantities that satisfy both $a\lesssim b$ and $b\lesssim a$, then we write $a\asymp b$. Lastly, we do not use the symbols $\lesssim$ or $\asymp$ when relating random variables.

Examples of sketching matrices.

Our theoretical results will deal with three common types of sketching matrices, reviewed below.

and leverage score sampling, for which further background may be found in the papers (Drineas et al. 2006b; Drineas et al. 2008; Drineas et al. 2012).

Subsampled randomized Hadamard transform (SRHT). Let $n$ be a power of $2$, and define the Walsh-Hadamard matrix ${\bf H}_{n}$ recursively The restriction that $n$ is a power of 2 can be relaxed with variants of SRHT matrices (Avron et al. 2010; Boutsidis and Gittens 2013).

is called an SRHT matrix. This type of sketching matrix was introduced in the seminal paper (Ailon and Chazelle 2006), and additional details regarding implementation may be found in the papers (Drineas et al. 2011; Wang 2015). (The factor $\textstyle\frac{1}{\sqrt{n}}$ is used so that $\textstyle\frac{1}{\sqrt{n}}{\bf H}_{n}$ is an orthogonal matrix.) An important property of SRHT matrices is that they can be multiplied with any $n\times d$ matrix in ${\mathcal{O}}(n\cdot d\cdot\log t)$ time (Ailon and Liberty 2009), which is faster than the $\mathcal{O}(n\cdot d\cdot t)$ time usually required for a dense sketching matrix.

Methodology

Before presenting our method in algorithmic form, we first explain the underlying intuition.

Furthermore, in the cases of length sampling and Gaussian projection, the matrices ${\bf D}_{1},\dots,{\bf D}_{t}$ are independent, and in the case of SRHT sketches, these matrices are “nearly” independent. So, in light of the central limit theorem, it is natural to suspect that the random matrix (9) will be well-approximated (in distribution) by a matrix with Gaussian entries. In particular, if we examine the $(j_{1},j_{2})$ entry, then we may expect that ${\bf e}_{j_{1}}^{T}\big({\bf A}^{T}{\bf S}^{T}{\bf S}{\bf B}-{\bf A}^{T}{\bf B}\big){\bf e}_{j_{2}}$ will approximately follow the distribution ${\mathcal{N}}(0,\textstyle\frac{1}{t}\sigma_{j_{1},j_{2}}^{2})$, where the unknown parameter $\sigma_{j_{1},j_{2}}^{2}$ can be estimated with

Based on these considerations, the idea of the proposed bootstrap method is to generate a random matrix whose $(j_{1},j_{2})$ entry is sampled from ${\mathcal{N}}(0,\textstyle\frac{1}{t}\hat{\sigma}_{j_{1},j_{2}}^{2})$. It turns out that an efficient way of generating such a matrix is to sample i.i.d. random variables $\xi_{1},\dots,\xi_{t}\sim{\mathcal{N}}(0,1)$, independent of ${\bf S}$, and then compute

then the bootstrap algorithm will generate i.i.d. samples of $\varepsilon_{t}^{\star}$, conditionally on ${\bf S}$. In turn, the $(1-\alpha)$-quantile of the bootstrap samples, say $\hat{q}_{1-\alpha}(t)$, can be used to estimate $q_{1-\alpha}(t)$.

2 Multiplier bootstrap algorithm

3 Saving on computation with extrapolation

In its basic form, the cost of Algorithm 1 is ${\mathcal{O}}(B\cdot t\cdot d\cdot d^{\prime})$, which has the favorable property of being independent of the large dimension $n$. Also, the computation of the samples $\varepsilon^{\star}_{t,1},\dots,\varepsilon_{t,B}^{\star}$ is embarrassingly parallel, with the cost of each sample being ${\mathcal{O}}(t\cdot d\cdot d^{\prime})$. Moreover, due to the way that the quantile $q_{1-\alpha}(t)$ scales with $t$, it is possible to reduce the cost of Algorithm 1 even further — via the technique of extrapolation (also called Richardson extrapolation) (Sidi 2003; Brezinski and Zaglia 2013).

In order to take advantage of the theoretical scaling $q_{1-\alpha}(t)\approx\textstyle\frac{\kappa}{\sqrt{t}}$, we may use Algorithm 1 to compute $\hat{q}_{1-\alpha}(t_{0})$ with an initial sketch size $t_{0}$, and then approximate the value $q_{1-\alpha}(t)$ for $t\gg t_{0}$ with the following extrapolated estimator

Hence, if the user would like to determine a sketch size $t$ so that $q_{1-\alpha}(t)\leq\epsilon$, for some tolerance $\epsilon$, then $t$ should be selected so that $\hat{q}_{1-\alpha}^{\text{ ext}}(t)\leq\epsilon$, which is equivalent to

In our experiments in Section 5, we illustrate some examples where an accurate estimate of $q_{1-\alpha}(t)$ at $t=10,\!000$ can be obtained from the rule (13) using an initial sketch size $t_{0}\approx 500$, yielding a roughly 20-fold speedup on the basic version of Algorithm 1.

Given that the purpose of Algorithm 1 is to enhance sketching methods, it is important to understand how the added cost of the bootstrap compares to the cost of running sketching methods in the standard way. As a point of reference, we compare with the cost of computing ${\bf A}^{T}{\bf S}^{T}{\bf S}{\bf B}$ when ${\bf S}$ is chosen to be an SRHT matrix, since this is one of the most efficient sketching methods. If we temporarily assume for simplicity that ${\bf A}$ and ${\bf B}$ are both of size $n\times d$, then it follows from Section 2 that computing ${\bf A}^{T}{\bf S}^{T}{\bf S}{\bf B}$ has a cost of order $\mathcal{O}(t\cdot d^{2}+n\cdot d\cdot\log(t))$. Meanwhile, the cost of running Algorithm 1 with the extrapolation speedup based on an initial sketch size $t_{0}$ is $\mathcal{O}(B\cdot t_{0}\cdot d^{2})$. Consequently, the extra cost of the bootstrap does not exceed the stated cost of sketching when the number of bootstrap samples satisfies

and in fact, this could be improved further if parallelization of Algorithm 1 is taken into account. It is also important to note that rather small values of $B$ are shown to work well in our experiments, such as $B=20$. Hence, as long $t_{0}$ remains fairly small compared to $t$, then the condition (14) may be expected to hold, and this is borne out in our experiments. The same reasoning also applies when $n\log(t)\gg d\cdot t_{0}$, which conforms with the fact that sketching methods are intended to handle situations where $n$ is very large.

4 Relation with the non-parametric bootstrap

For readers who are more familiar with the “non-parametric bootstrap” (based on sampling with replacement), the purpose of this short subsection is to explain the relationship with the multiplier bootstrap in Algorithm 1. Indeed, an understanding of this relationship may be helpful, since the non-parametric bootstrap might be viewed as more intuitive, and perhaps easier to generalize to more complex situations. However, it turns out that Algorithm 1 is technically more convenient to analyze, and that is why the paper focuses primarily on Algorithm 1. Meanwhile, from a practical point of view, there is little difference between the two approaches, since both have the same order of computational cost, and in our experience, we have observed essentially the same performance in simulations. Also, the extrapolation technique can be applied to both algorithms in the same way.

Main results

Our main results quantify how well the estimate $\hat{q}_{1-\alpha}(t)$ from Algorithm 1 approximates the true value $q_{1-\alpha}(t)$, and this will be done by analyzing how well the distribution of a bootstrap sample $\varepsilon_{t,1}^{\star}$ approximates the distribution of $\varepsilon_{t}$. For the purposes of comparing distributions, we will use the Lévy-Prohorov metric, defined below.

The $d_{\text{LP}}$ metric is a standard tool for comparing distributions, due to the fact that convergence with respect to $d_{\text{LP}}$ is equivalent to convergence in distribution (Huber and Ronchetti 2009, Theorem 2.9).

Approximating quantiles.

An important property of the $d_{\text{LP}}$ metric is that if two distributions are close in this metric, then their quantiles are close in the following sense. Recall that if $F_{U}$ is the distribution function of a random variable $U$, then the $(1-\alpha)$-quantile of $U$ is the same as the generalized inverse $F_{U}^{-1}(1-\alpha):=\inf\{q\in[0,\infty)\,|\,F_{U}(q)\geq 1-\alpha\}$. Next, suppose that two random variables $U$ and $V$ satisfy

for some $\epsilon\in(0,\alpha)$ with $\alpha\in(0,1/2)$. Then, the quantiles of $U$ and $V$ are close in the sense that

where the function $\psi_{\alpha}(\epsilon):=F_{U}^{-1}(1-\alpha+\epsilon)-F_{U}^{-1}(1-\alpha-\epsilon)+\epsilon$ is strictly monotone, and satisfies $\psi_{\alpha}(0)=0$. (For a proof, see Lemma 15 of Appendix F.) In light of this fact, it will be more convenient to express our results for approximating $q_{1-\alpha}(t)$ in terms of the $d_{\text{LP}}$ metric.

1 Statements of results

Our main assumption involves three separate cases, corresponding to different choices of the sketching matrix ${\bf S}$.

The dimensions $d$ and $d^{\prime}$ satisfy $d\asymp d^{\prime}$. Also, there is a positive absolute constant $\kappa\geq 1$ such that $d^{1/\kappa}\lesssim t\lesssim d^{\kappa}$, which is to say that neither $d$ nor $t$ grows exponentially with the other. In addition, one of the following sets of conditions holds, involving the parameter $\nu({\bf A},{\bf B}):=\sqrt{\|{\bf A}^{T}{\bf A}\|_{\infty}\|{\bf B}^{T}{\bf B}\|_{\infty}}$.

(Length sampling case). The matrix ${\bf S}$ is generated by length sampling, with the probabilities in equation (5), and also, $t\;\gtrsim\;(\|{\bf A}\|_{F}\|{\bf B}\|_{F})^{2/3}(\log d)^{5}$.

(SRHT case). The matrix ${\bf S}$ is an SRHT matrix as defined in equation (6), and also, $t\gtrsim\nu({\bf A},{\bf B})^{2/3}(\log n)^{2}(\log d)^{5}$.

With regard to the original problem of estimating the quantile $q_{1-\alpha}(t)$ for $\varepsilon_{t}$, this rescaling makes no essential difference, since quantiles are homogenous with respect to scaling, and in particular, the $(1-\alpha)$-quantile of $Z_{t}$ is simply $\sqrt{t}q_{1-\alpha}(t)$.

As a second clarification, recall that the bootstrap method generates samples $\varepsilon_{t}^{\star}$ based upon a particular realization of ${\bf S}$. For this reason, the bootstrap approximation to $\mathcal{L}(Z_{t})$ is the conditional distribution $\mathcal{L}(Z_{t}^{\star}|{\bf S})$. Consequently, it should be noted that $\mathcal{L}(Z_{t}^{\star}|{\bf S})$ is a random probability measure, and $d_{\text{LP}}(\mathcal{L}(Z_{t})\,,\>\mathcal{L}(Z_{t}^{\star}|{\bf S}))$ is a random variable, since they both depend on the random matrix ${\bf S}$.

Let $h(x)=x^{1/2}+x^{3/4}$ for $x\geq 0$. If Assumption 1 (a) holds, then there is an absolute constant $c>0$ such that the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$,

If Assumption 1 (b) holds, then there is an absolute constant $c>0$ such that the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$,

Remarks.

A noteworthy property of the bounds is that they are dimension-free with respect to the large dimension $n$. Also, they have a very mild logarithmic dependence on $d$. With regard to the dependence on $t$, there are two other important factors to keep in mind. First, the practical performance of the bootstrap method (shown in Section 5) is much better than what the $t^{-1/8}$ rate suggests. Second, the problem of finding the optimal rates of approximation for multiplier bootstrap methods is a largely open problem — even in the simpler setting of bootstrapping the coordinate-wise maximum of vectors (rather than matrices). In the vector context, the literature has focused primarily on the Kolmogorov metric (rather than the LP metric), and some quite recent improvements beyond the $t^{-1/8}$ rate have been developed in Chernozhukov et al. 2017 and Lopes et al. 2018a. However, these works also rely on model assumptions that would lead to additional restrictions on the matrices $\mathbf{A}$ and $\mathbf{B}$ in our setup. Likewise, the problem of extending our results to achieve faster rates or handle other metrics is a natural direction for future work.

The SRHT case.

For the case of SRHT matrices, the analogue of Theorem 1 needs to be stated in a slightly different way for technical reasons. From a qualitative standpoint, the results for SRHT and sub-Gaussian matrices turn out to be similar.

Let $h(x)=x^{1/2}+x^{3/4}$ for $x\geq 0$. If Assumption 1 (c) holds, then there is an absolute constant $c>0$ such that the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}-\textstyle\frac{c}{n}$,

Remarks.

Up to a factor involving $\log(n)$, the bound for SRHT matrices matches that for sub-Gaussian matrices. Meanwhile, from a more practical standpoint, our empirical results will show that the bootstrap’s performance for SRHT matrices is generally similar to that for both sub-Gaussian and length-sampling matrices.

Further discussion of results.

To comment on the role of $\nu({\bf A},{\bf B})$ and $\|{\bf A}\|_{F}\|{\bf B}\|_{F}$ in Theorems 1 and 2, it is possible to interpret them as problem-specific “scale parameters”. Indeed, it is natural that the bounds on $d_{\text{LP}}$ should increase with the scale of ${\bf A}$ and ${\bf B}$ for the following reason. Namely, if ${\bf A}$ or ${\bf B}$ is multiplied by a scale factor $\kappa>0$, then it can be checked that the quantile error $|\hat{q}_{1-\alpha}(t)-q_{1-\alpha}(t)|$ will also change by a factor of $\kappa$, and furthermore, the inequality (15) demonstrates a monotone relationship between the sizes of the quantile error and the $d_{\text{LP}}$ error. For this reason, the bootstrap may still perform well in relation to the scale of the problem when the magnitudes of the parameters $\nu({\bf A},{\bf B})$ and $\|{\bf A}\|_{F}\|{\bf B}\|_{F}$ are large. Alternatively, this idea can be seen by noting that the $d_{\text{LP}}$ bounds can be made arbitrarily small by simply changing the units used to measure the entries of ${\bf A}$ and ${\bf B}$.

Beyond these considerations, it is still of interest to compare the results for different sketching matrices once a particular scaling has been fixed. For concreteness, consider a scaling where the spectral norms of ${\bf A}$ and ${\bf B}$ satisfy $\|\mathbf{A}\|_{2}\asymp\|\mathbf{B}\|_{2}\asymp 1$. (As an example, if we view $\mathbf{A}^{\top}\mathbf{A}$ as a sample covariance matrix, then the condition $\|\mathbf{A}\|_{2}\asymp 1$ simply means that the largest principal component score is of order 1.) Under this scaling, it is simple to check that $\nu(\mathbf{A},\mathbf{B})=\mathcal{O}(1)$, and $\|\mathbf{A}\|_{F}\|\mathbf{B}\|_{F}=\mathcal{O}(\sqrt{\mathsf{r}(\mathbf{A})\mathsf{r}(\mathbf{B})})$, where $\mathsf{r}(\mathbf{A}):=\|\mathbf{A}\|_{F}^{2}/\|\mathbf{A}\|_{2}^{2}$ is the “stable rank”. In particular, note that if $\mathbf{A}$ and $\mathbf{B}$ are approximately low rank, as is common in applications, then $\mathsf{r}(\mathbf{A})\ll d$, and $\mathsf{r}(\mathbf{B})\ll d^{\prime}$. Accordingly, we may conclude that if the conditions of Theorems 1 and 2 hold, then bootstrap consistency occurs under the following limits

where we have used the simplifying assumption that $d\asymp d^{\prime}$.

Experiments

This section outlines a set of experiments for evaluating the performance of Algorithm 1 with the extrapolation speed-up described in Section 3.3. The experiments involved both synthetic and natural matrices, as described below.

Natural matrices.

We also conducted experiments on five natural data matrices ${\bf A}$ from the LIBSVM repository Chang and Lin 2011, named ‘Connect’, ‘DNA’, ‘MNIST’, ‘Mushrooms’, and ‘Protein’, with the same normalization that was used for the synthetic matrices. These datasets are briefly summarized in Table 1.

1 Design of experiments

Extrapolated estimates.

In order to illustrate the variability of the estimate $\hat{q}_{0.99}^{\,\text{ext}}(t)$ over the 1,000 realizations, we plot three different curves as a function of $t$. The blue curve represents the average value of $\hat{q}_{0.99}^{\,\text{ext}}(t)$, while the green and yellow curves respectively correspond to the estimates ranking 100th an 900th out of the 1,000 realizations.

2 Comments on numerical results

With attention to the extrapolation rule (12), there are two main points to note. First, the plots show that the extrapolation may be initiated at fairly low values of $t_{0}$, which are much less than the sketch sizes needed to achieve a small sketching error $\varepsilon_{t}$. Second, we see that $\hat{q}_{0.99}^{\,\text{ext}}(t)$ remains accurate for $t$ much larger than $t_{0}$, well up to $t=10,000$ and perhaps even farther. Consequently, the results show that the extrapolation technique is capable of saving quite a bit of computation without much detriment to statistical performance.

To consider the relationship between theory and practice, one basic observation is that all three types of sketching matrices obey roughly similar bounds in Theorems 1 and 2, and indeed, we also see generally similar numerical performance among the three types. At a more fine-grained level however, the Gaussian and SRHT sketching matrices tend to produce estimates $\hat{q}^{\,\text{ext}}_{0.99}(t)$ with somewhat higher variance than in the case of length sampling. Another difference between theory and simulation, is that the actual performance of the method seems to be better than what the theory suggests — since the estimates are accurate at values of $t_{0}$ that are much smaller than what would be expected from the rates in Theorems 1 and 2.

Conclusions and extensions

In this paper, we have focused on estimating the quantile $q_{1-\alpha}(t)$ as a way of addressing two fundamental issues in randomized matrix multiplication: (1) knowing how accurate a given sketched product is, and (2) knowing how much computation is needed to achieve a specified degree of accuracy. With regard to methodology, our approach is relatively novel in that it uses the statistical technique of bootstrapping to serve a computational purpose — by quantifying the error of a randomized sketching algorithm. A second important component of our method is the extrapolation technique, which ensures that the cost of estimating $q_{1-\alpha}(t)$ does not substantially increase the overall cost of standard sketching methods. Furthermore, our numerical results show that the extrapolated estimate is quite accurate in a variety of different situations, suggesting that our method may offer a general way to enhance sketching algorithms in practice.

More generally, the problems we have addressed for randomized matrix multiplication arise for many other large-scale matrix computations. Hence, it is natural to consider extensions of our approach to more complex settings, and in the remainder of this section, we briefly mention a few possibilities for future study.

At a high level, each of the applications below deals with an object, say $\Theta$, that is difficult to compute, as well as a randomized approximation, say $\widetilde{\Theta}$, that is built from a sketching matrix ${\bf S}$ with $t$ rows. Next, if we consider the random error variable

which has ${\mathcal{O}}(nd^{2}+ndd^{\prime})$ cost. In the case where $\max\{d,d^{\prime}\}\ll n$, the matrix multiplications are a computational bottleneck, and an approximate solution can be obtained via

Appendices

Outline of appendices. Appendix A explains the main conceptual ideas underlying the proofs of Theorems 1 and 2. In particular, the proofs of these theorems will be decomposed into two main results: Propositions 3 and 4, which are given in Appendix A.

Appendix B will prove the sub-Gaussian case of Proposition 3, and Appendix C will prove the sub-Gaussian case of Proposition 4. Later on, Appendices D and E, will explain how the arguments can be changed to handle the length-sampling and SRHT cases.

Conventions used in proofs. If either of the matrices ${\bf A}$ or ${\bf B}$ are ${\bf 0}$, then $\varepsilon_{t}$ has a trivial point-mass distribution at 0. In this degenerate case, it is simple to check that the bootstrap produces an exact approximation. So, without loss of generality, all proofs are written under the assumption that ${\bf A}$ and ${\bf B}$ are non-zero. Next, since Assumption 1 is formulated using the $\lesssim$ notation, there is no loss of generality in carrying out calculations under the assumption that all the numbers $t,n,d,d^{\prime}$ are at least 8, which will ensure that quantities such as $\log(d)$ are greater than 2. Lastly, if a numbered lemma is invoked in the middle of a proof, the lemma may be found in Appendix F.

Appendix A Gaussian and bootstrap approximations

Section A.1 introduces some notation that helps us to analyze the rescaled sketching error $Z_{t}$ from the viewpoint of empirical processes. Next, in Section A.2, Theorem 1 will be decomposed into two propositions that compare $Z_{t}$ and $Z_{t}^{\star}$ with the maximum of a suitable Gaussian process. The proofs of these propositions may be found in Appendices B and C.

The main idea of our analysis is to view $Z_{t}$ as the maximum of an empirical process, which we now define. Recall the notation

For future reference, we also define the corresponding bootstrap process

where $\xi_{1},\dots,\xi_{t}$ are i.i.d. ${\mathcal{N}}(0,1)$ and independent of ${\bf S}$.

Clearly, $\mathsf{card}({\mathscr{F}})=2dd^{\prime}$. Under this definition, it is simple to check that $Z_{t}$ and $Z_{t}^{\star}$, defined in equations (16) and (17), can be expressed as

A.2 Statements of the approximation results

for all ${\bf j},{\bf k}\in{\mathscr{J}}$. In turn, define the following random variable as the the maximum of this Gaussian process,

We are now in position to state the approximation results.

Under Assumption 1 (a), the following bound holds,

Under Assumption 1 (b), the following bound holds,

Under Assumption 1 (c), the following bound holds with probability at least $1-c/n$

If Assumption 1 (a) holds, then the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$,

If Assumption 1 (b) holds, then the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$,

If Assumption 1 (c) holds, then the following bound holds with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}-\textstyle\frac{c}{n}$,

Appendix B Proof of Proposition 3, part (a)

where we define the following non-random quantities

The remainder of the proof consists in bounding each of these quantities, and we will establish the following two bounds for all $\delta>0$,

Recall also that $\text{card}({\mathscr{F}})=2dd^{\prime}$, and $d\asymp d^{\prime}$ under Assumption 1.

For the moment, we set aside the task of proving these bounds, and consider the choice of $\delta$. There are two constraints that we would like $\delta$ to satisfy. First, we would like to choose $\delta$ so that the bounds on $L_{t}$ and $(K_{t}(\delta)+J_{t}(\delta))$ are of the same order. In particular, we desire

Second, with regard to line (24) we would like $\delta$ to solve the equation

so that the second term in line (24) is of order $\delta$. The idea is that if $\delta$ satisfies both of the conditions (30) and (31), then the definition of the $d_{\text{LP}}$ metric and line (24) imply

which clearly satisfies line (31). Futhermore, it can be checked that $\delta_{0}$ also satisfies the constraint (30) under Assumption 1 (a). (The details of verifying this are somewhat tedious and are given in Lemma 16 in Appendix F.)

To finish the proof, it remains to establish the bounds (28) and (29). To handle $L_{t}$, note that In this step, we use the assumption that $\|\sqrt{t}S_{i,j}\|_{\psi_{2}}\leq c$ for all $i$ and $j$.

which proves the claimed bound in line (28).

Next, regarding $K_{t}(\delta)$, let us consider the random variable

It follows from Lemma 9 (part 4) and Lemma 13 in Appendix F that $K_{t}(\delta)$ can be bounded in terms of the Orlicz norm $\|\eta\|_{\psi_{1}}$,

To handle $\|\eta\|_{\psi_{1}}$, it follows from Lemma 9 (part 3), that

Furthermore, due to the earlier calculation starting at line (32) above,

Combining the last few steps, we conclude that

Lastly, we turn to bounding $J_{t}(\delta)$. Fortunately, much of the argument for bounding $K_{t}(\delta)$ can be carried over. Specifically, consider the random variable

Lemma 13 in Appendix F shows that $J_{t}(\delta)$ can be bounded in terms of $\|\zeta\|_{\psi_{1}}$,

Proceeding in a way that is similar to the bound for $K_{t}(\delta)$, it follows from part (3) of Lemma 9 that

Furthermore, for every $f_{{\bf j}}\in{\mathscr{F}}$, the facts in Lemma 9 imply

Appendix C Proof of Proposition 4, part (a)

If we set $\delta$ to the particular choice $\delta_{0}:=t^{-1/8}\sqrt{\nu({\bf A},{\bf B})\cdot\,\log(\mathsf{card}({\mathscr{F}}))}$, then $\delta_{0}$ solves the equation

Consequently, by the definition of the $d_{\text{LP}}$ metric, this implies that whenever the event $\{{\bf S}\in{\mathscr{S}}\}$ occurs, we have

and this implies the statement of Proposition 4.

where we define the following function of ${\bf S}$,

Based on this definition, it is simple to check that the proof is reduced to showing that the event $\{{\bf S}\in{\mathscr{S}}\}$ occurs with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$. This is guaranteed by the lemma below.

Suppose Assumption 1 (a) holds. Then, the event

occurs with probability at least $1-\textstyle\frac{1}{t}-\textstyle\frac{1}{dd^{\prime}}$.

We begin by bounding $\Delta_{t}({\bf S})$ with two other quantities (to be denoted $\Delta_{t}^{\prime}({\bf S})$, $\Delta_{t}^{\prime\prime}({\bf S})$) that are easier to bound. Using the fact that $\textstyle\frac{1}{t}\sum_{i=1}^{t}{\bf D}_{i}={\bf A}^{T}{\bf S}^{T}{\bf S}{\bf B}$ it can be checked that

From looking at the last two lines, it is natural to define the following zero-mean random variables for any triple $i,{\bf j},{\bf k}$, Note that $Q_{i,{\bf j},{\bf k}}$ is a multivariate polynomial of degree-4 in the variables $S_{i,j}$, and so techniques based on moment generating functions, like Chernoff bounds, are not generally applicable to controlling $Q_{i,{\bf j},{\bf k}}$. For instance, if $X\sim{\mathcal{N}}(0,1)$, then the variable $X^{4}$ does not have a moment generating function. Handling this obstacle is a notable aspect of our analysis.

Suppose Assumption 1 (a) holds. Then, the event

occurs with probability at least $1-\frac{1}{t}$, and the event

occurs with probability at least $1-\frac{1}{dd^{\prime}}$.

Let $p>2$. Due to part (3) of Lemma 9 in Appendix F, we have

Note that each variable $Q_{i,{\bf j},{\bf k}}$ has moments of all orders, and when ${\bf j}$ and ${\bf k}$ are held fixed, the sequence $\{Q_{i,{\bf j},{\bf k}}\}_{1\leq i\leq t}$ is i.i.d. For this reason, it is natural to use Rosenthal’s inequality to bound the $L_{p}$ norm of the right side of the previous line. Specifically, the version of Rosenthal’s inequality Here we are using the version of Rosenthal’s inequality with the optimal dependence on $p$. It is a notable aspect of our argument that it makes essential use of this scaling in $p$. stated in Lemma 10 in Appendix F leads to

The $L_{2}$ norm on the right side of Rosenthal’s inequality (45) satisfies the bound

where the last step follows from the fact

obtained in the bounds (32) through (33).

Next, to handle the $L_{p}$ norms in the bound (45), observe that

Hence, the second term in the Rosenthal bound (45) satisfies

and as long as the first term in the Rosenthal bound dominates Under the choice of $p=\log(\mathsf{card}({\mathscr{F}}))=\log(2dd^{\prime})$ that will be made at the end of this argument, it is straightforward to check that the condition (47) holds under Assumption 1., i.e.

then we conclude that for any ${\bf j}$ and ${\bf k}$,

Since the previous bound does not depend on ${\bf j}$ or ${\bf k}$, combining it with the first step in line (44) leads to

Next, we convert this norm bound into a tail bound. Specifically, if we consider the value

and noting that $\mathsf{card}({\mathscr{F}})^{1/p}=e$, it follows that under this choice of $p$,

Moreover, as long as $t\lesssim\mathsf{card}({\mathscr{F}})^{\kappa}$ for some absolute constant $\kappa\geq 1$ (which holds under Assumption 1), then the last factor on the right satisfies

So, combining the last few steps, there is an absolute constant $c$ such that

Proof of Lemma 6 (ii).

Note that for each $i\in[t]$ and ${\bf j}\in{\mathscr{J}}$, we have

which is a centered sub-Gaussian quadratic form. Due to the bound (35), we have

Furthermore, this can be combined with a standard concentration bound for sums of independent sub-exponential random variables (Lemma 12) to show that for any $r\geq 0$,

Hence, taking a union bound over all ${\bf j}$ gives

Regarding the choice of $r$, note that by Assumption 1, we have $\textstyle\frac{1}{\sqrt{t}}\sqrt{\log(\mathsf{card}({\mathscr{F}}))}\lesssim 1$. It follows that there is a sufficiently large absolute constant $c_{1}>0$ such that if we put

where $c$ is the same as in the bound (51). In turn, this implies

Appendix D Proof of Propositions 3 and 4 in case (b) (length sampling)

Hence, it remains to show that $\|f_{{\bf j}}({\bf D}_{1})\|_{\psi_{1}}\leq c\|{\bf A}\|_{F}\|{\bf B}\|_{F}$, which is the content of Lemma 7 below.

If ${\bf S}$ is generated by length sampling with the probabilities in line (5), then for any ${\bf j}\in\mathscr{J}$, we have the bound

Hence, if we take $r=2\|{\bf A}\|_{F}\|{\bf B}\|_{F}$, then the right hand side is at most $e^{1/2}\leq 2$.

Appendix E Proof of Propositions 3 and 4 in case (c) (SRHT)

Since we are not aware of a standard notation for a conditional Orlicz norm, we define

which is a random variable, since it is a function of ${\bf D}_{n}^{\circ}$. The following lemma provides a bound on this quantity, which turns out to be of order $\log(n)\,\nu({\bf A},{\bf B})$. For this reason, the SRHT case (c) of Propositions 3 and 4 will have the same form as case (a), but with $\log(n)\,\nu({\bf A},{\bf B})$ replacing $\nu({\bf A},{\bf B})$.

If ${\bf S}$ is an SRHT matrix, then the following bound holds with probability at least $1-c/n$,

For an SRHT matrix ${\bf S}={\bf P}\textstyle\frac{1}{\sqrt{n}}{\bf H}_{n}{\bf D}_{n}^{\circ}$, recall that the rows of $\sqrt{t}{\bf P}$ are sampled uniformly at random from the set $\{\textstyle\frac{1}{\sqrt{1/n}}{\bf e}_{1},\dots,\textstyle\frac{1}{\sqrt{1/n}}{\bf e}_{n}\}$. It follows that

Finally, this means that if we take $r=4c\log(n)\nu({\bf A},{\bf B})$ in the bound (56), then the event

holds with probability at least $1-\textstyle\frac{c}{n}$, which completes the proof, since $e^{1/2}\leq 2$. ∎

Appendix F Technical Lemmas

Orlicz norms have the following properties, where $c,c_{1},$ and $c_{2}$ are positive absolute constants.

For any random variable $X$, and any $p\geq 1$,

If $X\sim{\mathcal{N}}(0,\sigma^{2})$, then $\|X\|_{\psi_{2}}\leq c\,\sigma$.

Let $p\geq 1$. For any sequence of random variables $X_{1},\dots,X_{d}$,

Let $X$ be any random variable. Then, for any $x>0$ and $p\geq 1$, we have

In part 1, line (59) follows from line 5.11 of Vershynin 2012, line (60) follows from definition 5.13 of Vershynin 2012, and line (61) follows from p.94 of van der Vaart and Wellner 1996. Next, part 2 follows from the definition of the $\psi_{2}$-Orlicz norm and the moment generating function for ${\mathcal{N}}(0,\sigma^{2})$. Part 3 is due to Lemma 2.2.2 of van der Vaart and Wellner 1996. Lastly, part 4 follows from Markov’s inequality and line 5.14 of Vershynin 2012. ∎

See the paper Johnson et al. 1985. The statement above differs slightly from the Theorem 4.1 in the paper Johnson et al. 1985, which requires symmetric random variables, but the remark on p.247 of that paper explains why the variables $Y_{1},\dots,Y_{t}$ need not be symmetric as long as they have mean 0. ∎

See the paper Rudelson and Vershynin 2013. ∎

If $\eta$ is a non-negative random variable, and there are numbers $a,b>0$ such that

for all $x>0$, then the following bound holds for all $r>0$,

See Lemma 6.6 in Chernozhukov et al. 2016. ∎

The following lemma may be of independent interest, since it provides an explicit bound on the $\psi_{1}$-Orlicz norm of a centered sub-Gaussian quadratic form. Although this bound follows from the Hanson-Wright inequality, we have not seen it stated in the literature.

Next, we employ the Hanson-Wright inequality (Lemma 11). By considering the “threshold” $u^{*}:=\kappa^{2}\textstyle\frac{\|{\bf H}\|_{F}^{2}}{\|{\bf H}\|_{2}}$, it is helpful to note that the quantities in the exponent of the Hanson-Wright inequality satisfy $\textstyle\frac{u^{2}}{\kappa^{4}\|{\bf H}\|_{F}^{2}}\leq\textstyle\frac{u}{\kappa^{2}\|{\bf H}\|_{2}}$ if and only if $u\leq u^{*}$. Hence,

Evaluating the last two integrals directly, if we let $C^{\prime}:=\textstyle\frac{c}{\kappa^{2}\|{\bf H}\|_{2}}-\textstyle\frac{1}{r}$ and choose $r$ so that $C^{\prime}>0$, then

Note that the condition $C^{\prime}>0$ means that it is necessary to have $r>\textstyle\frac{1}{c}\kappa^{2}\|{\bf H}\|_{2}$. To finish the argument, we further require that $r$ is large enough so that (say)

as desired. Note that the constraints (65) are the same as

Remark.

Fix $\alpha\in(0,1/2)$ and suppose there is some $\epsilon\in(0,\alpha)$ such that random variables $U$ and $V$ satisfy

Then, the quantiles of $U$ and $V$ satisfy

It is a fact that this metric is always dominated by the $d_{\text{LP}}$ metric in the sense that

for all scalar random variables $U$ and $V$ (Huber and Ronchetti 2009, p.36). Based on the definition of the $d_{\text{L}}$ metric, it is straightforward to check that the following inequalities hold under the assumption of the lemma,

(Specifically, consider the choices $x=F_{V}^{-1}(1-\alpha)$ and $x=F_{U}^{-1}(1-\alpha+\epsilon)+\epsilon$.) Next, if we subtract $F_{U}^{-1}(1-\alpha)$ from each side of the inequalities above, and note that $F_{U}^{-1}(\cdot)$ is non-decreasing, it follows that if we put $a=F_{U}^{-1}(1-\alpha+\epsilon)-F_{U}^{-1}(1-\alpha)+\epsilon$ and $b=F_{U}^{-1}(1-\alpha)-F_{U}^{-1}(1-\alpha-\epsilon)+\epsilon$, then

where $r\geq 1$ is a free parameter to be adjusted. Based on the bound (29), it is easy to check that plugging $\delta_{1}(r)$ into $K_{t}(\cdot)$ and $J_{t}(\cdot)$ leads to

and if we take $r\geq c\log(\log(d)^{4})$, then

as desired in (30). Hence, as long as there is a choice of $r$ satisfying

then $\delta_{1}(r)$ will satisfy both of the desired constraints (30) and (31). Solving the equation $\delta_{1}(r)=\delta_{0}$ gives

and then the condition $r\geq c\log(\log(d)^{4})$ is the same as

References