Random Projections For Large-Scale Regression

Introduction

Using a random projection, the data can be “compressed” either row- or column-wise. Row-wise compression was proposed and discussed in Zhou et al., 2007; Dhillon et al., 2013a; McWilliams et al., 2014b. These approaches replace the least-squares estimator

Column-wise compression addresses this later issue by reducing the problem to a $d$-dimensional optimization with $d\ll p$ by replacing the least-squares estimator

Random projections have also been considered under the aspect of preserving privacy (Blocki et al., 2012). By pre-multiplication with a random projection matrix as in (1) no observation in the resulting matrix can be identified with one of the original data points. Similarly, post-multiplication as in (2) produces new variables that do not reveal the realized values of the original variables.

In many applications the random projection used in practice falls under the class of Fast Johnson-Lindenstrauss Transforms (FJLT) (Ailon and Chazelle, 2006). One instance of such a fast projection is the Subsampled Randomized Hadamard Transform (SRHT) (Tropp, 2010). Due to its recursive definition, the matrix-vector product has a complexity of $\mathcal{O}(p\log(p))$, reducing the cost of the projection to $\mathcal{O}(np\log(p))$. Other proposals that lead to speedups compared to a Gaussian random projection matrix include random sign or sparse random projection matrices (Achlioptas, 2003). Notably, if the data matrix is sparse, using a sparse random projection can exploit sparse matrix operations. Depending on the number of non-zero elements in $\vec{X}$, one might prefer using a sparse random projection over a FJLT that cannot exploit sparsity in the data. Importantly, using $\vec{X}\vec{\phi}$ instead of $\vec{X}$ in our regression algorithm of choice can be disadvantageous if $\vec{X}$ is extremely sparse and $d$ cannot be chosen to be much smaller than $p$. (The projection dimension $d$ can be chosen by cross validation.) As the multiplication by $\vec{\phi}$ “densifies” the design matrix used in the learning algorithm the potential computational benefit of sparse data is not preserved.

For OLS and row-wise compression as in (1), where $n$ is very large and $p<m<n$, the SRHT (and similar FJLTs) can be understood as a subsampling algorithm. It preconditions the design matrix by rotating the observations to a basis where all points have approximately uniform leverage (Dhillon et al., 2013a). This justifies uniform subsampling in the projected space which is applied subsequent to the rotation in order to reduce the computational costs of the OLS estimation. Related ideas can be found in the way columns and rows of $\vec{X}$ are sampled in a CUR-matrix decomposition (Mahoney and Drineas, 2009). While the approach in Dhillon et al., 2013a focuses on the concept of leverage, McWilliams et al., 2014b propose an alternative scheme that allows for outliers in the data and makes use of the concept of influence (Cook, 1977). Here, random projections are used to approximate the influence of each observation which is then used in the subsampling scheme to determine which observations to include in the subsample.

Lastly, random projections have been used as an auxiliary tool. As an example, the goal of McWilliams et al., 2014a is to distribute ridge regression across variables with an algorithm called Loco. The design matrix is split across variables and the variables are distributed over processing units (workers). Random projections are used to preserve the dependencies between all variables in that each worker uses a randomly projected version of the variables residing on the other workers in addition to the set of variables assigned to itself. It then solves a ridge regression using this local design matrix. The solution is the concatenation of the coefficients found from each worker and the solution vector lies in the original space so that the coefficients are interpretable. Empirically, this scheme achieves large speedups while retaining good predictive accuracy. Using some of the ideas and results outlined in the current manuscript, one can show that the difference between the full solution and the coefficients returned by Loco is bounded.

Clearly, row- and column-wise compression can also be applied simultaneously or column-wise compression can be used together with subsampling of the data instead of row-wise compression. In the remaining sections, we will focus on the column-wise compression as it poses more difficult challenges in terms of statistical performance guarantees. While row-wise compression just reduces the effective sample size and can be expected to work in general settings as long as the compressed dimension $m<n$ is not too small (Zhou et al., 2007), column-wise compression can only work well if certain conditions on the data are satisfied and we will give an overview of these results. If not mentioned otherwise, we will refer with compressed regression and random projections to the column-wise compression.

The structure of the manuscript is as follows: We will give an overview of bounds on the estimation accuracy in the following section 2, including both known results and new contributions in the form of tighter bounds. In Section 3 we will discuss the possibility and properties of variance-reducing averaging schemes, where estimators based on different realized random projections are aggregated. Finally, Section 4 concludes the manuscript with a short discussion.

Theoretical Results

We will discuss in the following the properties of the column-wise compressed estimator as in (2), which is defined as

where we assume that $\phi$ has i.i.d. $\mathcal{N}(0,1/d)$ entries. This estimator will be referred to as the compressed least squares estimator (CLSE) in the following. We will focus on the unpenalized form as in (3) but note that similar results also apply to estimators that put an additional penalty on the coefficients $\beta$ or $\gamma$. Due to the isotropy of the random projection, a ridge-type penalty as in Lu et al., 2013 and McWilliams et al., 2014a is perhaps a natural choice. An interesting summary of the bounds on random projections is on the other hand that the random projection as in (3) already acts as a regularization and the theoretical properties of (3) are very much related to the properties of a ridge-type estimator of the coefficient vector in the absence of random projections.

We will restrict discussion of the properties mostly to the mean squared error (MSE)

The following result of (Kabán, 2014) gives a bound on the MSE for fixed design.

Compared with (Maillard and Munos, 2009), the result removes an unnecessary $\mathcal{O}(\log(n))$ term and demonstrates the $\mathcal{O}(1/d)$ behaviour of the bias. The result also illustrates the tradeoffs when choosing a suitable dimension $d$ for the projection. Increasing $d$ will lead to a $1/d$ reduction in the bias terms but lead to a linear increase in the estimation error (which is proportional to the dimension in which the least-squares estimation is performed). An optimal bound can only be achieved with a value of $d$ hat depends on the unknown signal and in practice one would typically use cross-validation to make the choice of the dimension of the projection.

One issue with the bound in Theorem 2.1 is that the bound on the bias term in the noiseless case ($Y=\vec{X}\beta$)

is usually weaker than the trivial bound (by setting $\hat{\beta}_{d}^{\vec{\phi}}=0$) of

for most values of $d<p$. By improving the bound, it is also possible to point out the similarities between ridge regression and compressed least squares.

The improvement in the bound rests on a small modification in the original proof in (Kabán, 2014). The idea is to bound the bias term of (4) by optimizing over the upper bound given in the foregoing theorem. Specifically, one can use the inequality

To simplify the exposition we will from now on always assume we have rotated the design matrix to an orthogonal design so that the Gram matrix is diagonal:

The $w_{i}$ are shrinkage factors. By defining the proportion of the total variance observed in the direction of the $i$-th principal component as

we can rewrite the shrinkage factors in the foregoing theorem as

Analyzing this term shows that the shrinkage is stronger in directions of high variance compared to directions of low variance. To explain this relation in a bit more detail we compare it to ridge regression. The MSE of ridge regression with penalty term $\lambda\|\beta\|_{2}^{2}$ is given by

Imagine that the signal lives on the space spanned by the first $q$ principal directions, that is $\beta_{i}=0$ for $i>q$. The best MSE we could then achieve is $\sigma^{2}q$ by running a regression on the first $q$ first principal directions. For random projections, we can see that we can indeed reduce the bias term to nearly zero by forcing $w_{i}\approx 0$ for $i=1,\ldots,q$. This requires $d\gg q$ as the bias factors will then vanish like $1/d$. Ridge regression on the other hand requires that the penalty $\lambda$ is smaller than the $q$-th largest eigenvalue $\lambda_{q}$ (to reduce the bias on the first $q$ directions) but large enough to render the variance factor $\lambda_{i}/(\lambda_{i}+\lambda)$ very small for $i>q$. The tradeoff in choosing the penalty $\lambda$ in ridge regression and choosing the dimension $d$ for random projections is thus very similar. The number of directions for which the eigenvalue $\lambda_{i}$ is larger than the penalty $\lambda$ in ridge corresponds to the effective dimension and will yield the same variance bound as in random projections. The analogy between the MSE bounds (9) for random projections and (13) for ridge regression illustrates thus a close relationship between compressed least squares and ridge regression or principal component regression, similar to Dhillon et al., 2013b.

Instead of an upper bound for the MSE of CLSE as in Maillard and Munos, 2009 and Kabán, 2014, we will in the following try to derive explicit expressions for the MSE, following the ideas in Kabán, 2014 and Marzetta et al., 2011 and we give a closed form MSE in the case of orthonormal predictors. The derivation will make use of the following notation:

The next Lemma (Marzetta et al., 2011) summarizes the main properties of $\vec{\phi}_{d}^{\vec{X}}$ and $T_{d}^{\vec{\phi}}$.

$\,\,\,(\vec{\phi}_{d}^{\vec{X}})^{\prime}=\vec{\phi}_{d}^{\vec{X}}$ (symmetric),

$\,\,\,\vec{\phi}_{d}^{\vec{X}}\vec{X}^{\prime}\vec{X}\vec{\phi}_{d}^{\vec{X}}=\vec{\phi}_{d}^{\vec{X}}$ (projection),

if $\Sigma=\vec{X}^{\prime}\vec{X}$ is diagonal $\Rightarrow$ $T_{d}^{\vec{\phi}}$ is diagonal.

The important point of this lemma is that when we assume orthogonal design then $T_{d}^{\vec{\phi}}$ is diagonal. We will denote this by

where the terms $\eta_{i}$ are well defined but without an explicit representation.

A quick calculation reveals the following theorem:

By comparing coefficients in Theorems 2 and 3, we obtain the following corollary.

As already mentioned in general we cannot give a closed-form expression for the terms $\eta_{i}$ in general. However, for some special cases (26) can help us to get to an exact form of the MSE of CLSE. If we assume orthonormal design ($\Sigma=CI_{p\times p}$) then we have that $\lambda_{i}/\eta_{i}$ is a constant for all $i$ and and thus, by (26), we have $\eta_{i}=Cp/d$. This gives

and thus we end up with a closed form MSE for this special case.

Providing the exact mean-squared errors allows us to quantify the conservativeness of the upper bounds. The upper bound has been shown to give a good approximation for small dimensions $d$ of the projection and for the signal contained in the larger eigenvalues.

Averaged Compressed Least Squares

We have so far looked only into compressed least squares estimator with one single random projection. An issue in practice of the compressed least squares estimator is its variance due to the random projection as an additional source of randomness. This variance can be reduced by averaging multiple compressed least squares estimates coming from different random projections. In this section we will show some properties of the averaged compressed least squares (ACLSE) estimator and discuss its advantage over the CLSE.

One major advantage of this estimator is that it can be calculated in parallel with the minimal number of two communications, one to send the data and one to receive the result. This means that the asymptotic computational cost of $\hat{\beta}_{d}^{K}$ is equal to the cost of $\hat{\beta}_{d}^{\vec{\phi}}$ if calculations are done on $K$ different processors. To investigate the MSE of $\hat{\beta}_{d}^{K}$, we restrict ourselves for simplicity to the limit case

and instead only investigate $\hat{\beta}_{d}$. The reasoning being that for large enough values of $K$ (say $K>100$) the behaviour of $\hat{\beta}_{d}$ is very similar to $\hat{\beta}_{d}^{K}$. The exact form of the MSE in terms of the $\eta_{i}$’s is given in Kabán, 2014. Here we build on these results and give an explicit upper bound for the MSE.

The MSE of $\hat{\beta}_{d}$ can be bounded from above by

where the $w_{i}$’s are given (as in Theorem 2.1) by

Comparing averaging to the case where we only have one single estimator we see that there are two differences: First the variance due to the model noise $\varepsilon$ turns into $\sigma^{2}\tau$ with $\tau\in[d^{2}/p,d]$, thus $\tau\leq d$. Secondly the shrinkage factors $w_{i}$ in the bias are now squared, which in total means that the MSE of $\hat{\beta}_{d}$ is always smaller or equal to the MSE of a single estimator $\hat{\beta}_{d}^{\vec{\phi}}$.

We investigate the behavior of $\tau$ as a function of $d$ in three different situations (Figure 2). We first look at two extreme cases of covariance matrices for which the respective upper and lower bounds $[d^{2}/p,d]$ for $\tau$ are achieved. For the lower bound, let $\Sigma=I_{p\times p}$ be orthonormal. Then $\lambda_{i}/\eta_{i}=c$ for all $i$, as above. From

we get $\lambda_{i}/\eta_{i}=d/p$. This leads to

We will not be able to reproduce the upper bound exactly for all $d\leq p$. But we can show that for any $d$ there exists a covariance matrix $\Sigma$, such that the upper bound is reached. The idea is to consider a covariance matrix that has equal variance in the first $d$ direction and almost zero in the remaining $p-d$. Define the diagonal covariance matrix

The same way we define $\beta_{d}$, $\beta_{r}$, $\vec{X}_{d}$ and $\vec{X}_{r}$. Now we bound the approximation error of $\hat{\beta}_{d}^{\Phi}$ to extract information about $\lambda_{i}/\eta_{i}$. Assume a squared data matrix ($n=p$) $\vec{X}=\sqrt{\Sigma}$, then

where $C$ is independent of $\epsilon$ and bounded since the expectation of the smallest and largest singular values of a random projection is bounded. This means that the approximation error decreases to zero as we let $\epsilon\rightarrow 0$. Applying this to the closed form for the MSE of $\hat{\beta}_{d}^{\Phi}$ we have that

has to go to zero as $\epsilon\rightarrow 0$, which in turn implies

and thus $\lim_{\epsilon\rightarrow 0}\;\lambda_{i}/\eta_{i}=1$ for all $i\in\{1,...,d\}$. This finally yields a limit

This illustrates that the lower bound $d^{2}/p$ and upper bound $d$ for the variance factor $\tau$ can both be attained. Simulations suggest that $\tau$ is usually close to the lower bound, where the variance of the estimator is reduced by a factor $d/p$ compared to a single iteration of a compressed least-squares estimator, which is on top of the reduction in the bias error term. This shows, perhaps unsurprisingly, that averaging over random projection estimators improves the mean-squared error in a Rao-Blackwellization sense. We have quantified the improvement. In practice, one would have to decide whether to run multiple versions of a compressed least-squares regression in parallel or run a single random projection with a perhaps larger embedding dimension. The computational effort and statistical error tradeoffs will depend on the implementation but the bounds above will give a good basis for a decision.

Discussion

We discussed some known results about the properties of compressed least-squares estimation and proposed possible tighter bounds and exact results for the mean-squared error. While the exact results do not have an explicit representation, they allow nevertheless to quantify the conservative nature of the upper bounds on the error. Moreover, the shown results allow to show a strong similarity of the error of compressed least squares, ridge and principal component regression. We also discussed the advantages of a form of Rao-Blackwellization, where multiple compressed least-square estimators are averaged over multiple random projections. The latter averaging procedure also allows to compute the estimator trivially in a distributed way and is thus often better suited for large-scale regression analysis. The averaging methodology also motivates the use of compressed least squares in the high dimensional setting where it performs similar to ridge regression and the use of multiple random projection will reduce the variance and result in a non random estimator in the limit, which presents a computationally attractive alternative to ridge regression.

References

Appendix

Finally a rather lengthy but straightforward calculation leads to

This last expression we can calculate following the same path as in Theorem 1:

where $\Sigma=X^{\prime}X$. Next we minimize the above expression w.r.t $v$. For this we take the derivative w.r.t. $v$ and then we zero the whole expression. This yields

Define the notation $s=\operatorname{trace}(\Sigma)$. We now plug this back into the original expression and get

by combining the summands we get for $w_{i}$ the expression mentioned in the theorem. ∎

The first term in the last line equals $\sum_{i=1}^{p}\beta_{i}^{2}\lambda_{i}^{2}/\eta_{i}$. The second can be calculated in two ways, both relying on the shuffling property of the trace operator:

Adding the first version to the expectation from above we get the exact expected mean squared error. Setting both versions equal we get the equation

First a simple calculation (Kabán, 2014) using the closed form solution gives the following equation:

Now using the corollary from the last section we can bound the second term the following way:

Now note that since $\lambda_{i}/\eta_{i}\leq 1$ we have

The problem is symmetric in each coordinate and thus $t_{i}=c$. Plugging this into the linear sum gives $c=d/p$ and we calculate the quadratic term to give the result claimed in the theorem. ∎