Orthogonal Random Features

Introduction

Kernel methods are widely used in nonlinear learning , but they are computationally expensive for large datasets. Kernel approximation is a powerful technique to make kernel methods scalable, by mapping input features into a new space where dot products approximate the kernel well . With accurate kernel approximation, efficient linear classifiers can be trained in the transformed space while retaining the expressive power of nonlinear methods .

This means that one can treat $p(\mathbf{w})$ as a density function and use Monte-Carlo sampling to derive the following nonlinear map for a real-valued kernel:

where $\mathbf{w}_{i}$ is sampled i.i.d. from a probability distribution with density $p(\mathbf{w})$. Let \mathbf{W}=\big{[}\mathbf{w}_{1},\cdots,\mathbf{w}_{D}\big{]}^{T}. The linear transformation $\mathbf{W}\mathbf{x}$ is central to the above computation since,

The choice of matrix $\mathbf{W}$ determines how well the estimated kernel converges to the actual kernel;

The computation of $\mathbf{W}\mathbf{x}$ has space and time costs of $\mathcal{O}(Dd)$. This is expensive for high-dimensional data, especially since $D$ is often required to be larger than $d$ to achieve low approximation error.

In this work, we address both of the above issues. We first show an intriguing discovery (Figure 1(c)): by enforcing orthogonality on the rows of $\mathbf{W}$, the kernel approximation error can be significantly reduced. We call this method Orthogonal Random Features (ORF). Section 3 describes the method and provides theoretical explanation for the improved performance.

Since both generating a $d\times d$ orthogonal matrix ($\mathcal{O}(d^{3})$ time and $\mathcal{O}(d^{2})$ space) and computing the transformation ($\mathcal{O}(d^{2})$ time and space) are prohibitively expensive for high-dimensional data, we further propose Structured Orthogonal Random Features (SORF) in Section 4. The idea is to replace random orthogonal matrices by a class of special structured matrices consisting of products of binary diagonal matrices and Walsh-Hadamard matrices. SORF has fast computation time, $\mathcal{O}(D\log d)$, and almost no extra memory cost (with efficient in-place implementation). We show extensive experiments in Section 5. We also provide theoretical discussions in Section 6 of applying the structured matrices in a broader range of applications where random Gaussian matrix is used.

Related Works

Explicit nonlinear random feature maps have been constructed for many types of kernels, such as intersection kernels , generalized RBF kernels , skewed multiplicative histogram kernels , additive kernels , and polynomial kernels . In this paper, we focus on approximating Gaussian kernels following the seminal Random Fourier Features (RFF) framework , which has been extensively studied both theoretically and empirically .

Key to the RFF technique is Monte-Carlo sampling. It is well known that the convergence of Monte-Carlo can be largely improved by carefully choosing a deterministic sequence instead of random samples . Following this line of reasoning, Yang et al. proposed to use low-displacement rank sequences in RFF. Yu et al. studied optimizing the sequences in a data-dependent fashion to achieve more compact maps. In contrast to the above works, this paper is motivated by an intriguing new discovery that using orthogonal random samples provides much faster convergence. Compared to , the proposed SORF method achieves both lower kernel approximation error and greatly reduced computation and memory costs. Furthermore, unlike , the results in this paper are data independent.

Structured matrices have been used for speeding up dimensionality reduction , binary embedding , deep neural networks and kernel approximation . For the kernel approximation works, in particular, the “structured randomness” leads to a minor loss of accuracy, but allows faster computation since the structured matrices enable the use of FFT-like algorithms. Furthermore, these matrices provide substantial model compression since they require subquadratic (usually only linear) space. In comparison with the above works, our proposed methods SORF and ORF are more effective than RFF. In particular SORF demonstrates both lower approximation error and better efficiency than RFF. Table 1 compares the space and time costs of different techniques.

Orthogonal Random Features

Our goal is to approximate a Gaussian kernel of the form

Recall that the linear transformation matrix of RFF can be written as

The idea of Orthogonal Random Features (ORF) is to impose orthogonality on the matrix on the linear transformation matrix $\mathbf{G}$. Note that one cannot achieve unbiased kernel estimation by simply replacing $\mathbf{G}$ by an orthogonal matrix, since the norms of the rows of $\mathbf{G}$ follow the $\chi$-distribution, while rows of an orthogonal matrix have the unit norm. The linear transformation matrix of ORF has the following form

Denote the approximate kernel based on the above $\mathbf{W}_{\textrm{ORF}}$ as $K_{\text{ORF}}(\mathbf{x},\mathbf{y})$. The following shows that $K_{\text{ORF}}(\mathbf{x},\mathbf{y})$ is an unbiased estimator of the kernel, and it has lower variance in comparison to RFF.

$K_{\text{ORF}}(\mathbf{x},\mathbf{y})$ is an unbiased estimator of the Gaussian kernel, i.e.,

Let $D\leq d$, and $z=||\mathbf{x}-\mathbf{y}||/\sigma$. There exists a function $f$ such that for all $z$, the variance of $K_{\text{ORF}}(\mathbf{x},\mathbf{y})$ is bounded by

We now show a proof sketch of the variance. Suppose, $a_{i}=\cos(\mathbf{w}^{T}_{i}\mathbf{z})$.

(Appendix A.3) There is a function $f$ such that for any $z$,

Therefore, for a large $d$, and $D\leq d$, the ratio of the variance of ORF and RFF is

Figure 2(a) shows the ratio of the variance of ORF to that of RFF when $D=d$ and $d$ is large. First notice that this ratio is always smaller than 1, and hence ORF always provides improvement over the conventional RFF. Interestingly, we gain significantly for small values of $z$. In fact, when $z\rightarrow 0$ and $d\rightarrow\infty$, the ratio is roughly $z^{2}$ (note $e^{x}\approx 1+x$ when $x\rightarrow 0$), and ORF exhibits infinitely lower error relative to RFF. Figure 2(b) shows empirical simulations of this ratio. We can see that the variance ratio is close to that of $d=\infty$ (3), even when $d=32$, a fairly low-dimensional setting in real-world cases.

Recall that $z=||\mathbf{x}-\mathbf{y}||/\sigma$. This means that ORF preserves the kernel value especially well for data points that are close, thereby retaining the local structure of the dataset. Furthermore, empirically $\sigma$ is typically not set too small in order to prevent overfitting—a common rule of thumb is to set $\sigma$ to be the average distance of 50th-nearest neighbors in a dataset. In Figure 2(c), we plot the distribution of $z$ for several datasets with this choice of $\sigma$. These distributions are all concentrated in the regime where ORF yields substantial variance reduction.

The above analysis is under the assumption that $D\leq d$. Empirically, for RFF, $D$ needs to be larger than $d$ in order to achieve low approximation error. In that case, we independently generate and apply the transformation (2) multiple times. The next lemma bounds the variance for this case.

Let $D=m\cdot d$, for an integer $m$ and $z=||\mathbf{x}-\mathbf{y}||/\sigma$. There exists a function $f$ such that for all $z$, the variance of $K_{\text{ORF}}(\mathbf{x},\mathbf{y})$ is bounded by

Structured Orthogonal Random Features

In the previous section, we presented Orthogonal Random Features (ORF) and provided a theoretical explanation for their effectiveness. Since generating orthogonal matrices in high dimensions can be expensive, here we propose a fast version of ORF by imposing structure on the orthogonal matrices. This method can provide drastic memory and time savings with minimal compromise on kernel approximation quality. Note that the previous works on fast kernel approximation using structured matrices do not use structured orthogonal matrices .

Let us first introduce a simplified version of ORF: replace $\mathbf{S}$ in (2) by a scalar $\sqrt{d}$. Let us call this method ORF´. The transformation matrix thus has the following form:

(Appendix B) Let $K_{\text{ORF}^{\prime}}(\mathbf{x},\mathbf{y})$ be the approximate kernel computed with linear transformation matrix (4). Let $D\leq d$ and $z=||\mathbf{x}-\mathbf{y}||/\sigma$. There exists a function $f$ such that the bias of $K_{\text{ORF}^{\prime}}(\mathbf{x},\mathbf{y})$ satisfies

The above implies that when $d$ is large $K_{\text{ORF}^{\prime}}(\mathbf{x},\mathbf{y})$ is a good estimation of the kernel with low variance. Figure 3(a) shows that even for relatively small $d$, the estimation is almost unbiased. Figure 3(c) shows that when $d\geq 32$, the variance ratio is very close to that of $d=\infty$. We find empirically that ORF´also provides very similar MSE in comparison with ORF in real-world datasets.

We now introduce Structured Orthogonal Random Features (SORF). It replaces the random orthogonal matrix $\mathbf{Q}$ of ORF´in (4) by a special type of structured matrix $\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}\mathbf{H}\mathbf{D}_{3}$:

Computing $\mathbf{W}_{\text{SORF}}\mathbf{x}$ has the time cost $\mathcal{O}(d\log d)$, since multiplication with $\mathbf{D}$ takes $\mathcal{O}(d)$ time and multiplication with $\mathbf{H}$ takes $\mathcal{O}(d\log d)$ time using fast Hadamard transformation. The computation of SORF can also be carried out with almost no extra memory due to the fact that both sign flipping and the Walsh-Hadamard transformation can be efficiently implemented as in-place operations .

Figures 3(b)(d) show the bias and variance of SORF. Note that although the curves for small $d$ are different from those of ORF, when $d$ is large ($d>32$ in practice), the kernel estimation is almost unbiased, and the variance ratio converges to that of ORF. In other words, it is clear that SORF can provide almost identical kernel approximation quality as that of ORF. This is also confirmed by the experiments in Section 5. In Section 6, we provide theoretical discussions to show that the structure of (5) can also be generally applied to many scenarios where random Gaussian matrices are used.

Experiments

Kernel Approximation. We first show kernel approximation performance on six datasets. The input feature dimension $d$ is set to be power of 2 by padding zeros or subsampling. Figure 4 compares the mean squared error (MSE) of all methods. For fixed $D$, the kernel approximation MSE exhibits the following ordering:

By imposing orthogonality on the linear transformation matrix, Orthogonal Random Features (ORF) achieves significantly lower approximation error than Random Fourier Features (RFF). The Structured Orthogonal Random Features (SORF) have almost identical MSE to that of ORF. All other fast kernel approximation methods, such as circulant and FastFood have higher MSE. We also include DigitalNet, the best performing method among Quasi-Monte Carlo techniques . Its MSE is lower than that of RFF, but still higher than that of ORF and SORF. The order of time cost for a fixed $D$ is

Remarkably, SORF has both better computational efficiency and higher kernel approximation quality compared to other methods.

We also apply ORF and SORF on classification tasks. Table 2 shows classification accuracy for different kernel approximation techniques with a (linear) SVM classifier. SORF is competitive with or better than RFF, and has greatly reduced time and space costs.

Simplifying SORF. The SORF transformation consists of three Hadamard-Diagonal blocks. A natural question is whether using fewer computations and randomness can achieve similar empirical performance. Figure 5(c) shows that reducing the number of blocks to two (HDHD) provides similar performance, while reducing to one block (HD) leads to large error.

Analysis and General Applicability of the Hadamard-Diagonal Structure

We provide theoretical discussions of SORF in this section. We first show that for large $d$, SORF is an unbiased estimator of the Gaussian kernel.

(Appendix C) Let $K_{\text{SORF}}(\mathbf{x},\mathbf{y})$ be the approximate kernel computed with linear transformation matrix $\sqrt{d}\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}\mathbf{H}\mathbf{D}_{3}$. Let $z=||\mathbf{x}-\mathbf{y}||/\sigma$. Then

Even though SORF is nearly-unbiased, proving tight variance and concentration guarantees similar to ORF remains an open question. The following discussion provides a sketch in that direction. We first show a lemma of RFF.

Let $\mathbf{W}$ be a random Gaussian matrix as in RFF, for a given $\mathbf{z}$, the distribution of $\mathbf{W}\mathbf{z}$ is $N(0,||z||_{2}\mathbf{I}_{d})$.

Note that $\mathbf{W}\mathbf{z}$ in RFF can be written as $\mathbf{R}\mathbf{g}$, where $\mathbf{R}$ is a scaled orthogonal matrix such that each row has norm $||z||_{2}$ and $\mathbf{g}$ is distributed according to $N(0,\mathbf{I}_{d})$. Hence the distribution of $\mathbf{R}\mathbf{g}$ is $N(0,||z||_{2}\mathbf{I}_{d})$, identical to $\mathbf{W}\mathbf{z}$. The concentration results of RFF use the fact that the projections of a Gaussian vector $\mathbf{g}$ onto orthogonal directions $\mathbf{R}$ are independent.

The above result can also be applied to settings not limited to kernel approximation. In the appendix, we show empirically that the same scheme can be successfully applied to angle estimation where the nonlinear map $f$ is a non-smooth $\operatornamewithlimits{sign}(\cdot)$ function . We note that the $\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}\mathbf{H}\mathbf{D}_{3}$ structure has also been recently used in fast cross-polytope LSH .

Conclusions

We have demonstrated that imposing orthogonality on the transformation matrix can greatly reduce the kernel approximation MSE of Random Fourier Features when approximating Gaussian kernels. We further proposed a type of structured orthogonal matrices with substantially lower computation and memory cost. We provided theoretical insights indicating that the Hadamard-Diagonal block structure can be generally used to replace random Gaussian matrices in a broader range of applications. Our method can also be generalized to other types of kernels such as general shift-invariant kernels and polynomial kernels based on Schoenberg’s characterization as in .

References

Appendix A Variance Reduction via Orthogonal Random Features

Let $\mathbf{z}=\frac{\mathbf{x}-\mathbf{y}}{\sigma}$, and $z=||\mathbf{z}||$. For a vector $\mathbf{y}$, let $y(i)$ denote its $i^{\text{th}}$ coordinate. Let $n!!$ be the double factorial of $n$, i.e., the product of every number from n to 1 that has the same parity as n.

A.2 Proof of Lemma 1

Let $\mathbf{z}=(\mathbf{x}-\mathbf{y})/\sigma$. Recall that in RFF, we compute the Kernel approximation as

where each $\mathbf{w}_{i}$ is a $d$ dimensional vector distributed $N(0,I_{d})$. Let $\mathbf{w}$ be a $d$ dimensional vector distributed $N(0,I_{d})$. By Bochner’s theorem,

and hence RFF yields an unbiased estimate.

We now compute the variance of RFF approximation. Observe that

If we take $D$ such independent random variables $\mathbf{w}_{1},\mathbf{w}_{2},\ldots\mathbf{w}_{D}$, since variance of the sum is sum of variances,

A.3 Proof of Lemma 2

For a set of non-negative values $\alpha_{1},\alpha_{2},\ldots\alpha_{k}$ and $\beta_{1},\beta_{2},\ldots\beta_{k}$ such that for all $i$, $\beta_{i}\leq\alpha_{i}$,

Combining the above two equations results in the first part of the lemma. For the second part observe that

Combining the above two equations yields the second part of the lemma. ∎

Since the problem is rotation invariant, instead of projecting a vector $\mathbf{z}$ onto a randomly chosen two orthogonal vectors $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$, we can choose a vector $\mathbf{y}$ that is uniformly distributed on a sphere of radius $z$ and project it on to the first two dimensions. Thus,

The $k^{\text{th}}$ term in the Taylor’s series expansion of sum of above two terms is

A way to compute a uniformly distributed random variable on a sphere with radius $z$ is to generate $d$ independent random variables $\mathbf{x}=(x(1),x(2),\ldots,x(d))$ each distributed $N(0,1)$ and setting $y(i)=zx(i)/||\mathbf{x}||$. Hence,

$(a)$ follows from linearity of expectation and the observation above. $(b)$ follows from the independence of $s_{1}$, $s_{2}$, and $\mathbf{x}$. $(d)$ follows from substituting the moments of chi and Gaussian distributions. $(e)$ follows from numerical simplification. We now describe the reasoning behind $(c)$. Let $\mathbf{z}=\frac{\mathbf{x}||\mathbf{y}||}{||\mathbf{x}||}$, where $\mathbf{y}$ and $\mathbf{x}$ are independent $N(0,I_{d})$ random variables. By the properties of the Gaussian random variables $\mathbf{z}$ is also a $N(0,I_{d})$ random variable. Thus,

and hence $(c)$. Substituting the above equation in the cosine expansion, we get that the expectation is

where $|c_{k,d}|\leq\frac{k^{4}}{4d^{2}}+\frac{k^{2}(k-1)}{2d^{3}}$. Thus,

Appendix B Proof of Theorem 2

The proof of the theorem is similar to that of Lemma 2 and we outline some key steps. We first bound the bias in Lemma 5 and then the variance in Lemma 6.

If $\mathbf{w}=\sqrt{d}\mathbf{y}$, where $y$ is distributed uniformly on a unit sphere, then

Without loss of generality, we can assume $\mathbf{z}$ is along the first coordinate and hence $\mathbf{w}^{T}\mathbf{z}=\sqrt{d}zy(1)$. A way to compute a uniformly distributed random variable on a sphere with radius $z$ is to generate $d$ independent random variables $\mathbf{x}=(x(1),x(2),\ldots,x(d))$ each distributed $N(0,1)$ and setting $y(i)=zx(i)/||\mathbf{x}||$. hence,

The $k^{\text{th}}$ term in the Taylor’s series expansion of cosine in the above equation is

Similar to the proof of Lemma 2, it can be shown that the expectation of this term is

where $|c^{\prime}_{k,d}|\leq\left(\frac{k(k-1)}{d}\right)^{2}$. Hence,

Let $D\leq d$. If $\mathbf{W}=\sqrt{d}\mathbf{Q}$, where $\mathbf{Q}$ is a uniformly chosen random rotation, then

Let $a_{i}=\cos(\mathbf{w}^{T}_{i}\mathbf{z})$. Expanding the variance we have,

For the first term, rewriting $\cos^{2}(\mathbf{w}^{T}\mathbf{z})=\frac{1+\cos(2\mathbf{w}^{T}\mathbf{z})}{2}$, similar to the proof of Lemma 5 it can be shown that

Second term can be bounded similar to Lemma 2 and here we just sketch an outline. Similar to the proof of Lemma 2, the variance boils down to computing the expectation of $\cos(\mathbf{w}^{T}_{1}\mathbf{z}+\mathbf{w}^{T}_{2}\mathbf{z})$. Using Lemma 4 and summing Taylor’s series we get

Substituting the above bound and the expectation from Lemma 5, we get

Appendix C Proof of Theorem 3

The proof follows from the following two technical lemmas.

Let $z^{\prime}$ be distributed according to $N(0,||\mathbf{x}||^{2}_{2})$ and $y^{\prime}=\sum^{d}_{i=1}x(i)d_{i}$, where $d_{i}$s are independent Rademacher random variables. For any function $g$ such that $|g^{\prime}|\leq 1$ and $|g|\leq 1$,

Let $z=z^{\prime}/||\mathbf{x}||_{2}$, $y=y^{\prime}/||\mathbf{x}||_{2}$, and $h(x)=g(||\mathbf{x}||_{2}x)$, for all $x$. Hence $h(z)=g(z^{\prime})$ and $h(y)=g(y^{\prime})$. By a lemma due to Stein ,

We now bound the term on the right hand side by classic Stein-type arguments.

Let $y_{i}=y-\frac{x(i)d_{i}}{||\mathbf{x}||_{2}}$. Observe that

where the first equality follows from the fact that $y_{i}$ and $d_{i}$ are independent and $d_{i}$ has zero mean. By Taylor series approximation, the first term is bounded by

Combining the above four equations, we get

Combining the above two equations, we get

Substituting the bound on the second moment of $f$ yields the result. ∎

Let $\mathbf{G}$ be a random matrix with i.i.d. $N(0,1)$ entries as before. Using the above lemma we show that $\sqrt{d}\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}$ behaves like $\mathbf{G}$ while computing the bias.

For a given $\mathbf{x}$, let $\mathbf{z}=\mathbf{G}\mathbf{x}$ and $\mathbf{y}=\sqrt{d}\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}\mathbf{x}$. For any function $g$ such that $|g^{\prime}|\leq 1$ and $|g|\leq 1$,

Let $\mathbf{u}=\mathbf{H}\mathbf{D}_{2}\mathbf{x}$. Then for every $i$, $y(i)=\sum_{j}H(i,j)D_{2}(j)u(j)$. Hence by Lemma 7, we can relate expectation under $y$ to the expectation under Gaussian distribution:

where the last equality follows from the fact that $\mathbf{H}\mathbf{D}_{2}$ does not change rotation and $(a)$ follows from the law of total expectation. By Cauchy-Schwartz inequality, for each $i$

Summing over all the indices yields the lemma. ∎

Theorem 3 follows from the Bochner’s theorem and the fact that $\cos(\cdot)$ satisfies requirements for the above lemma. We note that Theorem 3 holds for the matrix $\sqrt{d}\mathbf{H}\mathbf{D}_{1}\mathbf{H}\mathbf{D}_{2}$ itself and the third component $\mathbf{H}\mathbf{D}_{3}$ is not necessary to bound the bias.

Appendix D Proof of Theorem 4

To prove Theorem 4, we use the Hanson-Wright Inequality.

for some universal positive constant $c>0$.

For a vector $\mathbf{u}$, let $\text{diag}(\mathbf{u})$ denote the diagonal matrix whose entries correspond to the entries of $\mathbf{u}$. For a diagonal matrix $\mathbf{D}$, let $\text{vec}(\mathbf{D})$ denote the vector corresponding to the diagonal entries of $\mathbf{D}$. Let $\mathbf{v}=\mathbf{H}\mathbf{D}_{3}\mathbf{z}$ and $\mathbf{u}=\mathbf{H}\mathbf{D}_{2}\mathbf{v}=\mathbf{H}\text{diag}(\mathbf{v})\text{vec}(\mathbf{D}_{2})$. Observe that

is small. We first show that the expectation of the above quantity is and then use the Hanson-Wright inequality to prove concentration. Let $\mathbf{A}$ be a diagonal matrix with $k^{th}$ entry being $\sqrt{d}H(i,k)H(j,k)$. The above equation can be rewriten as

Observe that the $(l,l)$ entry of the $\mathbf{H}^{T}\mathbf{A}\mathbf{H}$ is

where the last equality follows from observing that the rows of $\mathbf{H}$ are orthogonal to each other. Together with the fact that elements of $\mathbf{D}_{2}$ are independent of each other, we get

To prove the concentration result, observe that the entries of $\text{vec}(\mathbf{D}_{2})$ are independent and sub-Gaussian, and hence we can use the Hanson-Wright inequality. To this end, we bound the Frobenius and the spectral norm of the underlying matrix. For the Frobenius norm, observe that

where $(a)$ follows by observing that each $\text{diag}(\mathbf{v})$ changes the Frobenius norm by at most $||\mathbf{v}||^{2}_{\infty}$, $(b)$ follows from the fact that $\mathbf{H}$ does not change the Frobenius norm, and $(c)$ follows by substituting $\mathbf{A}$.

where $(a)$ follows by observing that each $\text{diag}(\mathbf{v})$ changes the spectral norm by at most $||\mathbf{v}||_{\infty}$, $(b)$ follows from the fact that rotation does not change the spectral norm, and $(c)$ follows by substituting $\mathbf{A}$. Since $\mathbf{v}=\mathbf{H}\mathbf{D}_{3}\mathbf{z}$, by McDiarmid’s inequality, it can be shown that with probability $\geq 1-2de^{-d\epsilon^{2}/2}$, $||\mathbf{v}||_{\infty}\leq\epsilon||\mathbf{z}||_{2}$. Hence, by the Hanson-Wright inequality, we get

where $c$ is a constant. Choosing $\epsilon=({t/d})^{1/3}$ results in the theorem. ∎

Appendix E Discrete Hadamard-Diagonal Structure in Binary Embedding

We compare random projection based Locality Sensitive Hashing (LSH) , Circulant Binary Embedding (CBE) and Kronecker Binary Embedding (KBE) . We closely follow the experimental settings of . We choose to compare with because it proposed to use another type of structured random orthogonal matrix (Kronecker product of orthogonal matrices). As shown in Figure 6, our result (HDHDHD) provides higher recall and lower angular MSE in comparison with other methods.