Error Estimation for Randomized Least-Squares Algorithms via the Bootstrap

Introduction

Randomized sketching algorithms have been intensively studied in recent years as a general approach to computing fast approximate solutions to large-scale least-squares (LS) problems (Drineas et al. 2006, Rokhlin and Tygert 2008, Avron et al. 2010, Drineas et al. 2011, Mahoney 2011, Drineas et al. 2012, Clarkson and Woodruff 2013, Woodruff 2014, Ma et al. 2014, Meng et al. 2014, Pilanci and Wainwright 2015, Pilanci and Wainwright 2016). During this time, much progress has been made in analyzing the performance of these algorithms, and existing theory provides a good qualitative description of approximation error (relative to the exact solution) in terms of various problem parameters. However, in practice, the user rarely knows the actual error of a randomized solution, or how much extra computation may be needed to achieve a desired level of accuracy.

A basic source of this problem is that it is difficult to translate theoretical error bounds into numerical error bounds that are tight enough to be quantitatively meaningful. For instance, theoretical bounds are often formulated to hold for the worst-case input among a large class of possible inputs. Consequently, they are often pessimistic for “generic” problems, and they may not account for the structure that is unique to the input at hand. Another practical issue is that these bounds typically involve constants that are either conservative, unspecified, or dependent on unknown parameters.

In contrast with worst-case error bounds, we are interested in “a posteriori” error estimates. By this, we mean error bounds that can be estimated numerically in terms of the computed solution or other observable information. Although methods for obtaining a posteriori error estimates are well-developed in some areas of computer science and applied mathematics, there has been very little development for randomized sketching algorithms (cf. Section 1.4). (For brevity, we will usually omit the qualifier ‘a posteriori’ from now on when referring to error estimation.)

The main purpose of this paper is to show that it is possible to directly estimate the error of randomized LS solutions in a way that is both practical and theoretically justified. Accordingly, we propose a flexible estimation method that can enhance existing sketching algorithms in a variety of ways. In particular, we will explain how error estimation can help the user to (1) select the “sketch size” parameter, (2) assess the convergence of iterative sketching algorithms, and (3) measure error in a wider range of metrics than can be handled by existing theory.

Classic Sketch (CS). For a given sketching matrix $S$, this type of algorithm produces a solution

and chronologically, this was the first type of sketching algorithm for LS (Drineas et al. 2006).

Hessian Sketch (HS). The HS algorithm modifies the objective function in the problem (1) so that its Hessian is easier to compute (Pilanci and Wainwright 2016, Becker et al. 2017), leading to a solution

This algorithm is also called “partial sketching”.

Remark. If the initial point for IHS is chosen as $\hat{x}_{0}=0$, then the first iterate $\hat{x}_{1}$ is equivalent to the HS solution $\breve{x}$ in equation (3). Consequently, HS may be viewed as a special case of IHS, and so we will restrict our discussion to CS and IHS for simplicity.

To briefly review the computational benefits of sketching algorithms, first recall that the cost of solving the full least-squares problem (1) by standard methods is $\mathcal{O}(nd^{2})$ (Golub and Van Loan 2012). On the other hand, if the cost of computing the matrix product $SA$ is denoted $C_{\text{sketch}}$, and if a standard method is used to solve the sketched problem (2), then the total cost of CS is $\mathcal{O}(md^{2}+C_{\text{sketch}})$. Similarly, the total cost of IHS with $t$ iterations is $\mathcal{O}(t(md^{2}+C_{\text{sketch}}))$. Regarding the sketching cost $C_{\text{sketch}}$, it depends substantially on the choice of $S$, but there are many types that improve upon the naive $\mathcal{O}(mnd)$ cost of unstructured matrix multiplication. For instance, if $S$ is chosen to be a Sub-sampled Randomized Hadamard Transform (SRHT), then $C_{\text{sketch}}=\mathcal{O}(nd\log(m))$ (Ailon and Chazelle 2006, Sarlós 2006, Ailon and Liberty 2009). Based on these considerations, sketching algorithms can be more efficient than traditional LS algorithms when $md^{2}+nd\log(m)\ll nd^{2}$.

2 Problem Formulation

3 Main Contributions

At a high level, a distinguishing feature of our approach is that it applies inferential ideas from statistics in order to enhance large-scale computations. To be more specific, the novelty of this approach is that it differs from the traditional framework of using bootstrap methods to quantify uncertainty arising from data (Davison and Hinkley 1997). Instead, we are using these methods to quantify uncertainty in the outputs of randomized algorithms — and there do not seem to be many works that have looked at the bootstrap from this angle. From a more theoretical standpoint, another main contribution is that we offer the first guarantees for a posteriori error estimation involving the CS and IHS algorithms. (As a clarification, it should be noted that these results appeared in a conference version of the current work (Lopes et al. 2018), but the proofs here have not previously been published.)

Looking beyond the present setting, there may be further opportunities for using bootstrap methods to estimate the errors of other randomized algorithms. In concurrent work, we have taken this approach in the distinct settings of randomized matrix multiplication, and randomized ensemble classifiers (Lopes et al. 2017, Lopes 2018).

4 Related work

Method

The proposed bootstrap method is outlined in Sections 2.1 and 2.2 for the cases of CS and IHS respectively. The formal analysis can be found in the proof of Theorem 1 in the appendices. Later on, in Section 2.3, we discuss computational cost and speedups.

2 Error Estimation for IHS

At first sight, it might seem that applying the bootstrap to IHS would be substantially different than in the case of CS — given that IHS is an iterative algorithm, whereas CS is a “one-shot” algorithm. However, the bootstrap only needs to be modified slightly. Furthermore, the bootstrap relies on just the final two iterations of a single run of IHS.

Remark. The ideas underlying the IHS version of the bootstrap are broadly similar to those discussed for the CS version. However, the details of this argument are much more involved than in the CS case, owing to the iterative nature of IHS. A formal analysis may be found in the proof of Theorem 1 in the appendices.

3 Computational Cost and Speedups

Cost of error estimation is independent of $n$. The inputs to Algorithms 1 and 2 consist of pre-computed matrices of size $m\times d$, or pre-computed vectors of dimension $d$. Consequently, both algorithms are highly scalable in the sense that their costs do not depend on the large dimension $n$. As a point of comparison, it should be noted that sketching algorithms for LS generally have costs that scale linearly with $n$.

4 Extrapolating with respect to mm for CS

5 Extrapolating with respect to tt for IHS

where $c>0$ and $\eta\in(0,1)$ are unknown parameters that do not depend on $i$.

The simple form of this bound lends itself to extrapolation. Namely, if estimates $\hat{c}$ and $\hat{\eta}$ can be obtained after the first 2 iterations of IHS, then the user can construct the extrapolated error estimate

which predicts how the error will decrease at all subsequent iterations $i\geq 3$. As a result, the user can adaptively determine how many extra iterations (if any) are needed for a specified error tolerance. Furthermore, with the help of Algorithm 2, it is straightforward to estimate $c$ and $\eta$. Indeed, from looking at the condition (19), we desire estimates $\hat{c}$ and $\hat{\eta}$ that solve the two equations

and direct inspection shows that the choices $\hat{\eta}:=\textstyle\frac{\hat{\varepsilon}_{2}(\alpha)}{\hat{\varepsilon}_{1}(\alpha)}$ and $\hat{c}:=\textstyle\frac{\hat{\varepsilon}_{1}(\alpha)}{\hat{\eta}}$ serve this purpose. In Section 4, our experiments show that this simple extrapolation procedure works remarkably well.

Main Result

Since sketching algorithms are most commonly used when $d\ll n$, our results will treat $d$ as fixed while $n\to\infty$. Also, the sketch size $m$ is often selected as a large multiple of $d$, and so we treat $m=m(n)$ as diverging simultaneously with $n$. However, we make no restriction on the size of the ratio $m/n$, which may tend to 0 at any rate. In the same way, the number of bootstrap samples $B=B(n)$ is assumed to diverge with $n$, and the ratio $B/n$ may tend to 0 at any rate. With regard to the number of iterations $t$, its dependence on $n$ is completely unrestricted, and $t=t(n)$ is allowed to remain fixed or diverge with $n$. (The fixed case with $t=1$ is of interest since it describes the HS algorithm.) Apart from these scaling conditions, we use the following two assumptions on $A$ and $b$, as well as the sketching matrices.

In essence, this assumption ensures that the sequence of LS problems is “asymptotically stable”, in the sense that the optimal solution $x_{\textup{opt}}$ does not change erratically from $n$ to $n+1$.

Remarks. Although this result can be stated in a concise form, the proof is actually quite involved. Perhaps the most significant technical obstacle is the sequential nature of the IHS algorithm. To handle the dependence of $\hat{x}_{t}$ on the previous iterates, it is natural to analyze $\hat{x}_{t}$ conditionally on them. However, because the set of previous iterates can grow with $n$, it seems necessary to establish distributional limits that hold “uniformly” over those iterates — and this need for uniformity creates difficulties when applying standard arguments.

Experiments

In this section, we present experimental results in the contexts of CS and IHS. At a high level, there are two main takeaways: (1) The extrapolation rules accurately predict how estimation error depends on $m$ or $t$, and this is shown in a range of conditions. (2) In all of the experiments, the algorithms are implemented with only $B=20$ bootstrap samples. The fact that favorable results can be obtained with so few samples underscores the point that the method incurs only modest cost in exchange for an accuracy guarantee.

Our numerical results are based on four linear regression datasets; two natural, and two synthetic. The natural datasets ‘YearPredictionMSD’, $n=463,\!715$, $d=90$, abbrev. MSD), and ‘cpusmall’ ($n=8,192$, $d=12$, abbrev. CPU) are available at the LIBSVM repository (Chang and Lin 2011).

2 Experiments for CS.

Comments on results for IHS. At a glance, Figure 2 shows that the extrapolated estimate stays on track with the ideal benchmark, and is a nearly unbiased estimate of $\varepsilon_{\text{IHS},i}(.05)$, for $i=3,\dots,10$. An interesting feature of the plots is how much the convergence rate of IHS depends on $m$. Specifically, we see that after 10 iterations, the choice of $m=10d$ versus $m=50d$ can lead to a difference in accuracy that is 4 or 5 orders of magnitude. This sensitivity to $m$ illustrates why selecting $t$ is a non-trivial issue in practice, and why the extrapolated estimate can provide a valuable source of extra information.

Conclusion

We have proposed a systematic approach to answer a very practical question that arises for randomized LS algorithms: “How accurate is a given solution?” A distinctive aspect of the method is that it leverages the bootstrap — a tool ordinarily used for statistical inference — in order to serve a computational purpose. To our knowledge, it is also the first error estimation method for randomized LS that is supported theoretical guarantees. Furthermore, the method does not add much cost to an underlying sketching algorithm, and it has been shown to perform well on several examples.

Lopes is partially supported by NSF grant DMS-1613218.

Outline of Appendices

The proof of Theorem 1 is decomposed into two parts, with the bounds for IHS and CS being handled in Appendices A and B respectively.

A Proof of Theorem 1 for Iterative Hessian Sketch

To make the structure of the proof clearer, the main ingredients are combined in Appendix A.1. The lower-level arguments are given in Appendix A.2.

Since we will often need to condition on the sketching matrices $S_{1},\dots,S_{t}$ in the IHS algorithm, we define $\mathcal{S}_{k}:=\{S_{1},\dots,S_{k}\}$ for any $1\leq k\leq t$, and put $\mathcal{S}_{0}:=\emptyset$.

A.1 High-level proof of the bound (23)

Next, let $\varepsilon_{1}^{*},\dots,\varepsilon_{B}^{*}$ be the samples generated by Algorithm 2, and define the empirical distribution function

In Proposition 2 below, we show that as $n\to\infty$,

Establishing this limit is the most difficult part of the proof. Next, for any number $p\in(0,1)$, and any distribution function $G$, define the quantile function $G^{-1}(p)=\inf\{\tau:G(\tau)\geq p\}$. Using this definition, as well as the limit (24), it follows that for any fixed $\delta\in(0,1-\alpha)$, the event

and in the last step we have used the basic fact $F_{n}(F_{n}^{-1}(p))\geq p$ for any $p\in(0,1)$. So, by taking the complement of the event $\{\|\hat{x}_{t}-x_{\textup{opt}}\|_{\circ}>\hat{\varepsilon}_{t}(\alpha)\}$, the previous bounds give

Taking the expectation of both sides with respect to $\mathcal{S}_{t-1}$ leads to

Since the left side above does not depend on the arbitrarily small number $\delta$, it follows that

and this implies the inequality (23). $\square$

Due to the fact that $\hat{\varepsilon}_{t}(\alpha)$ can be expressed as $\hat{F}_{n,B}^{-1}(1-\alpha)$, we have the basic inequality $\hat{F}_{n,B}(\hat{\varepsilon}_{t}(\alpha))\geq 1-\alpha$. Also, if we define $\hat{\delta}_{n}:=|\hat{F}_{n,B}(\hat{\varepsilon}_{t}(\alpha))-F_{n}(\hat{\varepsilon}_{t}(\alpha))|$, then

If the conditions of Theorem 1 hold, then the limit (24) holds.

(Note that this holds regardless of the rate at which $B$ diverges, and so no conditions on the relative sizes of $B$ and $n$ are needed.) So, due to the simple inequality

and this is the core aspect of the proof. This limit follows directly from Lemma 6, which can be found at the end of the next subsection. (Prior to Lemma 6, there are three other lemmas that assemble the main arguments.) $\square$

A.2 Lemmas supporting the proof of Proposition 2

In this section we will use some specialized notation. In addition, our proofs will rely on the convergence of conditional distributions, as reviewed below.

Define the normalized matrix $\bar{A}_{n}:=\frac{1}{\sqrt{n}}A$, as well as the following analogues of $H_{n}=\frac{1}{n}A^{\top}A$,

Lastly, when referring to the rows of $\sqrt{m}S_{t}$, we will omit the dependence on $t$ and write simply $s_{1},\dots,s_{m}$ for ease of notation.

Convergence of conditional distributions.

If a sequence of random vectors $V_{n}$ converges in distribution to a random vector $V$, we write $\mathcal{L}(V_{n})\xrightarrow{\ d\ }\mathcal{L}(V)$. In some situations, we will also need to discuss convergence of conditional distributions. To review the meaning of this notion, let $d_{\text{LP}}(\mathcal{L}(V_{n}),\mathcal{L}(V))$ denote the Lévy-Prohorov distance (Dudley 2002, p. 394) between the distributions $\mathcal{L}(V_{n})$ and $\mathcal{L}(V)$, and note the basic fact that $d_{\text{LP}}(\mathcal{L}(V_{n}),\mathcal{L}(V))\to 0$ if and only if $\mathcal{L}(V_{n})\xrightarrow{\ d\ }\mathcal{L}(V)$. Now, suppose $U_{n}$ is another sequence of random vectors, and let $d_{\text{LP}}(\mathcal{L}(V_{n}|U_{n}),\mathcal{L}(V|U_{n}))$ denote the $d_{\text{LP}}$ distance between $\mathcal{L}(V_{n}|U_{n})$ and $\mathcal{L}(V|U_{n})$, which are random probability distributions. Likewise, the sequence $\{d_{\text{LP}}(\mathcal{L}(V_{n}|U_{n}),\mathcal{L}(V|U_{n}))\}_{n=1}^{\infty}$ may be regarded as a sequence of scalar random variables, and if it happens that this sequence converges to 0 in probability, then we say ‘$\mathcal{L}(V_{n}|U_{n})\xrightarrow{\ d\ }\mathcal{L}(V|U_{n})\text{ in probability}$’.

The rest of this subsection consists of the four lemmas needed to prove Proposition 2.

As a preparatory step towards applying the central limit theorem, we now show that $\mathsf{var}(\xi_{1,n})$ converges to a positive limit. Because each vector $s_{i}$ is composed of i.i.d. random variables, we may use an exact formula for the variance of quadratic forms (Bai and Silverstein 2004, eqn. 1.15), which leads to

Now that we have shown $\mathsf{var}(\xi_{1,n})$ converges to a limit, we verify that this limit is positive. Since we assume $\kappa>1$, it is clear that $\kappa-3>{\epsilon}_{0}-2$ for some fixed ${\epsilon}_{0}\in(0,1)$. Also, since the second term in line (39) represents the sum of the squares of the diagonal entries of $\bar{A}_{n}C\bar{A}_{n}^{\top}$, the sum of the two terms must be at least ${\epsilon}_{0}\|\bar{A}_{n}C\bar{A}_{n}^{\top}\|_{F}^{2}$. Therefore, the limit of $\mathsf{var}(\xi_{1,n})$ is lower-bounded by ${\epsilon}_{0}\|H_{\infty}^{1/2}CH_{\infty}^{1/2}\|_{F}^{2}$, and because $H_{\infty}$ is positive definite, this lower bound is positive when $C\neq 0.$

Suppose the conditions of Theorem 1 hold, and let $V$ be the random vector in statement of Lemma 3. Then, as $n\to\infty$,

Remark.

where $\sigma^{2}(C)$ is as defined beneath line (42), and also

To verify the limit (45), note that because $\xi_{1,n}^{*},\dots,\xi_{m,n}^{*}$ can be viewed as samples with replacement from the set $\{\xi_{1,n},\dots,\xi_{m,n}\}$, it follows that

where $\mu_{4,n}$ is the fourth central moment of $\xi_{1,n}$, i.e.

Using a general bound for the moments of quadratic forms (Bai and Silverstein 2010, Lemma B.26), this quantity can be bounded as

where $M_{n}:=(\bar{A}_{n}C\bar{A}_{n}^{\top})^{2}$. Since both of the traces above can be expressed in terms of the matrix $C\bar{A}_{n}^{\top}\bar{A}_{n}$, which converges to $CH_{\infty}$, it follows that $\mu_{4,n}=\mathcal{O}(1)$. Also, it was shown in the proof of Lemma 3 that $\mathsf{var}(\xi_{1,n})$ has a positive limit. Altogether, this completes the work needed to prove the limit (45).

Consequently, the argument based on the bound (41) in the proof of Lemma 3 may be re-used to show that the right hand side above tends to 0. $\square$

Remarks on notation.

Proof Recall that we assume $\sqrt{m}(H_{n}-H_{\infty})\to 0$, where $H_{\infty}$ is positive definite. Since the map $\phi$ is differentiable, it follows from the delta method (van der Vaart 1998, Theorem 3.1) and Lemma 3 that

where we define $W:=\phi^{\prime}_{0}(V),$ with $V$ being the random vector in Lemma 3, and $\phi^{\prime}_{0}$ denoting the differential of $\phi$ at the point $\textup{uvec}(H_{\infty})$.

Suppose the conditions of Theorem 1 hold, and let $Z=\textup{sym}(W)$, with $W$ being random vector in the statement of Lemma 5. Then, for almost every sequence of sets $\mathcal{S}_{t-1}$, the following limit holds as $n\to\infty$,

it straightforward to check that the following events are also equal

is upper bounded by the following supremum over $\mathcal{C}\in\mathscr{C}_{\text{convex}}$,

To conclude the proof of (53), it suffices to show that the previous expression converges to 0 as $n\to\infty$. For this purpose, we apply the general fact that if a sequence of random vectors $\zeta_{n}$ converges in distribution to a random vector $\zeta$, and if $\zeta$ has a multivariate normal distribution with a positive definite covariance matrix, then

B Proof of Theorem 1 for Classic Sketch

Also, letting $\varepsilon_{1}^{*},\dots,\varepsilon_{B}^{*}$ denote the samples generated by Algorithm 1, define

By using these functions in place of their IHS counterparts $F_{n}(\tau)$, $\hat{F}_{n}(\tau)$, and $\hat{F}_{n,B}(\tau)$, the argument at the beginning of Appendix A.1 can be essentially repeated to reach the conclusion

which implies the desired inequality (22). The only part of the argument that needs to be updated is to prove the analogue of Proposition 2 for the case of CS. In other words, it suffices to show that

Proving this limit will be handled with Proposition 8 below. $\square$

The proof of Proposition 8 relies on the following preliminary result. To introduce some notation, we will use the normalized gradient vector $\texttt{g}_{n}:=\textstyle\frac{1}{n}A^{\top}b$, and the analogues

Proof The proof of Lemmas 3 and 4 can be adapted to show that the following joint limits hold

as well as the following relations, which are straightforward to verify

These relations can be written in terms of the map $\Phi$ as

Next, recall the assumptions $\sqrt{m}(H_{n}-H_{\infty})\to 0$ and $\sqrt{m}(\texttt{g}_{n}-\texttt{g}_{\infty})\to 0$, and note that the map $\Phi$ is differentiable. Consequently, it follows from the delta method (van der Vaart 1998, Theorem 3.1), as well as the limit (67) that

where $\Phi_{0}^{\prime}$ denotes the differential of $\Phi$ evaluated at the point $(\textup{uvec}(H_{\infty}),\texttt{g}_{\infty})$. Furthermore, since $\Phi$ and $\Phi^{-1}$ are differentiable, it follows that $\Phi_{0}^{\prime}$ is an invertible linear map, which implies that the random vector $\Phi_{0}^{\prime}(V,U)$ has a non-zero covariance matrix (since $(V,U)$ does). Similarly, the delta method can be applied to the limit (68) to obtain

Finally, letting $Y=\Phi_{0}^{\prime}(V,U)$ completes the proof. $\square$

If the conditions of Theorem 1 hold, then the limit (63) holds

Since the random vector $Y$ has a multivariate normal distribution with a non-zero covariance matrix, it is straightforward to show that the random variable $\|Y\|_{\circ}$ has a continuous distribution function. So, it follows from Polya’s theorem (Bickel and Doksum 2007, Theorem B.7.7) that

which implies the limit (63) by the triangle inequality. $\square$

References