Improved analysis of the subsampled randomized Hadamard transform

Introduction

Dimension reduction is an elegant idea from computer science that has found applications in numerical linear algebra. Here is the basic concept: it is often more efficient to solve a computational problem presented in a high-dimensional space if we first transport the problem instance to a lower-dimensional space while preserving its essential structure. Researchers have found that randomness provides an extraordinarily effective way to construct these dimension-reduction maps. This approach is usually traced to the celebrated paper of Johnson and Lindenstrauss [JL84].

In randomized algorithms for matrix approximation, the goal of dimension reduction is to find a low-dimensional subspace that captures most of the action of an input matrix. One way to accomplish this task is to multiply the input matrix $\bm{A}$ by a relatively small dimension-reduction matrix $\bm{\Omega}$ to obtain $\bm{Y}=\bm{A\Omega}$. We then perform a QR factorization $\bm{Y}=\bm{QR}$ to identify the range of the reduced matrix $\bm{Y}$. For an appropriately designed dimension-reduction map, it can be shown that $\bm{A}\approx\bm{QQ}^{*}\bm{A}$ with high probability. In words, we can approximate the input matrix by compressing it to the range of the reduced matrix. It is then possible to compute standard matrix decompositions of $\bm{A}$ by manipulating the low-rank approximation $\bm{QQ}^{*}\bm{A}$. See the recent survey [HMT11] for a comprehensive treatment of these ideas and an extensive bibliography.

In linear algebra applications, the cost of multiplication by an unstructured random matrix $\bm{\Omega}$, such as a Gaussian matrix, can sometimes be prohibitive. In that case, we may prefer to draw the random matrix $\bm{\Omega}$ from a highly structured distribution that allows us to multiply $\bm{\Omega}$ into the input matrix substantially faster. Sarlós is credited with bringing structured dimension reduction to numerical linear algebra [Sar06]; see also [WLRT08].

The subsampled randomized Hadamard transform (SRHT) is a type of structured dimension-reduction map that is based on the Walsh–Hadamard matrix. We will prove that the SRHT preserves the geometry of an entire subspace of vectors, which is the essential ingredient required to show that the SRHT can be used in algorithms for randomized linear algebra. See the discussion in [HMT11, Sec. 11] for details.

The literature already contains a number of papers, including [AC09, Lib09, NDT09, HMT11, AL11, KW11], that study the behavior of the SRHT and related dimension reduction maps. The current treatment differs in several regards. Here, the main technical difficulties are addressed using a version of the matrix Chernoff inequality [AW02, Tro11]. As a consequence of the simple proof schema, we are able—for the first time—to obtain optimal constants in our bounds. This improvement can be valuable in numerical applications that require concrete performance guarantees.

This document is adapted from Appendix B of the technical report [HMT09], but much of the proof is new. This paper should be viewed as a codicil to the published version [HMT11] of the technical report, which uses the results presented here.

$\bm{D}$ is a random $n\times n$ diagonal matrix whose entries are independent random signs, i.e., random variables uniformly distributed on $\{\pm 1\}$;

$\bm{H}$ is an $n\times n$ Walsh–Hadamard matrix, scaled by $n^{-1/2}$ so it is an orthogonal matrix;

Our analysis relies on two basic properties of the Walsh–Hadamard matrix: The derived matrix $\bm{H}$ is orthogonal, and its entries all have magnitude $\pm\,n^{-1/2}$. Walsh–Hadamard matrices exist for each $n=2^{p}$ where $p=1,2,3,\dots$.

2. Intuition

The purpose of the matrix product $\bm{HD}$ is to flatten out input vectors before we sample. To see why it achieves this goal, fix a unit vector $\bm{x}$, and examine the first component of $\bm{HD}\bm{x}$.

where $h_{ij}$ are the components of the matrix $\bm{H}$. This random sum clearly has zero mean. Since the entries of $\bm{H}$ have magnitude $n^{-1/2}$, the variance of the sum is $n^{-1}$. Hoeffding’s inequality [Hoe63] shows that

In words, the magnitude of the first component of $\bm{HD}\bm{x}$ is typically about $n^{-1/2}$. The same argument applies to the remaining entries. Therefore, it is unlikely that any one of the $n$ components of $\bm{HD}\bm{x}$ is larger than $\sqrt{n^{-1}\log(n)}$.

This discussion suggests that the precise form of the SRHT is not particularly important. Indeed, we can replace $\bm{H}$ by any unitary matrix whose entries are uniformly small and which is equipped with a fast matrix–vector multiply. We can also draw the diagonal entries of $\bm{D}$ from other subgaussian distributions. These changes complicate the analysis somewhat.

3. Main Results

It turns out that an appropriately designed SRHT also preserves the geometry of an entire subspace of vectors. For this type of structured map, the embedding requires $k\log(k)$ dimensions because of certain phenomena connected with random sampling (Section 3.3). It also becomes necessary to use more sophisticated proof techniques. We have the following result.

The symbol $\sigma_{j}(\cdot)$ denotes the $j$th largest singular value of a matrix.

The proof of Theorem 1.3 appears below. This result replaces Theorem B.4 of [HMT09]. Some precedents include [RV07, Thm. 3.1] and [Tro08, Sec. 9]; see also [NDT09]. Earlier results have the same structure as Theorem 1.3, but the constants are either exorbitant or absent.

For large problems, we have been able to obtain optimal numerical constants. Suppose that $\iota$ is a sufficiently small positive number. If $k\gg\log(n)$, then sampling

coordinates is sufficient to ensure that $\bm{\Phi V}$ has constant condition number. See Theorem 3.2 for a more precise statement. The discussion in Section 3.3 indicates that there are cases where it does not suffice to draw $(1-\iota)\cdot k\log(k)$ samples.

Technical Background

In preparation for the main argument, we present our notation and some probability inequalities. These inequalities encapsulate all the difficulty in the proof.

A Rademacher random variable takes the values $\pm 1$ with equal probability. We reserve the letter $\varepsilon$ for a Rademacher variable, and we often write $\bm{\varepsilon}$ for a vector whose entries are independent Rademacher variables.

2. Probability Inequalities

We delegate the hard work to some probability inequalities that describe the large-deviation behavior of specific types of random variables. First, we describe a tail bound for a convex function of Rademacher variates. This result was established by Ledoux [Led96, Eqn. (1.9)]; see also [Led01, §5.2] for a discussion of concentration in product spaces.

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$,

Next, we present a matrix analog of the well-known Chernoff inequality. The proof is based on the matrix Laplace transform method proposed in an influential paper of Ahlswede–Winter [AW02]. We derive the result using more recent ideas [Tro11], which deliver an essential improvement on the earlier work. The other key tool is a method [GN10], ultimately due to Hoeffding [Hoe63], for transferring results from the model where we sample with replacement to the model where we sample without replacement.

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

We establish only the upper bound; the lower bound is established by applying a similar method to $-\sum\nolimits_{j}\bm{X}_{j}$. By homogeneity, take $B=1$. We use the matrix Laplace transform method [Tro11, Prop. 3.1] to bound the probability of a large deviation:

The Laplace transform bound (2.1) is essentially due to Ahlswede and Winter [AW02].

Lemma 5.8 of [Tro11] gives a semidefinite bound for the mgf of $\bm{Y}_{j}$. (See also [AW02, Thm. 19].)

On the right-hand side of this bound, we have exploited the fact that $\bm{Y}_{j}\sim\bm{Y}_{1}\sim\bm{X}_{1}$. Lemma 3.4 of [Tro11] allows us to use the mgf bound (2.3) to bound the mgf of the entire sum:

Substitute (2.4) into (2.2) and introduce the resulting inequality into the probability bound (2.1). The infimum is achieved at $\theta=\log(1+\delta)$. ∎

The SRHT Preserves Geometry

In this section, we establish a slightly more specific version of our main result, Theorem 1.3.

We present the main line of reasoning here, postponing the proofs of the lemmata. The first step is to show that the matrix $\bm{HD}$ equilibrates row norms.

Let $\bm{V}$ be an $n\times k$ matrix with orthonormal columns. Then $\bm{HDV}$ is an $n\times k$ matrix with orthonormal columns, and

When $k\gg\log n$, the second term in the bound is negligible, in which case the row norms are essentially as small as possible. On the other hand, when $k<\log n$, small sample effects can make some row norms large.

The next result states that randomly sampling rows from a matrix with orthonormal columns results in a well-conditioned matrix. The minimum size of the sample depends primarily on the row norms, and the sampling procedure is most efficient when the matrix has uniformly small rows. We state this result in detail because it may have independent interest.

Let $\bm{W}$ be an $n\times k$ matrix with orthonormal columns, and define the quantity M:=n\cdot\max_{j=1,\dots,n}{\bigl{\|}{\mathbf{e}_{j}^{*}\bm{W}}\bigr{\|}}^{2}. For a positive parameter $\alpha$, select the sample size

Define the matrix $\bm{W}=\bm{HDV}$. Lemma 3.3 with $\beta=k$ establishes that $\bm{W}$ is a matrix with orthonormal columns whose largest row norm is essentially as small as possible, so that

except with probability $k^{-1}$. Next, we apply Lemma 3.4 with $\alpha=4$, with $\delta=5/6$, and with $\eta=7/6$. A numerical reckoning shows that the additional probability of failure is at most $2k^{-1}$. Altogether, we obtain the following result. For

We can establish a slightly stronger sample bound in Theorem 3.2 by choosing the parameter $\iota=A/\log(k)$ for a sufficiently large number $A$. With this selection, however, we do not obtain a constant bound for the lower singular value.

2. Proofs of Supporting Lemmas

It remains to check that the underlying results are true. We begin with the claim that $\bm{HD}$ balances row norms.

because $\bm{D},\bm{H}$ are orthogonal matrices.

Fix a row index $j\in\{1,2,\dots,n\}$, and define the function

We have written $\bm{E}:=\operatorname{diag}(\mathbf{e}_{j}^{*}\bm{H})$ for the diagonal matrix constructed from the $j$th row of $\bm{H}$; observe that each entry of $\bm{E}$ has magnitude $n^{-1/2}$. The function $f$ is convex, and we quickly determine its Lipschitz constant:

We may use the function $f$ to study the variation of the rows norms of $\bm{HDV}$. Recall that $\bm{D}:=\operatorname{diag}(\bm{\varepsilon})$ for a Rademacher vector $\bm{\varepsilon}$, and consider the random variable

Apply the Rademacher tail bound, Proposition 2.1, with $t=\sqrt{8\log(\beta n)}$ to reach

This estimate holds for each row index $j=1,2,\dots,n$. Finally, take a union bound over these $n$ events to reach the advertised conclusion. ∎

Now we establish the result on row sampling.

Let $\bm{w}_{j}^{*}$ denote the $j$th row of $\bm{W}$, and define $M:=n\cdot\max_{j}\left\|{\bm{w}_{j}}\right\|^{2}$.

We can control the extreme singular values of the random matrix $\bm{R}_{T}\bm{W}$ by bounding the extreme eigenvalues of the $k\times k$ Gram matrix

The matrix Chernoff bound, Theorem 2.2, allows us to obtain large deviation bounds for the extreme eigenvalues of $\bm{Y}$. We just need to determine the parameters involved in the statement. First, note that

We easily compute the expectation of the first sample using the fact that the columns of $\bm{W}$ are orthonormal:

Simplify the formulae to complete the proof. ∎

3. Collecting Coupons

The logarithmic factor in these results is necessary, as we show by example. This discussion is extracted from [HMT11, Sec. 11].

Fix an integer $k$, and set $n=k^{2}$. Form an $n\times k$ orthonormal matrix $\bm{W}$ by regular decimation of the $n\times n$ identity matrix. More precisely, $\bm{W}$ is the matrix whose $j$th row has a unit entry in column $1+(j-1)/k$ when $j\equiv 1\pmod{k}$ and is zero otherwise. To see why this type of matrix is inconvenient, it is helpful to consider an auxiliary matrix $\bm{V}=\bm{HD}\bm{W}$. Observe that, up to scaling and modulation of rows, $\bm{V}$ consists of $k$ copies of a $k\times k$ Walsh–Hadamard transform stacked vertically.

References