Unbiased estimators for random design regression

Michał Dereziński, Manfred K. Warmuth, Daniel Hsu

Introduction

In particular, the properties we want of the estimator $\widehat{\mathbf{w}}$ are the following.

Throughout the introduction we give some intuition about our results by discussing the one dimensional case. For example, consider the following $2\times 1$ fixed design problem:

The bias in least squares estimators is present even when each input component is drawn independently from a standard Gaussian. As an example, we let $d=5$ and set:

For this Gaussian setup we evaluate the bias of the least squares estimator produced for this problem by i.i.d. sampling of $k$ points. We do this by performing model averaging, i.e., producing many such estimators $\widehat{\mathbf{w}}_{1},\dots,\widehat{\mathbf{w}}_{T}$ independently, and looking at the estimation error of the average of those estimators $\widetilde{\mathbf{w}}\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{T}\sum_{t=1}^{T}\widehat{\mathbf{w}}_{t}$ :

Figure 1.1 (red curves) shows the experiment for several values of $k$ and a range of values of $T$ (each presented data point is an average over 50 runs). The i.i.d. sampled estimator is biased for any sample size (although the bias decreases with $k$ ), and therefore the averaged estimator clearly does not converge to the optimum. We next discuss how to construct an unbiased estimator (dashed blue curves), for which the estimation error of the averaged estimator exhibits $\frac{1}{T}$ convergence to zero (regardless of $k$ ). This type of convergence appears as a straight line on the log-log plot on Figure 1.1.

Even though i.i.d. sampling typically results in a biased least squares estimator, adding a volume-rescaled sample of size $d$ to the i.i.d. sample eliminates that bias altogether:

Perhaps surprisingly, volume-rescaled sampling may not lead to estimators with near-optimal loss guarantees: We show that for any $k\geq d$ there are distributions $D$ for which volume-rescaled sampling of size $k$ results in the linear least squares estimator having loss at least twice as large as the optimum loss (with probability at least 0.25). However, we remedy this bad behavior by composing a volume-rescaled sample of size $d$ with an i.i.d. leverage score sample of size $k-d$ . This composition achieves the following feat: It does not affect the unbiasedness of the estimator and, and it leads to good approximation properties. Specifically, in Theorem 3.1 we show that $k=O(d\log d+d/\epsilon)$ points are sufficient to construct an estimator $\widehat{\mathbf{w}}$ such that:

Note that an analogous loss bound is achievable for vanilla i.i.d. leverage score sampling, but (1) the estimators produced from leverage score sampling are biased, and (2) the expected loss bound holds only if we condition on a certain high-probability event (both of those are significant issues, e.g., in the context of model averaging). To show the expected loss bound that holds without conditioning and for an unbiased estimator, we break the analysis into two cases, depending on whether the high-probability event occurs. When it does not, then our analysis crucially relies on the expectation formulas we develop for volume-rescaled sampling. Note that the only expected loss bound previously developed for a volume-based sampling distribution was limited to fixed design, and required $d^{2}/\epsilon$ points to obtain an approximation factor of $1+\epsilon$ (Dereziński and Warmuth, 2018). To our knowledge, that analysis does not easily extend to $k>d$ , which is why our techniques are radically different.

Our work also leads to sampling algorithms which significantly improve on the state-of-the-art time complexity of volume-rescaled sampling, both in the fixed and random design settings, with further algorithmic implications for the broader class of determinantal point processes (see Section 1.2.3). We achieve this by introducing a new technique called distortion-free intermediate sampling: We first sample a larger pool of points based on approximate i.i.d. leverage scores and then down-sample from that pool to construct the volume-rescaled sample. We use rejection sampling for the down-sampling step to ensure exactness of the resulting overall sampling distribution. Surprisingly, this does not adversely affect the complexity because of the provably high acceptance rate during rejection sampling (see Theorem 5.6).

2 Applications of our results

While studying unbiased estimators for least squares regression is an old and classical problem, our new results have significant implications for modern data science, both from a computational and statistical perspective. We outline these implications below, along with some of the recent related work.

Model averaging is a standard technique for boosting the accuracy of a subsampled estimator by constructing multiple independent copies and then averaging them. This is particularly effective in parallel and distributed environments, where the computational cost of constructing multiple estimators is the same as the cost of computing one estimator. While model averaging has been proposed as a strategy for least squares regression (e.g., see Wang et al., 2017a), the bias which arises for commonly used estimators (e.g., based on i.i.d. sampling) constitutes a significant bottleneck for this approach.

2.2 Experimental design

In a recent follow-up work, Dereziński et al. (2019) used these ideas to develop a general framework for experimental design, which bridges the gap between the statistical perspective (linear response model) and the setting studied here (arbitrary responses), relying on our volume-rescaled sampling tool kit (in particular, Theorem 2.4). Furthermore, our strategy of combining volume-based sampling methods with i.i.d. importance sampling (e.g., leverage scores) has proven instrumental in developing randomized rounding methods for efficiently solving a range of experimental design problems (including A/C/D/V-optimal design, and Bayesian experimental design), drastically reducing their computational cost and improving the approximation quality, both for discrete (Nikolov et al., 2019; Dereziński et al., 2020b) and continuous domains (Poinas and Bardenet, 2020).

2.3 Determinantal point processes

3 Related work

A discrete variant of volume-rescaled sampling of size $k=d$ was introduced to computer science literature by Deshpande et al. (2006) for sampling from a finite set of $n$ vectors, with algorithms given later by Deshpande and Rademacher (2010); Guruswami and Sinop (2012). A first extension to samples of size $k>d$ is due to Avron and Boutsidis (2013), with algorithms by Li et al. (2017); Dereziński and Warmuth (2018); Dereziński et al. (2018), and additional applications in experimental design explored by Wang et al. (2017b); Nikolov et al. (2019); Mariet and Sra (2017). Prior to this work, the best known time complexity for this sampling method, called here discrete volume sampling, was $O(nd^{2})$ , as shown by Dereziński and Warmuth (2018). Here, we give an $O(nd\log n+d^{4}\log d)$ time algorithm.

As discussed in Section 1.2.3, volume-rescaled sampling of size $d$ is also known in the literature as a type of determinantal point process, called Projection DPP (to learn more, see Dereziński and Mahoney, 2021). Projection DPPs arise in many computational tasks outside of linear regression, such as dimensionality reduction (Belhadji et al., 2020), numerical integration (Bardenet and Hardy, 2020) and graph algorithms (Guenoche, 1983), therefore, efficient sampling algorithms for these distributions are of significant interest (Gautier et al., 2017). More broadly, determinantal point processes have found machine learning applications in recommendation systems (e.g., Gartrell et al., 2016), data summarization (e.g., Gong et al., 2014), stochastic optimization (e.g., Zhang et al., 2017; Mutný et al., 2020), and many others (see Kulesza and Taskar, 2012). The algorithmic technique of distortion-free intermediate sampling, introduced in this work, has already been applied beyond Projection DPPs (Dereziński et al., 2019; Calandriello et al., 2020), which makes it relevant to all of these applications.

The unbiasedness of least squares estimators under volume-based distributions was first explored in the context of sampling from finite datasets by Dereziński and Warmuth (2018), drawing on observations of Ben-Tal and Teboulle (1990). Focusing on small sample sizes, Dereziński and Warmuth (2018) proved multiplicative bounds for the expected loss under sample size $k=d$ with discrete volume sampling. Because the produced estimators are unbiased, averaging $d/\epsilon$ such estimators results in an unbiased estimator based on a sample of size $k=d^{2}/\epsilon$ with expected loss at most $1+\epsilon$ times the optimum at a total sampling cost of $O(nd^{2}\cdot d/\epsilon)$ . In contrast, our new techniques achieve an unbiased estimator with sample size $O(d\log d+d/\epsilon)$ and time complexity $O(nd\log n+d^{4}\log d+d^{3}/\epsilon)$ . Dereziński and Warmuth (2018) also showed additional variance bounds for discrete volume sampling under the assumption that the responses are linear functions of the input points plus white noise. We extend them here to arbitrary volume-rescaled sampling w.r.t. a distribution.

Our work greatly expands and generalizes the results of two conference papers: Dereziński et al. (2018, 2019). The first paper introduced the leverage score rescaling method in the limited context of discrete volume sampling, developed the new intermediate sampling algorithm, and proved the $O(d\log d+d/\epsilon)$ sample size bound for obtaining an unbiased estimator with a $(1+\epsilon)$ loss bound. Note that the original loss bound was shown to hold with a constant probability, as opposed to in expectation, which is a significant obstacle to using it in the context of model averaging. The second paper showed how to correct the bias of i.i.d. sampling using a small size $d$ volume-rescaled sample and refined the analysis of intermediate sampling. The current paper strengthens the loss bound of the first conference paper to the desired in-expectation form (this requires new technical tools such as Lemma 3.4), and generalizes it to the case of an arbitrary data distribution $D$ (Theorem 3.1). In the process, we develop new formulas for the expectation of the inverses and pseudoinverses of random matrices under volume-rescaled sampling (Theorems 2.8 and 2.9) and characterize the marginals of this distribution (Theorem 2.7). We also extend the decomposition property of volume-rescaled sampling given in the second conference paper (Theorem 2.4), thereby greatly simplifying our proofs. Finally, we give a new lower bound that complements our main results (Theorem 4.1).

Outline

In Section 6 we compare the performance of the algorithms discussed in this paper on some real datasets. We conclude with an overview and some open problems in Section 7.

Volume-rescaled sampling

In this section, we formally define volume-rescaled sampling and describe its basic properties. We then use it to introduce the central concept of this paper: an unbiased estimator for random design least squares regression.

For $k=d$ , this volume-rescaled sampling is a type of Determinantal Point Process known as Projection DPP (see Section 1.2.3; to learn more, see Dereziński and Mahoney, 2021). The case of $k>d$ can be viewed as an extension of that family of distributions.

Proof First, suppose that $k=d$ , in which case $\det(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{B})=\det(\mathbf{A})\det(\mathbf{B})$ . Recall that by definition the determinant can be written as:

which proves (2.1) for $k=d$ . The case of $k>d$ follows by induction via a standard determinantal formula:

where $(*)$ follows from the Cauchy-Binet formula. Finally, (2.2) can be derived from (2.1):

Proof For a full rank $d\times d$ matrix $\mathbf{A}$ we have $\mathbf{A}^{-1}=\mathbf{A}^{\dagger}$ and $\operatorname{\textnormal{adj}}(\mathbf{A})=\det(\mathbf{A})\mathbf{A}^{-1}$ . When $\mathbf{A}$ is not full rank but psd, then $\det(\mathbf{A})\mathbf{A}^{\dagger}=\mathbf{0}\preceq\operatorname{\textnormal{adj}}(\mathbf{A})$ . Thus Lemma 2.3 implies that

where $(*)$ becomes an equality if $\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}$ is full rank with probability 1.

2 Unbiased estimator for random design regression

where $\mathbf{X}\!\overset{i}{\leftarrow}\!\mathbf{y}$ is matrix $\mathbf{X}$ with column $i$ replaced by $\mathbf{y}$ . It follows that:

where we applied Lemma 2.3 to the pair of $d\times d$ matrices $\mathbf{A}=\mathbf{X}$ and $\mathbf{B}=\mathbf{X}\overset{i}{\leftarrow}\mathbf{y}$ . The case of $k>d$ follows by induction based on the following lemma shown by Dereziński and Warmuth (2018):

Proof of Theorem 2.8 The columns of $\bar{\mathbf{X}}^{\dagger}$ , equal $(\bar{\mathbf{X}}^{\scriptscriptstyle{\top}}\bar{\mathbf{X}})^{-1}\bar{\mathbf{x}}_{i}$ , are exchangeable, so

Loss bound for an unbiased estimator

Let $\bar{\mathbf{X}}$ and $\bar{\mathbf{y}}$ be distributed as in the theorem. Note that we can write the estimator $\widehat{\mathbf{w}}$ as follows:

Substituting $\widehat{\mathbf{w}}=\bar{\mathbf{X}}^{\dagger}\bar{\mathbf{y}}=(\bar{\mathbf{X}}^{\scriptscriptstyle{\top}}\bar{\mathbf{X}})^{-1}\bar{\mathbf{X}}^{\scriptscriptstyle{\top}}\bar{\mathbf{y}}$ for $\mathbf{w}$ , we additionally obtain:

We start by bounding the first term in (3.4), using a standard error decomposition (see Lemma 1 of Drineas et al., 2011):

where we used that $\|(\bar{\mathbf{X}}^{\scriptscriptstyle{\top}}\bar{\mathbf{X}})^{-1}\|\leq\|(\bar{\mathbf{X}}_{[s]^{c}}^{\scriptscriptstyle{\top}}\bar{\mathbf{X}}_{[s]^{c}})^{-1}\|\leq 4/k$ , when conditioned on $\mathcal{E}$ .

We apply Lemma 3.3 to the set $T=\{1,2\}$ and compute the determinant of a $2\times 2$ matrix:

Let us again use the notation of $\bar{\mathbf{r}}=\bar{\mathbf{y}}-\bar{\mathbf{X}}\mathbf{w}^{*}$ . To bound the second term in (3.4), we use a somewhat different decomposition of $\|\widehat{\mathbf{w}}-\mathbf{w}^{*}\|$ than we did in Part 1:

So, taking expectation, and noting that $\bar{\mathbf{X}}_{[s]}$ and $\bar{\mathbf{r}}_{[s]}$ are independent of $\mathcal{E}$ , we have:

To apply Theorem 2.9 again, we must disentangle the trace from $r_{s}^{2}$ , which is addressed in the following lemma proven at the end of the section.

Note that $\mathcal{E}^{\prime}$ implies $\mathcal{E}$ , and we can easily use Lemma 3.2 to bound its failure probability. Also, observe that, since the marginal distribution of each vector $\bar{\mathbf{x}}_{i}$ for $i\in[s]^{c}$ is the same, and the event $\mathcal{E}$ is invariant under permutation of the indices of these vectors, the marginal distributions of $\bar{r}_{i}^{2}=(\bar{y}_{i}-\bar{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}\mathbf{w}^{*})^{2}$ conditioned on $\neg\mathcal{E}$ are the same for each $i\in[s]^{c}$ , so:

where we used the fact that $\mathcal{E}^{\prime}$ is independent of $\bar{r}_{k}$ . Putting everything together, we conclude that:

The above result can also be achieved if we replace the exact leverage score sampling distribution with its approximation. As discussed in Section 5, producing samples from such approximation can be more practical in settings where exact leverage scores are too expensive to compute.

Proof of Lemma 3.4 Since the rows of $\bar{\mathbf{X}}$ are exchangeable, without loss of generality assume that $i=1$ . By definition of volume-rescaled sampling, we have:

where $\mathbf{x}^{-j}$ denotes vector $\mathbf{x}$ without the $j$ th entry and $\mathbf{X}^{-j}$ denotes matrix $\mathbf{X}$ without the $j$ th column, $(a)$ follows because $\det((\mathbf{X}_{-1}^{-j})^{\scriptscriptstyle{\top}}\mathbf{X}_{-1}^{-j})=0$ and $(b)$ comes from Lemma 2.3. Finally, to compute the trace, we sum up over $j$ :

Lower bounds

Here, we use the convention $0/0=0$ to handle the possibility of $S_{j}=\emptyset$ .

(Note that this is consistent with the case where $S_{j}=\emptyset$ .)

Let $A_{\bar{\mathbf{X}}}$ denote the event that there exists $j\in[d]$ such that no vector $\bar{\mathbf{x}}_{i}$ is equal to $\mathbf{e}_{j}$ . If $A_{\bar{\mathbf{X}}}$ holds then the $j$ th component of $\bar{\mathbf{X}}^{\dagger}\bar{\mathbf{y}}$ is so, setting $\gamma^{2}=\frac{\delta}{2d(1-\delta)}$ ,

Algorithms

For this theorem, given a p.d. matrix $\mathbf{A}$ , we use $\mathbf{A}^{\frac{1}{2}}$ to denote the unique lower triangular matrix with positive diagonal entries s.t. $\mathbf{A}^{\frac{1}{2}}(\mathbf{A}^{\frac{1}{2}})^{\scriptscriptstyle{\top}}=\mathbf{A}$ .

is matrix variate beta distributed, written as $\mathbf{U}\sim B_{d}(k_{1},k_{2})$ . The following was shown by Mitra (1970):

Now, Theorem 6.3.14 from Gupta and Nagar (1999) states that matrices $\mathbf{B}_{i}$ defined recursively as above can also be written as

In particular, we can construct them as $\mathbf{B}_{i}=\bar{\mathbf{x}}_{i}\bar{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}$ , where

2 Volume-rescaled sampling for arbitrary distributions

In this section, we present a general algorithm for volume-rescaled sampling which uses approximate leverage score sampling to generate a larger pool of points from which the smaller volume-rescaled sample can be drawn. The strategy introduced here, called distortion-free intermediate sampling, has since proven effective for sampling from other determinantal sampling distributions (Dereziński, 2019; Dereziński et al., 2019; Calandriello et al., 2020).

Next, we use the geometric-arithmetic mean inequality for the eigenvalues of matrix $\frac{1}{t}\widetilde{\mathbf{X}}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{X}}\widehat{\mathbf{\Sigma}}^{-1}$ to show that the Bernoulli sampling probability is bounded by 1:

Let $\widetilde{\mathbf{x}}^{\scriptscriptstyle{\top}}\sim D_{\cal\widetilde{\!X}}$ be distributed as a row vector of $\widetilde{\mathbf{X}}$ as sampled in line 3. The distribution of matrix $\widetilde{\mathbf{X}}$ returned by rejection sampling after exiting the repeat loop changes to:

So, using Lemma 2.3 on the matrix $\widetilde{\mathbf{X}}$ we obtain that:

since $\epsilon=\frac{1}{2\sqrt{d}}$ , where $(a)$ is the geometric-arithmetic mean inequality and $(b)$ is the Kantorovich inequality (Kantorovich, 1948) with $a=1-\epsilon$ and $b=1+\epsilon$ :

Now setting $t=2d^{2}$ we obtain the following lower bound for the acceptance probability:

3 Distributions with bounded support

We can lower bound the acceptance probability as follows:

4 Sampling from finite datasets

Experiments

The results confirm that our proposed leveraged volume sampling is as good or better than either of the baselines for any sample size $k$ . We can see that, in some of the examples, standard volume sampling exhibits bad behavior for larger sample sizes, as suggested by the lower bound of Theorem 4.2 (especially noticeable on bodyfat and cpusmall datasets). On the other hand, leverage score sampling exhibits poor performance for small sample sizes due to the coupon collector problem, which is most noticeable for abalone dataset, where we can see a very sharp transition after which leverage score sampling becomes effective. Neither of the variants of volume sampling suffers from this issue.

Conclusions

We showed that for any input distribution and $\epsilon>0$ , there is a random design consisting of $O(d\log d+d/\epsilon)$ points from which an unbiased estimator can be constructed whose expected square loss over the entire distribution is bounded by $1+\epsilon$ times the loss of the optimum. However, two main open problems remain. First, can the sample size bound be reduced to $O(d/\epsilon)$ ? This has already been done with a biased estimator by Chen and Price (2019), but finding an unbiased estimator of the smaller size remains open.

Second, the least squares estimator combined with i.i.d. leverage score sampling already achieves loss $1+\epsilon$ times the optimum with $O(d\log d+d/\epsilon)$ points. The resulting estimator is biased. However, in our preliminary experiments the bias of exact leverage score sampling is small and decreases rather quickly (unlike for uniform sampling, or even approximate leverage score sampling, where the bias can be significant). Thus, one of the key open problems is to quantify the bias of this method.

Michał Dereziński and Manfred K. Warmuth were supported by the NSF grant IIS-1619271. Michał Dereziński would also like to thank the NSF for funding via the NSF TRIPODS program. Daniel Hsu was supported by the NSF grant CCF-1740833 and a Sloan Research Fellowship. Part of this work was done while Manfred K. Warmuth was visiting Google Inc. in Zürich and Michał Dereziński was visiting the Simons Institute for the Theory of Computing. We would also like to acknowledge Eric Price for valuable discussions regarding this paper.

B Loss bound with approximate leverage scores

We now decompose the expectation into two terms depending on whether the event $\mathcal{E}$ occurs or not:

and the proof is divided into two parts, for handling the two terms.

We use the upper bound from (B.1). Event $\mathcal{E}$ implies that $\|(\bar{\mathbf{X}}^{\scriptscriptstyle{\top}}\bar{\mathbf{X}})^{-1}\|^{2}\leq 4^{2}/k^{2}$ . The second term in (B.1) is decomposed similarly as in (3.5), however bounding each of the obtained components will require a bit more care. Denoting $\bar{\mathbf{r}}=\bar{\mathbf{y}}-\bar{\mathbf{X}}\mathbf{w}^{*}$ , we have

This part follows identically as in the proof of Theorem 3.1, except that when applying Lemma 3.4, we use the fact that $\|\mathbf{x}\|^{2}\leq 3d$ , obtaining:

With the remaining steps same as in Theorem 3.1, this concludes the proof.