Quantization Algorithms for Random Fourier Features

Xiaoyun Li, Ping Li

Introduction

In recent years, machine learning with extremely large-scale (and high-dimensional) datasets has become increasingly important given the rapid development of modern technologies. Many industrial applications involve massive data collected from a wide range of sources, e.g., internet and mobile devices. Designing efficient large-scale learning algorithms and feature engineering techniques, in terms of both speed and memory, has been an important topic in the machine learning & data mining community. The method of Random Projection (RP) is a popular strategy to deal with massive data, for example, for efficient data processing, computations, storage, or transmissions. The theoretical merit of RP is highlighted by the celebrated Johnson-Lindenstrauss Lemma (Johnson and Lindenstrauss, 1984), which states that with high probability the Euclidean distance between data points is approximately preserved in the projected space provided that the number of projections is sufficiently large. In the past two decades or so, RP has been used extensively in dimensionality reduction, approximate near neighbor search, compressed sensing, computational biology, etc. See some examples of relatively early works on RP (Dasgupta, 2000; Bingham and Mannila, 2001; Buhler, 2001; Achlioptas, 2003; Fern and Brodley, 2003; Datar et al., 2004; Candès et al., 2006; Donoho, 2006; Li et al., 2006; Freund et al., 2007; Li, 2007). In this paper, we continue the line of research on random projections and focus on studying quantization schemes for using random Fourier features (RFF), which are nonlinear transformations of random projections, to accurately approximate the (nonlinear) Gaussian kernel.

There are two major general issues with large-scale nonlinear kernel learning (not limited to the Gaussian kernel). Firstly, storing/materializing a kernel matrix for a dataset of $n$ samples would need $n^{2}$ entries, which may not be realistic even just for medium datasets (e.g., $n=10^{6}$ ). To avoid this problem, the entries of the kernel matrix are computed on the fly from the original dataset. This however will increase the computation time, plus storing the original high-dimensional dataset for on-demand distance computations can also be costly. Secondly, the training procedure for nonlinear kernel algorithms is also well-known to be expensive (Platt, 1998; Bottou et al., 2007). Therefore, it has been an active area of research to speed up kernel machines, and using various types of random projections has become popular.

2 Random Projections (RP) and Random Fourier Features (RFF)

This is the basic idea of using random projections to approximate inner product. See Li et al. (2006) for the theoretical analysis (such as exact variance calculations) of this approximation scheme.

We can also use random projections to approximate the (nonlinear) Gaussian kernel with an additional step. The Random Fourier Feature (RFF) (Rudin, 1990; Rahimi and Recht, 2007) is defined as

where $\tau\sim uniform(0,\ 2\pi)$ , the uniform random variable. Some basic probability calculations reveal that

In other words, the inner product between the RFFs of two data samples provides an unbiased estimate of the Gaussian kernel. The simulations need to be repeated for a sufficient number of times in order to obtain reliable estimates. That is, we generate $m$ independent RFFs using i.i.d. $w_{1},...,w_{m}$ and $\tau_{1},...,\tau_{m}$ , and approximate the kernel $K(u,v)$ by the following unbiased estimator:

where $F_{i}$ denotes the RFF generated by $w_{i},\tau_{i}$ . Furthermore, Li (2017b) showed that one can actually reduce the estimation variances by normalizing the RFFs.

In large-scale learning, using above estimator simply requires taking the inner product between the RFF vectors of $u$ and $v$ . Therefore, feeding the RFFs into a linear machine will approximate training a non-linear kernel machine, known as kernel linearization, which may significantly accelerate training and alleviate memory burden for storing the kernel matrix. This strategy has become popular in the literature, e.g., (Raginsky and Lazebnik, 2009; Yang et al., 2012; Affandi et al., 2013; Hernández-Lobato et al., 2014; Dai et al., 2014; Yen et al., 2014; Hsieh et al., 2014; Shah and Ghahramani, 2015; Chwialkowski et al., 2015; Richard et al., 2015; Sutherland and Schneider, 2015; Li, 2017b; Avron et al., 2017; Sun et al., 2018; Tompkins and Ramos, 2018; Li et al., 2020).

3 Quantized Random Projections (QRP)

One can further compress the projected data by quantization, into discrete integer values, or even binary values in the extreme case. The so-called quantized random projection (QRP) has found useful in many problems, e.g., theory, similarity search, quantized compressed sensing, classification and regression (Goemans and Williamson, 1995; Charikar, 2002; Datar et al., 2004; Zymnis et al., 2010; Jacques et al., 2013; Leng et al., 2014; Li et al., 2014; Li and Slawski, 2017; Slawski and Li, 2018; Li and Li, 2019b, a). The motivation is straightforward. If one can represent each RP (or RFF) using (e.g.,) 4 bits and still achieve similar accuracy as using the full-precision (e.g., 32 or 64 bits), it is then a substantial saving in storage space. Typically, savings in storage can directly translate into savings in data transmissions and subsequent computations. In addition to space (computation) savings, there is another motivation for QRP. That is, quantization also provides the capability of indexing due to the integer nature of quantized data, which can be used to build hash tables for approximate near neighbor search (Indyk and Motwani, 1998).

The simplest quantization scheme is the 1-bit (sign) random projections, including sign Gaussian random projections (Goemans and Williamson, 1995; Charikar, 2002) and sign Cauchy random projections (Li et al., 2013) (for approximating the $\chi^{2}$ kernel). Basically, one only keeps the signs of projected data. Even though the 1-bit schemes appear to be overly crude and simplistic, in some cases 1-bit random projections can achieve better performance than full-precision RPs in similarity search and nearest neighbor classification tasks. Nevertheless, in general, one would need more than just 1-bit in order to achieve sufficient accuracy. For example, Li and Slawski (2017); Slawski and Li (2018); Li and Li (2019b) apply the (multi-bit) Lloyd-Max (LM) quantization (Max, 1960; Lloyd, 1982) on the projected data.

4 Summary of Contributions

Since each Lloyd-Max (LM) quantizer is associated with a specific random signal distribution, at the first glance, designing LM quantizers for the random Fourier features and the Gaussian kernel might appear challenging, due to the tuning parameter $\gamma$ , which is a crucial component of the Gaussian kernel. Initially, one might expect that a different LM quantizer would be needed for a different $\gamma$ value. In this paper, our contribution begins with an interesting finding that the marginal distribution of the RFF is actually free of the parameter $\gamma$ . This result greatly simplifies the design of LM quantization schemes for the RFF, because only one quantizer would be needed for all $\gamma$ values. Once we have derived the marginal distribution of the RFF, we incorporate the idea of distortion optimal quantization theory to nonlinear random feature compression by providing a thorough study on the theoretical properties and practical performance. Extensive simulations and machine learning experiments validate the effectiveness of the proposed LM quantization schemes for the RFF.

The Probability Distributions of RFF

We start the introduction to our proposed method by providing analysis on the probability distribution of RFF (4), which is key to the design of quantization schemes in Section 3. First, we introduce some notations.

Throughout the paper, we will use the following two definitions for $\phi_{\sigma}(t)$ and $\Phi(t)$ :

That is, $\Phi(t)$ is the cumulative distribution function (cdf) of the standard normal $N(0,1)$ and $\phi_{\sigma}(t)$ is the probability density function (pdf) of $N(0,\sigma^{2})$ .

We first consider the marginal distributions of the RFF, which serve as the foundation of our proposed quantization schemes. The following Lemma is a result of the convolution of normal and uniform distributions.

Suppose $X\sim N(0,1)$ and $\tau\sim uniform(0,2\pi)$ are independent, $\gamma>0$ . Then

In the following, we formally give the distribution of the RFF.

Let $X\sim N(0,1)$ , $\tau\sim uniform(0,2\pi)$ be independent. Denote $Z=\cos(\gamma X+\tau)$ , and $Z_{2}=\cos^{2}(\gamma X+\tau)$ . We have the probability density functions

for any $\gamma>0$ . In particular, $Z\overset{d}{\sim}\cos(\tau)$ in distribution.

The density plots can be found in Figure 1. Theorem 2.2 says that for any kernel parameter $\gamma$ , the (unscaled) RFF follows the same distribution as the cosine of the uniform noise itself. Intuitively, this is because cosine is a $2\pi$ -periodic function and normal distribution is symmetric. As will be introduced in Section 3, each Lloyd-Max (LM) quantizer is associated with a signal distribution. We will characterize two LM-type quantizers w.r.t. density (6) and (7), respectively. This interesting result is favorable for our purpose as it implies we only need to construct one LM quantizer, which covers all the Gaussian kernels with different $\gamma$ value. Thus, the design of LM quantizer for RFF is convenient.

In Theorem 2.2 we consider $X\sim N(0,1)$ because we assume data samples are normalized for conciseness. It is easy to see that this result also holds without data normalization (i.e., $X$ is Gaussian with arbitrary variance) since we can offset the variance of $X$ by altering $\gamma$ . Therefore, Theorem 2.2 is a universal result implying that the LM quantizer also works without data normalization.

2 Joint Distribution

In the sequel, we analyze the joint distribution of RFFs of two data samples with correlation $\rho$ ’s. The joint distribution will play an important role in later theoretical analysis. The following Lemma 2.3 leads to the desired result presented in Theorem 2.4.

Denote $z_{x}=\gamma X+\tau$ , $z_{y}=\gamma Y+\tau$ with $(X,Y)\sim N\big{(}0,\begin{pmatrix}1&\rho\\ \rho&1\end{pmatrix}\big{)}$ , $\tau\sim uniform(0,2\pi)$ . We have the joint distribution

Denote $z_{x}=\cos(\gamma X+\tau)$ , $z_{y}=\cos(\gamma Y+\tau)$ where $X,Y,\tau$ are the same as Lemma 2.3. Then we have the joint density function for $(z_{x},z_{y})\in^{2}$ ,

where $a_{x}^{*}=\cos^{-1}(z_{x}),a_{y}^{*}=\cos^{-1}(z_{y})$ . In addition, $\big{(}\sin(\gamma X+\tau),\ \sin(\gamma Y+\tau)\big{)}$ follows the same distribution.

In Figure 2, we plot the joint density at several $\gamma$ values. We conclude several properties of the joint distribution. Firstly, it is obvious that $z_{x}$ and $z_{y}$ are exchangeable, i.e., $f(Z_{x},Z_{y})=f(Z_{y},Z_{x})$ . Secondly, it is symmetric which means $f(Z_{x},Z_{y})=f(-Z_{x},-Z_{y})$ . Moreover, we have the following important result, which is helpful for our analysis on the monotonicity and variance of quantized kernel estimators in Section 4.

Let the density function $f$ be defined as Theorem 2.4. If $\sqrt{2(1-\rho)}\gamma\leq\pi$ , then $f(z_{x},z_{y})>f(z_{x},-z_{y})$ for $\forall(z_{x},z_{y})\in(0,1]^{2}$ or $(z_{x},z_{y})\in[-1,0)^{2}$ .

In Proposition 2.5, the quantity $\sqrt{2(1-\rho)}\gamma$ will be reduced if we either increase $\rho$ or decrease $\gamma$ . In Figure 2, we see how this term characterizes the joint density of RFF. In particular, smaller $\sqrt{2(1-\rho)}\gamma$ reinforces the dependency between $z_{x}$ and $z_{y}$ . The density around $(1,1)$ and $(-1,-1)$ reaches the highest when $\rho=0.9$ and $\gamma=1$ . As $\rho$ decreases or $\gamma$ increases, the density is “flattened”.

Quantization Schemes for RFF

Quantization is, to a good extent, an ancient topic in information theory and signal processing (Widrow and Kollár, 2008). On the other hand, many interesting research works appear in the literature even very recently, for achieving better efficiency in data storage, data transmission, computation, and energy consumption. Quantization for random projections has been heavily studied. In this paper, we focus on developing quantization schemes for random Fourier features, in particular, based on the Lloyd-Max (LM) framework.

This article https://www.eetimes.com/an-introduction-to-different-rounding-algorithms/, published in EETimes 2006, provides a good summary of common quantization (rounding) schemes: “rounding in decimal”, “round-toward-nearest”,“round-half-up”, “round-half-down”,“round-half-even”, “Round-half-odd”, “round-alternate”, “round-random”, “round-ceiling”, “round-floor”, “round-toward-zero”, “round-away-from-zero”, “round-up”, “round-down”, “truncation”. For random projections, quantization methods such as Datar et al. (2004); Li et al. (2014); Li (2017a) belong to those categories. One should expect that, as long as we use a sufficient number (such as 32 or 64) of bits, any reasonable quantization scheme should achieve a good accuracy.

In this paper, we mainly focus on quantization for RFF with a small number (such as 1, 2, 3, or 4) bits. We consider general multi-bit quantizers, with 1-bit quantization as a special case, for the cosine feature in RFF bounded in $ $. Here, “$ b $-bit” means the quantizer has$ 2^{b} $levels. For simplicity, we will denote$ z=\cos(\gamma w^{T}u+\tau)$ as the item to be quantized. In particular, we focus on two algorithms: “round-random” (which we refer to as the “stochastic quantization (StocQ)”) and “Lloyd-Max (LM) quantization”.

A $b$ -bit StocQ quantizer, also known as stochastic rounding, splits $ $into ($ 2^{b}-1 $) intervals, with consecutive borders$ -1=t_{0}<... $. Let$

2 Lloyd-Max (LM) Quantization

aiming to keep most amount of information of the original signal. For the signal distribution, we consider two variants. First, it is natural to set the target distribution as the distribution of RFF itself (6). Consequently, the first LM quantizer is subject to the distortion: {smallequation} LM-RFF: D_1≜∫_ ( z-Q(z))^2 1π1-z2dz. Conceptually, optimizing (3.2) minimizes the average difference between RFF $z$ and $Q(z)$ . Alternatively, we may also choose to address more on the high similarity region (which is more important sometimes). Denote $z_{u}$ and $z_{v}$ as the RFFs of $u$ and $v$ . As $\rho\rightarrow 1$ , we have $z_{x}\cdot z_{y}\rightarrow z_{x}^{2}$ , with density given by (7). Thus, our second variant, LM2-RFF, is designed to approximate $z^{2}$ by $Q(z)^{2}$ with distortion {smallequation} LM2-RFF: D_2≜∫_ ( z^2-Q(z)^2)^2 1π1-z2dz. Minimizing (3.2) and (3.2) leads to our proposed two LM-type quantizers for RFF compression in this paper. Note that, in Eq. (3.2), $\int_{0}^{1}\frac{1}{\pi\sqrt{1-z^{2}}}dz=\int_{0}^{1}\frac{1}{\pi\sqrt{1-z}}dz^{1/2}=\frac{1}{2}\int_{0}^{1}\frac{1}{\pi\sqrt{z-z^{2}}}dz=\frac{1}{2}\int_{0}^{1}f_{Z_{2}}(z)dz$ , where $f_{Z_{2}}$ is the density (7).

3 Optimization for LM Quantizers

To solve the introduced optimization problems, we exploit classical Lloyd’s algorithm, which alternatively updates two parameters until convergence. By e.g., Wu (1992), the algorithm converges to the globally optimal solution since the squared loss is convex and symmetric. The algorithm terminates when the total absolute change in borders and reconstruction levels in two consecutive iterations is smaller than a given threshold (e.g., $10^{-5}$ ). We provide the concrete steps for LM-RFF and LM2-RFF in Algorithm 1 and Algorithm 2, respectively. For LM-RFF, we see that the procedure is standard (exactly the same as above derivation). Denote $z$ as the unscaled RFF, $z=\cos(\gamma X+\tau)$ , $X\sim N(0,1)$ and $\tau\sim uniform(0,2\pi)$ . For LM2-RFF, recall that our objective is to minimize (by change of random variable)

In Figure 3, we plot the $2$ -bit LM-RFF and LM2-RFF quantizer as an example, along with the distortions of LM quantizers and uniform stochastic quantization (StocQ), with various number of bits. We see that both LM methods give non-uniform quantization borders and codes. LM2-RFF “expands” more towards two ends since it tries to approximate $z^{2}$ . From the distortion plots, we validate that LM-RFF provides smallest $D_{1}$ and LM2-RFF gives smallest $D_{2}$ .

As a final remark before ending this subsection, we note that LM quantization is more convenient and faster compared with StocQ in practical implementation, for quantizing the full-precision RFFs. While LM is a fixed quantization approach, StocQ requires generating an extra random number for each sketch and each data point. For large datasets, producing these additional random numbers might be rather slow in practice.

4 Quantized Kernel Estimators

with $z_{u,i}=\cos(w_{i}^{T}u+\tau_{i})$ and $z_{v,i}=\cos(w_{i}^{T}v+\tau_{i})$ the $i$ -th unscaled RFF of $u$ and $v$ , respectively. Moreover, for the proposed LM quantizers, we consider normalized estimator,

This estimator can also be conveniently used, as we only need to normalize the quantized RFFs (per data point) before learning. We will use $\hat{K}_{Q,(2)}$ and $\hat{K}_{n,Q,(2)}$ to denote the corresponding estimators using LM2-RFF quantization. We will analyze and compare the estimators using different quantization methods, theoretically and practically, in the remaining sections of the paper.

Theoretical Analysis

In this section, we first analyze the mean, variance, and monotonicity property of the quantized kernel estimators, then discuss kernel matrix approximation property based on some new evaluation metrics that can well align with the generalization performance. The proofs are deferred to Appendix C.

We start this section by analyzing the stochastic rounding method for RFF. In Zhang et al. (2019), the exact variance of the kernel estimator is not provided. In the following, we establish the precise variance calculation in Theorem 2.4, which is in fact a more general result on any symmetric stochastic quantizer.

which is always greater than $Var[\hat{K}]$ defined in (5).

The important take-away messages are: 1) the StocQ kernel estimator is unbiased of the Gaussian kernel; 2) the variance is always larger than full-precision RFF estimate. Further, we have the following result for 1-bit StocQ, which is a straightforward consequence of Theorem 4.1

With 1-bit, $Var[\hat{K}_{StocQ}]=4-K(u,v)^{2}$ .

2 LM Estimators

In this subsection, we study the moments of the proposed LM kernel estimators. Since our following results generalize to both LM-RFF and LM2-RFF, we will unify the notation as $Q$ to denote a LM-type quantizer. First, we have the following formulation of the mean estimate of LM quantized estimator (8) based on Chebyshev functional approximation.

Next, we provide an asymptotic analysis on the normalized quantized kernel estimate (9) under LM scheme.

Under same setting as Theorem 4.3, as $m\rightarrow\infty$ ,

Validation. We plot the empirical bias of LM-RFF against Observations 4.1 and 4.2 in Figure 4. As we see, the proposed surrogates for bias align with true biases very well when $\rho$ is not very close to $1$ . The biases shrink to as $b$ increases (e.g., negligibly $\mathcal{O}(10^{-3})$ with $b=4$ ). As $\rho\rightarrow 1$ , at some "disjoint point" the absolute biases have sharp drops and quickly converge to the theoretical values (red dots) given in Theorem 4.3 and 4.4. As $b$ or $\gamma$ increases, the “disjoint point” gets closer to $\rho=1$ .

3 Variance Comparisons

Figure 5 provides variance comparisons, where full-precision estimator variances are plotted for reference. As $b$ gets larger, the variances of LM-type quantized estimators converges to those of full-precision estimators. The variance of StocQ is significantly larger than RFF and LM quantization, especially when $b=1,2$ . This to a good extent explains why StocQ performs poorly in approximate kernel learning with low bits (Section 5).

Variance of debiased kernel estimates. As shown previously, LM estimators are slightly biased which brings theoretical challenges on finding a method to “properly” compare their variances. In this paper, we investigate the concept of “debiased variance”, which refers to the estimator variance after bias corrections.

Note that, the debiasing step is only for analytical purpose. Intuitively, Definition 4.1 is reasonable in that it compares the variation of different estimation procedures given that they have same mean. It is worth mentioning that, DB-variance is invariant of linear scaling, i.e., $c\hat{K}$ and $\hat{K}$ have same DB-variance for $c>0$ . Classical metrics for estimation quality, such as the Mean Squared Error (MSE), might be largely affected by such simple scaling. Note that, the DB-variance of all 1-bit estimators (both simple and normalized) from fixed quantizers are essentially identical. This can be easily verified by writing every 1-bit quantizer as $Q(z)=sign(z)\cdot C_{Q}$ for some $C_{Q}>0$ and substituting it into (8) and (9). Thus, we will focus on multi-bit quantizers (i.e., $b\geq 2$ ).

LM-RFF v.s. LM2-RFF. In Figure 6, we provide the DB-variance ratio of LM2-RFF estimator against that of LM-RFF estimator in the 2-bit case. (The observed pattern is the same for more bits.) For the simple kernel estimator, we see that in general LM-RFF has smaller DB-variance. Yet, the DB-variance of LM2-RFF sharply drops towards 0 and beats LM-RFF as $\rho\rightarrow 1$ , i.e., in high similarity region, which meets the goal of LM2-RFF quantizer design (to favor high similarity region). However, for normalized estimators, $\hat{K}_{n,Q}$ has consistently smaller DB-variance than $\hat{K}_{n,Q,(2)}$ .

Benefit of normalization. Next we prove the theoretical merit of normalizing RFFs, in terms of DB-variance.

Suppose $u,v$ are two samples with correlation $\rho$ . Let the simple and normalized kernel estimator, $\hat{K}_{Q}$ and $\hat{K}_{n,Q}$ , be defined as (8) and (9), respectively, where $Q$ is any LM-type quantizer. Assume $\gamma\leq\pi/\sqrt{2}$ . Then, $Var^{db}[\hat{K}_{n,Q}]\leq Var^{db}[\hat{K}_{Q}]$ on $\rho\in$ as $m\rightarrow\infty$ .

Theorem 4.5 says that when $\gamma\leq\pi/\sqrt{2}\approx 2.2$ , normalization is guaranteed to reduce the DB-variance at any $\rho\in$ . In Figure 7, we plot the DB-variance ratio of $\frac{Var^{db}[\hat{K}_{n,Q}(x,y)]}{Var^{db}[\hat{K}_{Q}(x,y)]}$ at multiple $\gamma$ and $b$ , for LM-RFF and LM2-RFF respectively. We corroborate the advantage of normalized estimates over simple estimators in terms of DB-variance (ratio always $<1$ ), especially with large $\rho$ .

4 Monotonicity of Mean Kernel Estimation

Next, in the following theorem, we extend the above result to discrete functions, which include our proposed LM quantizers as special cases.

Numerical Experiments

We conduct experiments with compressed RFFs on approximate kernel SVM (KSVM) classification and kernel ridge regression (KRR) tasks. Our results illustrate the effectiveness of large-scale kernel learning with highly compressed RFFs, highlighting the superior advantage of the proposed LM-RFF quantization. Moreover, we also propose and evaluate robust kernel approximation error metrics to consolidate our claims.

For this task, we use four popular public datasets from UCI repositoryhttps://archive.ics.uci.edu/ml/index.php (Dua and Graff, 2017) and ASU databasehttps://jundongl.github.io/scikit-feature/datasets.html (Li et al., 2016). All the data samples are pre-processed by instance normalization, and we randomly split each dataset into 60% for training and 40% for testing. For each task and each quantization method, the best test tuned accuracy is reported, averaged over 10 independent runs.

To compare the learning power of different compression schemes, we provide the test accuracy vs. number of RFFs in the left two columns of Figure 9, with $b=1,2$ . We observe: 1) LM-RFF substantially outperforms StocQ on all datasets with low bits. In particular, 1-bit StocQ performs very poorly, while 1-bit LM-RFF achieves similar accuracy as full-precision RFF; 2) On all datasets, LM-RFF with $b=2$ already approaches the accuracy of full-precision RFF with moderate $m\approx 4000$ , indicating the superior learning capacity of LM-RFF under deep feature compression. As expected, with larger $b$ , the performance of StocQ approaches that of LM-RFF. In particular, when $b=4$ , LM-RFF and StocQ perform similarly on those datasets.

To characterize the memory efficiencyFor simplicity, we mainly consider the memory cost for (quantized) RFF storage, which dominates in large-scale learning., note that under $b$ -bit compression, each data sample requires $m\times b$ bits in total for storage. If we assume that each full-precision RFF is represented by 32 bits, then the storage cost per sample for full-precision RFF is $32m$ . This allows us to plot the test accuracy against the total memory cost per sample, as shown in the right two columns of Figure 9. A curve near upper-left corner is more desirable, which means that the method requires less memory to achieve some certain test accuracy.

We observe significant advantage of LM-RFF over full-precision RFF in terms of memory efficiency. For example, to achieve $95\%$ accuracy on Isolet, LM-RFF (both 1-bit and 2-bit) requires $\approx 2000$ bits per sample, while RFF needs $\approx 18000$ bits, leading to a 9x compression ratio. Similar comparison holds for all datasets, and in general the compression ratio of LM-RFF is around 10x.

When compared with StocQ, we see consistently advantage of LM-RFF in memory cost. In general, LM-RFF can further improve the compression ratio of StocQ by 2x $\sim$ 4x. Additionally, LM-RFF typically requires fewer-bit quantizers (smaller $b$ ) than StocQ to achieve satisfactory accuracy.

2 Kernel Ridge Regression (KRR)

We summarize KRR results in Figure 10. Again, with same $b$ and number of RFFs, LM-RFF consistently beats StocQ especially with low bits. We see that 1-bit LM-RFF even outperforms 2-bit StocQ, and when $b=4$ , we still observe considerable advantage of LM-RFF over StocQ. In the second row, we present the memory efficiency comparison. Note that, due to high-order terms in the true model, the test MSE of linear kernel is $20.8$ , while learning with full-precision RFF significantly reduces it to $3.5$ . With largest memory budget that is tested, 1-bit and 2-bit LM-RFF yield $5.9$ and $4.1$ test MSE respectively, which are already quite close to 3.5, while for 1-bit and 2-bit StocQ, the test losses are $14.5$ and $5.0$ respectively, much worse than those of LM-RFF. We again see significant storage saving of LM-RFF. For instance, to reach the same test MSE (e.g., 10), the compression ratio is about 5x for $b=4$ compared with full-precision RFF. Moreover, the advantage of LM-RFF over StocQ is also significant for this regression problem.

3 Scale-invariant Kernel Approximation Error

Recall the notation $U=[u_{1},...,u_{n}]^{T}$ as the data matrix. Let $\mathcal{K}$ be the $n\times n$ Gaussian kernel matrix, with $\mathcal{K}_{ij}=K(u_{i},u_{j})$ . Denote $\hat{\mathcal{K}}$ as the estimated kernel matrix by an approximation algorithm. Kernel Approximation Error (KAE) has been shown to play an important role in the generalization of learning with random features, including the norms (Cortes et al., 2010; Gittens and Mahoney, 2013; Sutherland and Schneider, 2015) of $\hat{\mathcal{K}}-\mathcal{K}$ and spectral approximations (Bach, 2013; Alaoui and Mahoney, 2015; Avron et al., 2017; Zhang et al., 2019). We investigate the KAEs to better justify the impressive generalization ability of LM-RFF from a theoretical aspect.

Let $\mathcal{K}$ be a kernel matrix and $\hat{\mathcal{K}}$ be its randomized approximation. We define

Denote the minimizers as $\beta_{2}^{*}$ and $\beta_{F}^{*}$ , respectively. Define

Our new KAE metrics are more general, adapted to the best scaling factor $\beta_{2}^{*}$ or $\beta_{F}^{*}$ of the estimated kernel. Since LM-RFF estimators are slightly biased (recall Observations 4.1 and 4.2), Definition 5.1 is important for appropriately evaluating our proposed LM-RFF kernel estimation approach. In Figure 11, we provide scale-invariant $\|\cdot\|_{2}^{*}$ , $\|\cdot\|_{F}^{*}$ and $\delta_{2}^{*}$ metricsZhang et al. (2019) found that for kernel approximation methods, $\delta_{2}$ is fairly predictive of the generalisation performance. on Isolet and BASEHOCK dataset as representatives. As we can see, LM-RFF always has smaller KAEs than StocQ with equal bits. In particular, with extreme 1-bit compression, StocQ has exceedingly large loss due to its large variance, while in many cases the KAEs of 1-bit LM-RFF are already quite small. The KAE comparison well aligns with, and to an extent explains, our experimental results in Section 5.1 and Section 5.2 that 1) LM-RFF consistently outperforms StocQ, and 2) 1-bit StocQ generalizes very poorly. Thus, it provides a general justification of the superior effectiveness of LM-RFF in machine learning.

Conclusion

The technique of random Fourier features (RFF) is a popular method to solve the computational bottleneck in large-scale (Gaussian) kernel learning tasks. In this paper, we study quantization methods to compress RFFs for substantial memory savings and efficient computations. In particular, we focus on developing quantization algorithms based on the Lloyd-Max (LM) framework and propose two methods named LM-RFF and LM2-RFF. In addition, we also analyze a method based on stochastic rounding (StocQ). Both theoretically and empirically, LM-RFF significantly outperforms StocQ on many tasks, especially when the number of bits is not large. Compared to full-precision (e.g., 32- or 64-bit) RFFs, the experiments imply that often a 2-bit LM-RFF quantizer achieve comparable performance with full-precision, at a substantial (e.g., 10x) saving in memory cost, which would be highly beneficial in practical applications.

References

Appendix A Lloyd-Max (LM) Quantization: Derivation and Properties

We provide a detailed derivation of Lloyd-Max (LM) quantization scheme and its properties, which would be useful to our analysis. Recall that our proposed LM-RFF quantizers minimize the distortion defined as

where $f(z)$ is the signal distribution. Also, our $b$ -bit fixed quantizer $Q$ has borders $t_{0}<...<t_{M}$ and reconstruction levels $\mu_{1}<...<\mu_{M}$ , with $M=2^{b}$ . Since the sine and cosine function are bounded within $ $, we have$ t_{0}=-1 $and$ t_{M}=1$. Thus the distortion is

Lloyd’s algorithm finds a stationary point of above system. By setting the derivative of $D_{Q}$ w.r.t. $\mu_{i}$ to 0

We do the same thing for $t_{i}$ (i.e., setting $\frac{\partial D_{Q}}{\partial t_{i}}=0$ ) and get

The following two useful properties hold for LM quantizers.

Appendix B More Analytical Figures in Section 4

In Figure 12, we present more figures on the bias of LM quantized estimators, corresponding to Theorem 4.3, Theorem 4.4. Same as in the main paper, we see that the proposed surrogates (Observations 4.1 and 4.2) align well with true biases. As $b$ increases, the bias vanishes towards .

In Figure 13, we provide more plots on variance of proposed LM-RFF estimators at more $\gamma$ levels. As we expect, the variances of LM-RFF quantized estimators converge to the corresponding full-precision estimators as the number of bits $b$ increases, i.e., $Var[\hat{K}_{Q}]\rightarrow Var[\hat{K}]$ , $Var[\hat{K}_{n,Q}]\rightarrow Var[\hat{K}_{n}]$ , as $b\rightarrow\infty$ .

Appendix C Proofs

(of Lemma 2.1) We have the convolution of uniform and Gaussian distribution as

(of Theorem 2.2) Denote $Y=\gamma X+\tau$ . We have

where $f(y)$ is given by Lemma 2.1. Let the density of Z be $g_{Z}$ , and denote $t^{*}=\cos^{-1}z$ . It follows that

To prove the last line, denote the term in the bracket as $\alpha_{k}$ . By cancellation, for any $k_{1},k_{2}$ , we have

which equals to $2$ in the limit $k_{1}\rightarrow-\infty,k_{2}\rightarrow\infty$ . Using a similar approach, we can show that Eq. (11) is exactly the density of the cosine of a uniform random variable on $[0,2\pi]$ . For $Z_{2}=\cos(\gamma X+\tau)=Z^{2}$ , we have

Taking the derivative we get the p.d.f. as

C.2 Lemma 2.3 & Theorem 2.4

(of Lemma 2.3) Similar to the proof of Lemma 2.1, we have

where $\phi_{\sqrt{2(1-\rho)}\gamma}$ is the density of $N(0,2(1-\rho)\gamma^{2})$ . ∎

(of Theorem 2.4) Denote $Z_{x}=\cos(t_{x}),Z_{y}=\cos(t_{y})$ . Let $a_{x}^{*}=\cos^{-1}(z_{x}),a_{y}^{*}=\cos^{-1}(z_{y})$ . Denote $\phi=\phi_{\sqrt{2(1-\rho)}\gamma}$ for simplicity. We have

where (a) is derived by writing the summations $\sum_{k_{x}=-\infty}^{\infty}\sum_{k_{y}=-\infty}^{\infty}\{\cdot\}$ into $\sum_{l=-\infty}^{\infty}\sum_{k_{x}=-\infty}^{\infty}\{\cdot\}$ with $l=k_{x}-k_{y}$ and canceling terms, along with the symmetry of $\phi(\cdot)$ . This gives the joint density of $z_{x}$ and $z_{y}$ .

For the sine counterpart, with some abuse of notation, let us denote $z_{x}=\sin(t_{x})$ and $z_{y}=\sin(t_{y})$ from now on. Using similar argument, we have

After simplification, we finally arrive at

Considering $Z_{x}=\sin(t_{x}),Z_{y}=\sin(t_{y})$ . Since $\sin^{-1}(x)=\frac{\pi}{2}-\cos^{-1}(x)$ , we can substitute into the density to derive

which is the same as the previous cosine transformation. This completes the proof. ∎

C.3 Proposition 2.5

Let us denote $\sigma=\sqrt{2(1-\rho)}\gamma$ for simplicity. By symmetry and exchangeability of $f$ , to prove the desired result, it suffices to consider the case where both $z_{x}$ and $z_{y}$ are positive, i.e., $(z_{x},z_{y})\in(0,1]^{2}$ . Define the notation $a_{x}^{*}=\sin^{-1}(z_{x})\geq 0,a_{y}^{*}=\sin^{-1}(z_{y})\geq 0$ . From (12), we deduct

where we let $d=a_{x}^{*}-a_{y}^{*}$ and $d=a_{x}^{*}+a_{y}^{*}$ , and we use the fact that $\phi_{\sigma}(-x)=\phi_{\sigma}(x)$ . Note that, we consider $z_{y}>0$ so that $d\neq s$ , since when $z_{y}=0$ we trivially have $f(z_{x},0)=f(z_{x},0)$ . For now, we assume that $z_{x}\geq z_{y}>0$ , such that $d$ and $s$ are defined on the domain $0<s\leq\pi$ and $0\leq d<\min\{s,\pi-s\}$ . Since

we know that $\phi_{\sigma}$ is piecewise concave on $(0,\sigma)$ and piecewise convex on $(\sigma,\infty)$ . Thus,

for any $\sigma\leq a\leq c$ and $g\geq 0$ . The equality holds only when $a=c$ or $g=0$ . Consequently, under the assumption that $\sigma\leq\pi$ , $M_{k}\geq 0$ for $k\geq 2$ since $2\pi-s\geq\sigma$ , where the equality holds only when $d=s$ , i.e., $z_{y}=0$ . Furthermore, the piecewise convexity of $\phi_{\sigma}(\cdot)$ and (14) imply that for $\sigma\leq a<c$ ,

Also note that the function $e^{-x}$ is convex on the real line, which gives for $\forall a<c$ ,

Now that $M_{k}>0$ for $k\geq 2$ , evaluating (13) we obtain

where (a) uses (15) and (16), and (b) is because $s\leq\pi-d$ . It is easy to verify that the ratio

for $\sigma\leq\pi$ and $0\leq d<\min\{s,\pi-s\}<\frac{\pi}{2}$ . Therefore, we have proved that $f(z_{x},z_{y})>f(z_{x},-z_{y})$ , for $z_{x}\geq z_{y}>0$ . Now, by exchangeability and symmetry of $f$ , we have

Therefore, our result also holds for $z_{y}\geq z_{x}>0$ . The proof is now complete.

C.4 Theorem 4.1

Denote the StocQ quantizer as $Q$ . For each RFF $z$ , assume $z\in[t_{i-1},t_{i}]$ for some $i$ . We can then write $Q(z)=z+\epsilon$ , where

For two data vectors $u,v$ , let $F^{StocQ}(u)=\sqrt{2}Q(z_{u})$ and $F^{StocQ}(v)=\sqrt{2}Q(z_{v})$ , where $z_{u}=\cos(w^{T}u+\tau)$ and $z_{v}=\cos(w^{T}v+\tau)$ follows the distribution $f$ given by Theorem 2.4. We can write $Q(z_{u})=z_{u}+\epsilon_{u}$ , $Q(z_{v})=z_{v}+\epsilon_{v}$ where $\epsilon_{u}$ and $\epsilon_{v}$ are independent. Let $\hat{=}F^{StocQ}(u)F^{StocQ}(v)$ . We have

implying that StocQ estimate is unbiased. The variance factor can be computed as

where $Var[\hat{K}]$ is the variance of full-precision RFF kernel estimator. Obviously, $A>0$ , thus StocQ estimator always has larger variance than full-precision RFF. Continuing our analysis,

where equation $(a)$ is due to the symmetry of density $f$ and the borders $t_{0}<...<t_{2^{b}-1}$ . Substituting above expressions into (17) and cancelling terms, we obtain

The proof is completed by noting that StocQ estimator is the average of i.i.d. Bernoulli random variables. ∎

C.5 Theorem 4.3

For simplicity, we prove the result specifically for LM-RFF quantization. Similar arguments holds for general quantizers. The Chebyshev polynomials [Borwein and Erdélyi, 1995] of the first kind are defined through trigonometric identities

$\{T_{0},T_{1},...\}$ forms an orthogonal basis of the function space on $ $with finite number of discontinuities. Precisely, define the inner product w.r.t. measure$ \frac{1}{\sqrt{1-x^{2}}}$ as

By Chebyshev functional decomposition, our LM quantizer can be written as

where $\alpha_{k}$ are computed through the inner products,

This proves the first part. There is an intrinsic constraint on $\alpha_{i}$ , $i=3,5,...$ . First, we can compute the cosine of $Q(x)$ and each $T_{i}(x)$ as

Since the Chebyshev polynomials form an orthogonal basis of function space on $ $, it holds that$ \sum_{i=0}^{\infty}c_{i}^{2}=1 $. Therefore, we have$ \sum_{i=0}^{\infty}\alpha_{i}^{2}=1-2D $. Now that$ \alpha_{i}=0 $when$ i $is even, and$ \alpha_{1}=1-2D $, we then have$ \sum_{i=3,odd}^{\infty}\alpha_{i}^{2}=1-2D-(1-2D)^{2}=2D(1-2D)$.

When $\rho=1$ ( $K(u,v)=1$ ), we have $\int_{-1}^{1}T_{i}(z_{x})T_{j}(z_{x})f(z_{x})dz_{x}=0$ for $i\neq j$ by orthogonality of Chebyshev polynomials, where $f(z_{x})$ is the marginal distribution of $z_{x}$ . It follows that

This completes the proof of the theorem. ∎

C.6 Theorem 4.4

Furthermore, we have the expectation of $\Lambda$ as

This completes the proof for asymptotic mean. With a little abuse of notation, let $\hat{K}_{n,Q}=\frac{a}{\sqrt{bc}}$ , with

The gradient vector at the expectations is

The theorem is proved by plugging in the expressions. ∎

C.7 Theorem 4.5

Thus, we can compute the debiased estimator variance as (after simplification)

where the inequality is due to the fact that $\zeta_{4,0}-\zeta_{2,0}^{2}=Var[Q^{2}(z_{x})]\geq 0$ . Here we denote $M$ as a function of $\rho$ . At $\rho=0$ , we have

so that $M(0)=0$ . At $\rho=1$ , it holds that

hence $M(1)>0$ . Notice that $Q(\cdot)$ and $Q^{3}(\cdot)$ are non-decreasing odd functions, and $Q^{2}(\cdot)$ is a even function. For $\rho\in$ , since $\sqrt{2(1-\rho)}\gamma\leq\sqrt{2}\gamma\leq\pi$ by assumption, it follows from Theorem 4.7 that $\zeta_{1,1}$ , $\zeta_{2,2}$ and $\zeta_{3,1}$ are all increasing in $\rho$ on $ $. Consequently,$ M(\rho)>0 $for any$ \rho\in$. The desired result thus follows. ∎

C.8 Lemma 4.6

(of Lemma 4.6) We use the technique of Gaussian interpolation and Stein’s Lemma. First, we formulate $Y=\gamma\rho X+\gamma\sqrt{1-\rho^{2}}Z$ where $Z\sim N(0,1)$ independent of $X$ . By continuity and boundedness of $g_{1}$ and $g_{2}$ , it holds that

We analyze two parts respectively. By Lemma C.1 and law of total expectation, we have

To prove the monotonicity, suppose that $g_{1}$ and $g_{2}$ are increasing odd or non-constant even functions. So, $g_{1}^{\prime}(-x)g_{2}^{\prime}(-x)=g_{1}^{\prime}(x)g_{2}^{\prime}(x)>0$ , $\forall x\in$ . Assume $\sqrt{2(1-\rho)}\gamma\leq\pi$ , and denote $f(x,y)$ as the joint density given by Theorem 2. We can write

where (a) is due to the symmetry of $f$ and $g$ , and (b) is a consequence of Proposition 2.5 that $f(z_{x},z_{y})>f(z_{x},-z_{y})$ for all $z_{x},z_{y}\in(0,1]^{2}$ , provided that $\sqrt{2(1-\rho)}\gamma\leq\pi$ . The proof is complete.

C.9 Theorem 4.7

Since $Q_{1}$ and $Q_{2}$ both are non-decreasing and have finite number of discontinuities, by Baire’s Characterization Theorem we know that each of them is the pointwise limit of a sequence of continuous increasing functions. Suppose that $\{g_{1,n}\}$ and $\{g_{2,n}\}$ are two sequences of continuous increasing functions such that as $n\rightarrow\infty$ , $g_{1,n}\rightarrow Q_{1}$ and $g_{2,n}\rightarrow Q_{2}$ with pointwise convergence. By dominated convergence theorem, we have

where Lemma 4.6 is adopted for continuous $g_{1,n}$ and $g_{2,n}$ functions. ∎