Improved matrix algorithms via the Subsampled Randomized Hadamard Transform

Introduction

Fix an integer $n=2^{p}$, for $p=1,2,3,...$. The (non-normalized) $n\times n$ matrix of the Walsh–Hadamard transform is defined recursively as,

Fix integers $r$ and $n=2^{p}$ with $r<n$ and $p=1,2,3,...$. An SRHT matrix is an $r\times n$ matrix of the form

The purpose of this article is to analyze the theoretical performance of an SRHT-based randomized low-rank approximation algorithm introduced in and analyzed in . Our analysis (see Theorem 4) provides sharper approximation bounds than those in .

Our study should also be viewed as followup to the work of Drineas et al. and on designing fast approximation algorithms for solving least-squares regression problems. One of the two algorithms presented in employs the SRHT to quickly reduce the dimension of the least squares problem and then solves the smaller problem with a direct least-squares solver, while use the SRHT to design a good preconditioner for an iterative method, which is then used to solve the regression problem. The results in this article along with the work in have implications in all these studies . We discuss these implications in Section 3.

Finally, notice that the SRHT is defined only when the matrix dimension is a power of two. An alternative option is to use other structured orthonormal randomized transforms such as the discrete cosine transform (DCT) or the discrete Hartley transform (DHT) , whose entries are on the order of $n^{-1/2}.$ All these transforms do not place any restrictions on the size of the matrix. The results of this paper - with minimal effort - can be extended unchanged to encompass these transforms. To see this, notice that Lemma 3.3 in remains unchanged for all these orthogonal transforms. Thus Lemma 6 in our work as well as all other results presented in this article are true for these orthogonal transforms as well.

2 Roadmap

This article is structured as follows. Section 1.3 introduces the notation. In Section 2, we present our main results on the quality of SRHT low-rank approximations and compare them to prior results in the literature. In Section 3, we discuss two approaches to least-squares regression involving SRHT dimensionality-reduction. Section 4 first recalls known facts on the application of SRHTs to orthogonal matrices and then presents new results on the application of SRHTs to general matrices and the approximation of matrix multiplication using SRHTs under the Frobenius norm. Section 5 contains the proofs of our two main theorems presented in Sections 2 and 3. We conclude the paper with an experimental evaluation of the SRHT low-rank approximation algorithm in Section 6.

3 Preliminaries

Low-rank matrix approximation using SRHTs

Similarly, the same setup ensures that with probability at least $1-5\delta,$ the following spectral norm bounds hold simultaneously:

The first two Frobenius norm bounds in this theorem (residual error analysis) are slightly stronger than the best bounds appearing in prior efforts . The spectral norm bounds on the residual error are significantly better than the bounds presented in prior work and shed light on an open question mentioned in . We do not, however, claim that the error bounds provided are the tightest possible. Certainly the specific constants ($22,6,$ etc.) in the error estimates are not optimized.

We now present a detailed comparison of the guarantees given in Theorem 4 with those available in the existing literature.

To put our result into perspective, we compare it to prior efforts at analyzing the SRHT algorithm introduced above. Halko et al. argue that if $r$ satisfies

then when $r$ is chosen according to Theorem 4 the quantity in Eqn. (4) is much smaller than

1.2 Nguyen et al. [37]

with probability of success at least $1-5\delta$, one requires

1.3 The subsampled randomized Fourier transform (SRFT)

then with probability at least $1-3/\beta$ ($\omega=\max\{m,n\}$),

1.4 Two alternative dimensionality-reduction algorithms

We use the notion of the stable rank of a matrix,

When $r>k+4,$ Theorem 10.7 and Corollary 10.9 in imply that, when using Gaussian sampling, with probability at least $1-2\cdot 32^{-(r-k)}-e^{{-(r-k+1)\over 2}}$,

and with probability at least $1-3e^{-(r-k)}$,

to ensure a $(6+\varepsilon)$ relative spectral error approximation. We note that the SRHT algorithm can be used to obtain relative spectral error approximations of matrices with arbitrary stable rank at the cost of increasing $r$ (the same is of course true for the Gaussian and random sign algorithms).

Least squares regression

We now show how one can use the SRHT to solve least squares problems of the form

During the last decade, researchers have developed several randomized algorithms that (approximately) solve the regression problem in less running time than the approaches mentioned above . We refer the reader to Section 3.3 in for a survey of these methods. The fastest non-iterative method is in while the fastest iterative algorithm is in . Both approaches employ the Subsampled Randomized Hadamard Transform.

shows that with probability at least $0.8$,

Below, we provide a novel analysis of this SRHT least squares algorithm which shows that one needs asymptotically fewer samples $r$. This immediately implies an improvement on the running time of the algorithm. Additionally, we show logarithmic dependence on the failure probability.

We prove this theorem in Section 5.3. Another possibility to obtain a better analysis of the method of Drineas et al. is to use Lemma 10 in this article, which was proved in and presents bounds for sampling without replacement. This analysis is not straightforward and is beyond the scope of this paper.

2 Iterative methods

The analysis of Blendenpik was recently improved in . More specifically, Corollary 3.11 in , along with Lemma 7 in our manuscript, which gives a bound on the coherence, show that if

then, with probability at least $1-2\delta$,

Then, with probability at least $1-2\delta$,

Finally, notice that we form the SRHT by uniform sampling without replacement while Blendenpik samples the columns of the randomized Hadamard matrix with replacement. A different sampling scheme - Bernoulli sampling - was analyzed in Theorem 6.1 in and Section 4 in .

then with probability at least $1-\delta$,

Matrix Computations with SRHT matrices

These perturbed orthogonal matrices have small norm precisely when their singular values are close to those of the original orthogonal matrices.

In this section, we collect known results on how the singular values of a matrix with orthonormal rows are affected by postmultiplication by an SRHT matrix.

It has recently been shown by Tropp that, if the SRHT matrix is of sufficiently large dimensions, post-multiplying a short-fat matrix with orthonormal rows with an SRHT matrix preserves the singular values of the orthonormal matrix, with high probability, up to a small multiplicative factor. The following lemma is essentially a restatement of Theorem 3.1 in , but we include a full proof (later in this subsection) for completeness.

Then, with probability at least $1-3\delta$, for all $i=1,...,k$,

To prove Lemma 6 we need one more result on uniform random sampling (without replacement) of rows from tall-thin matrices with orthonormal columns.

Replacing $\varepsilon$ with $\sqrt{\varepsilon}$ and using the bound on $r$ in Eqn. (7) concludes the proof.

1.2 SRHTs by uniform sampling with replacement

Lemma 6 and Lemma 8 analyze uniform random sampling without replacement. Below, we present the analogs of these two lemmas for uniform random sampling with replacement. Lemma 9 is essentially a restatement of Algorithm 2 (with the probabilities set to $1/m$) along with the third point in Remark 3.9 and Lemma 2.1 (with $\alpha=\sqrt{n/r}$) in .

Then, with probability at least $1-3\delta,$ for all $i=1,...,k$,

Replacing $\varepsilon$ with $\sqrt{\varepsilon}$ and using the bound on $r$ in Eqn. (8) concludes the proof.

2 SRHTs applied to general matrices

Our main tool is a generalization of Lemma 7 that states that the maximum column norm of a matrix to which an SRHT has been applied is, with high probability, not much larger than the root mean-squared average of the column norms of the original matrix.

Suppose $f$ is a convex function on vectors that satisfies the Lipschitz bound

Let $\bm{\varepsilon}$ be a Rademacher vector. For all $t\geq 0,$

As an interesting aside, we note that just as Lemma 6, which states that the SRHT essentially preserves the singular value of matrices with orthonormal rows and an aspect ratio of $k/n$, follows from Lemma 7, Lemma 11 implies that the SRHT essentially preserves the singular values of general rectangular matrices with the same aspect ratio. This can be shown using, e.g., the results on the effects of column sampling on the singular values of matrices from [24, Section 6].

2.2 SRHT preserves the spectral norm

The following lemma shows that even if the aspect ratio is larger than $k/n,$ the SRHT does not substantially increase the spectral norm of a matrix.

To establish Lemma 13, we use the following Chernoff bound for sampling matrices without replacement.

Let $\mathcal{X}$ be a finite set of positive-semidefinite matrices with dimension $k,$ and suppose that

When $E$ holds, for all $j=1,2,\ldots,n,$

Take the parameter $\nu$ in Lemma 14 to be

The second inequality holds because $\nu\geq 4$ implies that $(1+\nu)\ln(1+\nu)\geq\nu\cdot\tfrac{5}{4}\ln 5.$

2.3 SRHT preserves the Frobenius norm

Similarly, the SRHT is unlikely to substantially increase the Frobenius norm of a matrix.

Applying Lemma 14 conditioned on $E,$ we conclude that

2.4 SRHT preserves matrix multiplication

Finally, we prove a novel result on approximate matrix multiplication involving SRHT matrices.

Proof of Lemma 16

To prove the Lemma, we first develop a generic result for approximate matrix multiplication via uniform sampling (without replacement) of the columns and the rows of the two matrices involved in the product (see Lemma 18 below). Lemma 16 is a simple instance of this generic result. We mention that Lemma 3.2.8 in gives a similar result for approximate matrix multiplication, which, however gives a bound for the expected value of the error term, while our Lemma 16 gives a comparable bound which holds with high probability. To prove Lemma 18, we use the following vector Bernstein inequality for sampling without replacement in Banach spaces; this result follows directly from a similar inequality for sampling with replacement established by Gross in .

Let $\mathcal{V}$ be a collection of $n$ vectors in a normed space with norm $\mbox{}\left|\cdot\right|.$ Choose $\bm{V}_{1},\ldots,\bm{V}_{r}$ from $\mathcal{V}$ uniformly at random without replacement. Also choose $\bm{V}_{1}^{\prime},\ldots,\bm{V}_{r}^{\prime}$ from $\mathcal{V}$ uniformly at random with replacement. Let

We proceed by developing a bound on the moment generating function (mgf) of

This mgf is controlled by the mgf of a similar sum where the vectors are sampled with replacement. That is, for $\lambda\geq 0,$

This follows from a classical observation due to Hoeffding (see also for a more modern exposition) that for any convex $\R$-valued function $g,$

In the proof of Theorem 12 in , Gross establishes that any random variable $Z$ whose mgf is less than the righthand side of Eqn. (15) satisfies a tail inequality of the form

and $M$ almost surely bounds $\mbox{}\left|\bm{V}_{i}^{\prime}-{{\mathbf{E}}\left[\bm{V}_{1}^{\prime}\right]}\right|$ for all $i=1,\ldots,r.$ To apply this result, note that for all $i=1,\ldots,r,$

Also take $\bm{V}_{1}^{\prime\prime}$ to be an i.i.d. copy of $\bm{V}_{1}^{\prime}$ and observe that, by Jensen’s inequality,

The bound given in the statement of Lemma 17 follows from taking $s^{2}=\sigma^{2}$ and $M=B$ in Eqn. (16).

This vector Bernstein inequality gives us a tail bound on the Frobenius error of a simple approximate matrix multiplication scheme based upon column and row sampling.

have the same distribution, therefore any probabilistic bound developed for the latter holds for the former. The conclusion of the lemma follows from applying Lemma 17 to bound the second quantity.

We calculate the variance-like term in Lemma 17, $4r{{\mathbf{E}}\left[\mbox{}\|\bm{V}_{1}\|_{2}^{2}\right]}:$

In doing so, we will use the notation $\mathbf{E}_{A,B,\ldots}\left[{C}\right]$ to denote the conditional expectation of a random variable $C$ with respect to the random variables $A,B,\ldots.$ Recall that a Rademacher vector is a random vector whose entries are independent and take the values $\pm 1$ with equal probability. Let $\bm{\varepsilon}$ be a Rademacher vector of length $r$ and sample $\bm{V}_{1}^{\prime},\ldots,\bm{V}_{r}^{\prime}$ and $\bm{V}_{1}^{\prime\prime},\ldots,\bm{V}_{r}^{\prime\prime}$ uniformly at random from $\mathcal{V}$ with replacement. Now $\mu$ can be bounded as follows:

The first inequality is Jensen’s, and the following equality holds because the components of the sequence $\{\bm{V}_{i}^{\prime}-\bm{V}_{i}^{\prime\prime}\}$ are symmetric and independent. The next two manipulations are the triangle inequality and Jensen’s inequality. This stage of the estimate is concluded by conditioning and using the orthogonality of the Rademacher variables. Next, the triangle inequality and the fact that ${{\mathbf{E}}\left[\mbox{}\|\bm{V}_{1}^{\prime}\|_{2}^{2}\right]}={{\mathbf{E}}\left[\mbox{}\|\bm{V}_{1}\|_{2}^{2}\right]}$ allow us to further simplify the estimate of $\mu:$

The stipulated tail bound follows from applying Lemma 17 with our estimates for $B$, $\sigma^{2},$ and $\mu.$

Lemma 16 now follows from this result on matrix multiplication.

To apply Lemma 18, we first condition on the event that the SRHT equalizes the column norms of our matrices. Namely, we observe that, from Lemma 11, with probability at least $1-2\delta,$

Conditioning on these nice interactions, we choose the parameters $\sigma$ and $B$ in Lemma 18. We first take

so this choice of $\sigma$ satisfies the inequality stipulated in Lemma 18. Next we choose

For simplicity, let $\gamma=8\ln(n/\delta).$ With these choices for $\sigma^{2}$ and $B,$

Now, referring to Eqn. (18), identify the numerator as $\sqrt{r}\sigma$ to see that

Apply Lemma 18 to see that, when Eqns. (17) hold and $0\leq R\sigma\leq\sigma^{2}/B,$

From our lower bound on $\sigma^{2}/B,$ we know that the condition $R\sigma\leq\sigma^{2}/B$ is satisfied when

Also, we established above that Eqns. (17) hold with probability at least $1-2\delta.$ From these two facts, it follows that when $0\leq R\leq\sqrt{r}/(1+\sqrt{8\ln(n/\delta)}),$

The tail bound given in the statement of Lemma 16 follows from substituting our estimate of $\sigma.$

Proofs of our main Theorems

Lemma 20 is the analog of the Pythagoras theorem in the matrix setting. A proof of this lemma can be found in . Lemma 21 is an immediate corollary of Matrix-Pythogoras.

1.2 Low-rank matrix approximation based on projections

The low-rank matrix approximation algorithm investigated in this paper is an instance of a wider class of low-rank approximation schemes wherein a matrix is projected onto a subspace spanned by some linear combination of its columns. The problem of providing a general framework for studying the error of such projection schemes is well studied . The following result appeared as Lemma 7 in (see also Theorem 9.1 in ).

This lemma provides an upper bound for the residual error of the low-rank matrix approximation obtained via projections. We now prove a new result for the forward error.

1.3 Least squares regression based on projections

Similarly, one of the two SRHT least squares regression algorithms analyzed in this article is an instance of a wider class of approximation algorithms where the dimensions of the input matrix and the vector of the regression problem are reduced via pre-multiplication with a random matrix. Lemma 9 in provides a general framework for the analysis of such projection algorithms.

The following lemma is a restatement of Lemma 2, along with Eqn. (9) and Eqn. (11) in . It gives a bound on the forward error of the approximation of a least-squares problem that is obtained via projections. In the parameters $\alpha$ and $\beta$ are fixed to $\alpha=1/\sqrt{2}$ and $\beta=\varepsilon/2$, for some parameter $0<\varepsilon<1$. Showing the result for general $\alpha>0$ and $\beta>$ is straightforward, hence a detailed proof is omitted.

2 Proof of Theorem 4

We continue by bounding the second term in the right hand side of the above inequality,

Use the lower bound on $r$ to justify the estimate

The remaining estimates in the third inequality follow from applying Lemma 6 (keeping in mind our lower bound on $r$) to obtain

Combining (24) with the bound on $S$, we obtain

Taking the square-roots of both sides and using the fact that $\sqrt{1+22\varepsilon}\leq 1+22\varepsilon$ gives the bound

Eqn. (i) in the theorem follows directly from this:

Taking the square-roots of both sides and using the fact that $\sqrt{1+22\varepsilon}\leq 1+22\varepsilon$ gives Eqn. (iii).

where the first inequality follows by the triangle inequality and the second by using the bound obtained in Eqn. (ii) in the theorem.

The failure probability in the theorem follows from a union bound on all the probabilistic events involved in bounding $S$.

2.2 Spectral norm bounds

We now prove the spectral norm bounds in Theorem 4 (i.e Eqns. (v), (vi), (vii), and (viii)). Lemma 6 implies that, with this choice of $r,$

Also, the spectral norm bound in Lemma 19 implies

We now provide an upper bound for $\sqrt{Z}$ where $Z$ is the scalar

with probability at least $1-5\delta.$ Using that $\varepsilon<1/3$, we see that $(1-\sqrt{\varepsilon})^{-1}<3$, so

Use the subadditivity of the square-root function and rearrange the spectral and Frobenius norm terms to obtain that

Apply Eqn. (25) to arrive at Eqn. (v) in the theorem,

Take the square root of both sides of Eqn. (26), use the subadditivity of the square root function, and use the bound for $\sqrt{Z}$ to find Eqn. (vi):

Take the square-root of both sides of this inequality to obtain

and identify the righthand side as $\sqrt{Z}.$ Use the bound on $\sqrt{Z}$ to arrive at Eqn. (vii):

In conjunction with the previous inequality, this gives us the desired bound:

2.3 Running time Analysis

3 Proof of Theorem 5

To prove the first bound in the theorem (residual error analysis), we will use Lemma 24, which is the analog of Lemma 22 but for linear regression. Using this lemma, the proof of the first bound in Theorem 5 is similar to the proof of the first Frobenius norm bound of Theorem 4.

Manipulations similar to those used in proving the first Frobenius norm bound in Theorem 4 show that the lower bound on $r$ implies

The remaining estimates in the third inequality follow from applying Lemma 6 to obtain

Manipulations similar to those used in proving the first Frobenius norm bound in Theorem 4 show that the latter bound implies

Combining (27) with the bound on $S$, we obtain

Taking the square-root to both sides and using the fact that $\sqrt{1+22\varepsilon}\leq 1+22\varepsilon$ gives the bound in the theorem,

The failure probability in the theorem follows by a union bound on all the probabilistic events involved in the proof.

so $\alpha=1-\sqrt{\varepsilon}$ satisfies Eqn. (21). In the proof of the residual error bound in this theorem, we showed that

so $\beta=4\varepsilon$ satisfies Eqn. (22). With these choices of $\alpha$ and $\beta,$ Eqn. (23) in Lemma 25 gives the claimed forward error bound.

Experiments

Let $n=1024$ and consider the following three test matrices:

where $\bm{e}_{i}\in\R^{n+1}$ are the standard basis vectors.

2 Empirical comparison of the SRHT and Gaussian algorithms

3 Empirical evaluation of our error bounds

are used to construct the SRHT approximations, then the spectral norm residual error is no larger than

To bring perspective to this discussion, consider that even if one limits consideration to deterministic algorithms, the known error bounds for the Gu-Eisenstat rank-revealing QR—a popular and widely used algorithm for low-rank approximation—are quite pessimistic and do not reflect the excellent accuracy that is seen in practice . Regardless, we do not advocate using these approximation schemes for applications in which highly accurate low-rank approximations are needed. Rather, Theorem 4 and our numerical experiments suggest that they are appropriate in situations where one is willing to trade some accuracy for a gain in computational efficiency.

Acknowledgements

We would like to thank Joel Tropp and Mark Tygert for the initial suggestion that we attempt to sharpen the analysis of the SHRT low-rank approximation algorithm and for fruitful conversations on our approach. We are also grateful to an anonymous reviewer for pointing out the value in interpreting Lemma 16 as a relative error bound and to Malik Magdon-Ismail for providing the proof of Lemma 5.3.

Christos Boutsidis acknowledges the support from XDATA program of the Defense Advanced Research Projects Agency (DARPA), administered through Air Force Research Laboratory contract FA8750-12-C-0323. Alex Gittens was supported by ONR awards N00014-08-1-0883 and N00014-11-1002, AFOSR award FA9550-09-1-0643, DARPA award N66001-08-1-2065, and a Sloan Research Fellowship rewarded to Joel Tropp.

References