Random Features for Kernel Approximation: A Survey on Algorithms, Theory, and Beyond

Introduction

RFF spawns a new direction for kernel approximation, and the past ten years has witnessed a flurry of research papers devoted to this topic. On the algorithmic side, subsequent work has focused on improving the kernel approximation quality and decreasing the time and space complexities . Implementation of RFF has in fact been taken to the hardware level . On the theoretical side, a series of works aim to address the following two key questions:

Approximation: how many random features are needed to ensure high quality of kernel approximation?

Generalization: how many random features are needed to incur no loss in the expected risk of a learned estimator?

Here “no loss” means how large $s$ should be for the (approximated) kernel estimator with $s$ random features to be almost as good as the exact one. Much research effort has been devoted to this direction, including analyzing the kernel approximation error (the first question above) , and studying the risk and generalization properties (the second question above) . Increasingly refined and general results have been obtained over the years. In the Thirty-sixth International Conference on Machine Learning (ICML 2019), Li et al. were recognized by the Honorable Mentions (best paper finalist) for their unified theoretical analysis of RFF.

RFF has proved effective in a broad range of machine learning tasks. Given its remarkable empirical success and the rapid growth of the related literature, we believe it is desirable to have a comprehensive overview on this topic summarizing the progress in algorithm design and applications, and elucidating existing theoretical results and their underlying assumptions. With this goal in mind, in this survey we systematically review the work from the past ten years on the algorithms, theory and applications of random features methods. Figure 1 shows a schematic overview of the history of the work on random features in recent years. The main contributions of this survey include:

We provide an overview of a wide range of random features based algorithms, re-organize the formulation of representative approaches under a unifying framework for a direct understanding and comparison.

We summarize existing theoretical results on the kernel approximation error measured in various metrics, as well as results on generalization risk of kernel estimators. The underlying assumptions in these results are discussed in detail. In particular, we (partly) answer an open question in this topic: why good kernel approximation performance cannot lead to good generalization performance?

We systematically evaluate and compare the empirical performance of representative random features based algorithms under different experimental settings.

We discuss recent research trends on (high dimensional) random features in over-parameterized settings for understanding generalization properties of over-parameterized neural networks as well as the gaps in existing theoretical analysis. We view this topic as a promising research direction.

The remainder of this paper is organized as follows. Section 2 presents the preliminaries and a taxonomy of random features based algorithms. We review data-independent algorithms in Section 3 and data-dependent approaches in Section 4. In Section 5, we survey existing theoretical results on kernel approximation and generalization performance. Experimental comparisons of representative random features based methods are given in Section 6. In Section 7, we discuss recent results on random features in over-parameterized regimes. The paper is concluded in Section 8 with a discussion on future directions.

Preliminaries and Taxonomies

In this section, we introduce preliminaries on the problem setting and theoretical foundation of random features. We then present a taxonomy of existing random features based algorithms, which sets the stage for the subsequent discussion. A set of commonly used parameters is summarized in Table I.

By the representer theorem , the above problem can be rewritten as a finite-dimensional empirical risk minimization problem

2 Theoretical Foundation of Random Features

The theoretical foundation of RFF builds on Bochner’s celebrated characterization of positive definite functions.

where $\mu_{k}$ is a positive finite measure on the frequencies $\bm{\omega}$.

According to Bochner’s theorem, the spectral distribution $\mu_{k}$ of a stationary kernel $k$ is the finite measure induced by a Fourier transform. By setting $k(0)=1$, we may normalize $\mu_{k}$ to a probability density $p$ (the Fourier transform associated with $k$), hence

where the symbol $\bm{z}^{*}$ denotes the complex conjugate of $\bm{z}$. The kernels used in practice are typically real-valued and thus the imaginary part in Eq. (4) can be discarded. According to Eq. (4), RFF makes use of the standard Monte Carlo sampling scheme to approximate $k(\bm{x},\bm{x}^{\prime})$. In particular, one uses the approximation

where $\Lambda_{i}\geq 0$ are the Fourier coefficients, $Y_{i,j}$ is the spherical harmonics, and N(d,i)=\frac{2i+d-2}{i}\left(\begin{array}[]{c}i+d-3\\ d-2\end{array}\right).

Further, for a fully-connected deep neural network (more than two layers) and fixed random weights before the output layer, if the hidden layers are wide enough, one can still approach a kernel obtained by letting the widths tend to infinity . If both intermediate layers and the output layer are trained by (stochastic) gradient descent, for the network $f(\bm{x};\bm{\theta})$ with large enough $s$, the model remains close to its linearization around its random initialization throughout training, known as lazy training regime . Learning is then equivalent to a kernel method with another architecture-specific kernel, known as neural tangent kernel (NTK, ). Interestingly, NTK for two-layer ReLU networks can be constructed by arc-cosine kernels, i.e., $k\left(\bm{x},\bm{x}^{\prime}\right)=\|\bm{x}\|\|\bm{x}^{\prime}\|[u\kappa_{0}(u)+\kappa_{1}(u)]$. In fact, there is an interesting line of work showing insightful connections between kernel methods and (over-parameterized) neural networks, but this is out of scope of this survey on random features. We suggest the readers refer to some recent literature for details.

Further, if we consider the general non-stationary kernels , the spectral representation can be generalized by introducing two random variables $\bm{\omega}$ and $\bm{\omega}^{\prime}$.

() A non-stationary kernel $k$ is positive definite if and only if it admits

3 Commonly used kernels in Random Features

Random features based algorithms often consider the following kernels:

i) Gaussian kernel: Arguably the most important member of shift-invariant kernels, the Gaussian kernel is given by

where $\varsigma>0$ is the kernel width. The density (see Theorem 1 or Eq. (6)) associated with the Gaussian kernel is Gaussian $\bm{\omega}\sim\mathcal{N}(0,\varsigma^{-2}\bm{I}_{d})$.

ii) arc-cosine kernels: This class admits Eq. (6) by sampling from the Gaussian distribution $\mathcal{N}(0,\bm{I}_{d})$, that can be connected to a two-layer neural networks with various activation functions. Following , we define the $b$-order arc-cosine kernel by

where $\theta=\cos^{-1}\left(\frac{\bm{x}^{\!\top}\bm{x}^{\prime}}{\|\bm{x}\|_{2}\|\bm{x}^{\prime}\|_{2}}\right)$ and

Most common in practice are the zeroth order ($b=0$) and first order ($b=1$) arc-cosine kernels. The zeroth order kernel is given explicitly by

iii) Polynomial kernel: This is a widely used family of non-stationary kernels given by

where $b$ is the order of the polynomial.

4 Taxonomy of random features based algorithms

The key step in random features based algorithms is constructing the following random feature mapping

Data-independent random features based algorithms can be further categorized into three classes according to their sampling strategy:

iii) Quadrature based methods: Numerical integration techniques can be also used to approximate the integral representation in Eq. (4). These techniques may involve deterministic selection of the points and weights, e.g., by using Gaussian Quadrature (GQ) or Sparse Grids Quadrature (SGQ) over each dimension (their integration formulation can be found in the first blue box in Figure 2). The selection can also be randomized. For example, in the work , the $d$-dimensional integration in Eq. (4) is transformed to a double integral, and then approximated by using the Stochastic Spherical-Radial (SSR) rule (see the second blue box in Figure 2).

Data-dependent algorithms use the training data to guide the selection of points and weights in the random features for better approximation quality and/or generalization performance. These algorithms can be grouped into three classes according to how the random features are generated.

i) Leverage score sampling: Built upon the importance sampling framework, this class of algorithm replaces the original distribution $p(\bm{\omega})$ by a carefully chosen distribution $q(\bm{\omega})$ constructed using leverage scores (see the yellow box in Figure 2). The representative approach in this class is Leverage Score based RFF (LS-RFF) , and its accelerated version .

ii) Re-weighted random feature selection: Here the basic idea is to re-weight the random features by solving a constrained optimization problem. Examples of this approach include weighted RFF , weighted QMC , and weighted GQ . Note that these algorithms directly learn the weights of pre-given random features. Another line of methods re-weight the random features using a two-step procedure: i) “up-projection”: first generate a large set of random features $\{\bm{\omega}_{i}\}_{i=1}^{J}$ ; ii) “compression”: then reduce these features to a small number (e.g., $10^{2}\sim 10^{3}$) in a data-dependent manner, e.g., by using kernel alignment , kernel polarization , or compressed low-rank approximation .

iii) Kernel learning by random features: This class of methods aim to learn the spectral distribution of kernel from the data so as to achieve better similarity representation and prediction. Note that these methods learn both the weights and the distribution of the features, and hence differ from the other random features selection methods mentioned above, which assume that the candidate features are generated from a pre-given distribution and only learn the weights of these features. Representative approaches for kernel learning involve a one-stage or two-stage procedure . From a more general point of view, the aforementioned re-weighted random features selection methods can also be classified into this class. Since these methods belong to the broad area of kernel learning instead of kernel approximation, we do not detail them in this survey.

Besides the above three main categories, other data-dependent approaches include the following. i) Quantization random features : Given a memory budget, this method quantizes RFF for Gaussian kernel approximation. A key observation from this work is that random features achieve better generalization performance than Nyström approximation under the same memory space. ii) Doubly stochastic random features : This method uses two sources of stochasticity, one from sampling data points by stochastic gradient descent (SGD), and the other from using RFF to approximate the kernel. This scheme has been used for Kernel PCA approximation , and can be further extended to triply stochastic scheme for multiple kernel approximation .

Data-independent Algorithms

We describe several representative data-independent algorithms based on Monte Carlo sampling, using the Gaussian kernel $k(\bm{x},\bm{x}^{\prime})=k(\bm{\tau})=\exp(-{\|\bm{\tau}\|^{2}_{2}}/{2\varsigma^{2}})$ as an example. Note that these algorithms often apply to more general classes of kernels, as summarized in Table II.

RFF : For Gaussian kernels, RFF directly samples the random features from a Gaussian distribution (corresponds to the inverse Fourier transform): $\{\bm{\omega}\}_{i=1}^{s}\sim p(\bm{\omega})$. In particular, the corresponding transformation matrix is given by

Fastfood : By observing the similarity between the dense Gaussian matrix and Hadamard matrices with diagonal Gaussian matrices, Le et al. firstly introduce Hadamard and diagonal matrices to speed up the construction of dense Gaussian matrices in RFF, especially in high dimensions (e.g., $d\geq 1000$). In particular, $\bm{W}$ used in Eq. (8) is substituted by

where $\bm{H}$ is the Walsh-Hadamard matrix admitting fast multiplication in $\mathcal{O}(d\log d)$ time, and $\bm{\Gamma}\in\{0,1\}^{d\times d}$ is a permutation matrix that decorrelates the eigen-systems of two Hadamard matrices. The three diagonal random matrices $\bm{G}$, $\bm{B}_{1}$ and $\bm{B}_{2}$ are specified as follows: $\bm{G}$ has independent Gaussian entries drawn from $\mathcal{N}(0,1)$; $\bm{B}_{1}$ is a random scaling matrix with $(\bm{B}_{1})_{ii}=\|\bm{\omega}_{i}\|_{2}/\|\bm{G}\|_{\text{F}}$, which encodes the spectral properties of the associated kernel; $\bm{B}_{2}$ is a binary decorrelation matrix with independent random $\{\pm 1\}$ entries. FastFood is an unbiased estimator, but may have a larger variance than RFF:

which converges at an $\mathcal{O}(1/s)$ rate.

$\mathcal{P}$-model : A general version of Fastfood, the $\mathcal{P}$-model constructs the transformation matrix as

SCRF : It accelerates the construction of random features by using circulant matrices. The transformation matrix is

The above three approaches are designed to accelerate the computation of RFF. We next overview representative methods that aim for better approximation performance than RFF.

which is related to the normalized linear kernel . Albeit simple, NRFF is effective in variance reduction and in particular satisfies

ORF : It imposes orthogonality on random features for the Gaussian kernel and has the transformation matrix

where $\bm{Q}$ is a uniformly distributed random orthogonal matrix, and $\bm{S}$ is a diagonal matrix with diagonal entries sampled i.i.d from the $\chi$-distribution with $d$ degrees of freedom. This orthogonality constraint is useful in reducing the approximation error in random features. It is also considered in for unifying orthogonal Monte Carlo methods. ORF is unbiased and with variance bounded by

where we have $g(\tau)={e^{\tau^{2}}\left(\tau^{8}+6\tau^{6}+7\tau^{4}+\tau^{2}\right)}/{4}$ $+{e^{\tau^{2}}\tau^{4}\left(\tau^{6}+2\tau^{4}\right)}/{2d}$. It can be seen that the variance reduction property $\text{Var}[\text{ORF}]<\text{Var}[\text{RFF}]$ holds under some conditions, e.g., when $d$ is large and $\tau$ is small. For a large $d$, the ratio of the variances of ORF and RFF can be approximated by

Choromanski et al. further improve the variance bound to

SORF : It replaces the random orthogonal matrices used in ORF by a class of structured matrices akin to those in Fastfood. The transformation matrix of SORF is given by

where $\bm{H}$ can be the normalized Hadamard matrix or the Walsh matrix, and $\bm{D}_{i}$ is the Rademacher matrix as defined in SORF. Theoretical results in show that the ROM estimator achieves variance reduction compared to RFF. Interestingly, odd values of $t$ yield better results than even $t$. This provides an explanation for why SORF chooses $t=3$.

LP-RFF : It attempts to quantize RFF with the Gaussian kernel under a memory budget, i.e., mapping each $s$-dimensional random feature $z_{p}(\bm{x})=\sqrt{2/s}\cos(\bm{W}_{\text{RFF}}\bm{x})\in[-\sqrt{2/s},\sqrt{2/s}]$ to an $s$-dimensional low precision vector with $b$ bits via a stochastic rounding scheme. They divide the interval $[-\sqrt{2/s},\sqrt{2/s}]$ into $2^{b}-1$ equal-sized sub-intervals and randomly round each value $\sqrt{2/s}\cos(\bm{\omega}_{i}\bm{x})$ to either the top or bottom of the corresponding sub-interval. Strictly speaking, this method does not belong to data-independent algorithms. But we put it here for ease of description as this approach directly quantizes RFF. More importantly, a new insight demonstrated by this method is that, under the same memory budget, random features based algorithms achieve better generalization performance than Nyström approximation . Apart from the stochastic quantization scheme used in , the authors of employ Lloyd-Max quantization with a smaller number of bits.

From the above description, one can find that orthogonalization is a typical operation for variance reduction, e.g., ORF/SORF/ROM. Here we take the Gaussian kernel as an example to illustrate insights of such scheme. By sampling $\{\bm{\omega}_{i}\}_{i=1}^{s}\sim\mathcal{N}(\bm{0},\varsigma^{-2}\bm{I}_{d})$, the used Gaussian distribution is isotropic and only depends on the norm $\|\bm{\omega}\|_{2}$ instead of $\bm{\omega}$. The used orthogonal operator makes the direction of $\bm{\omega}_{i}$ orthogonal to each other (that means more uniform) while retaining its norm unchangedIn fact, while orthogonalization only makes the direction of $\{\bm{\omega}_{i}\}_{i=1}^{s}$ more uniform, one can make the length $\|\bm{\omega}_{i}\|_{2}$ uniform by sampling from the cumulative distribution function of $\|\bm{\omega}\|_{2}$., which leads to decrease the randomness in Monte Carlo sampling, and thus achieve variance reduction effect. If we attempt to directly decrease the randomness in Monte Carlo sampling, QMC is a powerful way to achieve this goal and can then be used to kernel approximation. This is another line of random features with variance reduction illustrated as below.

2 Quasi-Monte Carlo Sampling

Here we briefly review methods based on quasi-Monte Carlo sampling (QMC) , spherical structured feature (SSF) , and moment matching (MM) . These three methods achieve a lower variance or approximation error than RFF. Strictly speaking, the later two algorithms do not belong to the quasi-Monte Carlo sampling framework. However, SSF and MM share the same integration formulation with QMC and thus we introduce them here for simplicity.

Classical Monte Carlo sampling generates a sequence of samples randomly and independently, which may lead to an undesired clustering effect and empty spaces between the samples . Instead of fully random samples, QMC outputs low-discrepancy sequences. A typical QMC sequence has a hierarchical structure: the initial points are sampled on a coarse scale whereas the subsequent points are sampled more finely. For approximating a high-dimensional integral, QMC achieves an asymptotic error convergence rate of $\epsilon=\mathcal{O}((\log s)^{d}/s)$, which is faster than the $\mathcal{O}(s^{-1/2})$ rate of Monte Carlo. Note however that QMC often requires $s$ to be exponential in $d$ for the improvement to manifest.

MM : It also uses the transformation matrix in Eq. (14), but generates a $d$-dimensional uniform sampling sequence $\{\bm{t}_{i}\}_{i=1}^{s}$ by a moment matching scheme instead of using a low discrepancy sequence as in QMC. In particular, the transformation matrix is

3 Quadrature based Methods

Quadrature based methods build on a long line of work on numerical quadrature for estimating integrals. In quadrature methods, the weights are often non-uniform, and the points are usually selected using deterministic rules including Gaussian quadrature (GQ) and sparse grids quadrature (SGQ) . Deterministic rules can be extended to their stochastic versions. For example, Munkhoeva et al. explore the stochastic spherical-radial (SSR) rule in kernel approximation. Below we briefly review these methods.

GQ : It assumes that the kernel function $k$ factorizes with respect to the dimensions and the corresponding distribution $p({\bm{\omega}})=p([\omega^{(1)},\omega^{(2)},\ldots,\omega^{(d)}]^{\!\top})$ in Eq. (4) is sub-Gaussian. Therefore, the $d$-dimenionsal integral in Eq. (4) can be factorized as

Since each of the factors is a one-dimensional integral, we can approximate them using a one-dimensional quadrature rule. For example, one may use Gaussian quadrature with orthogonal polynomials:

In general, the univariate Gaussian quadrature with $L$ quadrature points is exact for polynomials up to $(2L-1)$ degrees. The multivariate Gaussian quadrature is exact for all polynomials of the form $\omega_{1}^{i_{1}}\omega_{2}^{i_{2}}\cdots\omega_{d}^{i_{d}}$ with $1\leq i_{j}\leq 2L-1$; however the total number of points $s=L^{d}$ scales exponentially with the dimension $d$ and thus this method suffers from the curse of dimensionality.

SGQ : To alleviate the curse of dimensionality, SGQ uses the Smolyak rule to decrease the needed number of points. Here we consider the third-degree SGQ using the symmetric univariate quadrature points $\left\{-\hat{p}_{1},0,\hat{p}_{1}\right\}$ with weights $(\hat{a}_{1},\hat{a}_{0},\hat{a}_{1})$:

where the function $g(\bm{\omega}):=\sigma(\bm{\omega}^{\top}\bm{x})\sigma(\bm{\omega}^{\!\top}\bm{x}^{\prime})$ is given by Eq. (6), and $\bm{e}_{i}$ is the $d$-dimensional standard basis vector with the $i$-th element being 1. The corresponding transformation matrix is

which leads to the explicit feature mapping

where $\bm{w}_{i}$ is the $i$-th row of $\bm{W}_{\text{SGQ}}$. Note that SGQ generates $2d+1$ points. To obtain a dimension-adaptive feature mapping, Dao et al. propose to subsample the points according to the distribution determined by their weights such that the mapping feature dimension is equal to $s$.

SSR : It transforms Eq. (6) (actually a $d$-dimensional integral) to a double integral over a hyper-sphere and the real line. Let ${\bm{\omega}}=r\bm{u}$ with $\bm{u}^{\!\top}\bm{u}=1$ for $r\in[0,\infty)$, we have

with $\bm{\vartheta}=[\vartheta_{1},\vartheta_{2},\cdots,\vartheta_{s}]$ and $\bm{V}=[\bm{v}_{1},\bm{v}_{2},\cdots,\bm{v}_{d+1}]$, where $\vartheta\sim\chi(d+2)$ and $\{\bm{v}_{i}\}_{i=1}^{d+1}$ are the vertices of a unit regular $d$-simplex, which is randomly rotated by $\bm{Q}$. To get $s$ features, one may stack $s/(2d+3)$ independent copies of $\bm{W}$ as suggested by . Finally, the feature mapping by SSR is given by

where $a_{0}=\sqrt{1-\sum_{j=1}^{d+1}\frac{d}{\rho_{j}^{2}}}$, $a_{j}=\frac{1}{\rho_{j}}\sqrt{\frac{d}{2(d+1)}}$ for $j\in[s]$, and $w_{j}$ is the $j$-th element of the stacked $\bm{W}$.

In general, according to Eq. (6), kernel approximation by random features is actually a $d$-dimensional integration approximation problem in mathematics. Sampling methods and quadrature based rules are two typical classes of approaches for high-dimensional integration approximation. Efforts on quadrature based methods focus on developing a high-accuracy, mesh-free, efficiency rule, e.g., . Note that, if the integrand $g(\bm{\omega}):=\sigma(\bm{\omega}^{\top}\bm{x})\sigma(\bm{\omega}^{\!\top}\bm{x}^{\prime})$ in the integration representation (6) belongs to a RKHS, the above quadrature rules can be termed as kernel-based quadrature, e.g., Bayesian quadrature and leverage-score quadrature . This approach is in essence different from the previously studied quadrature rules in functional spaces, model formulation, and scope of application.

Data-dependent algorithms

Data-dependent approaches aim to design/learn the random features using the training data so as to achieve better approximation quality or generalization performance. Based on how the random features are generated, we can group these algorithms into three classes: leverage score sampling, random features selection, and kernel learning by random features.

Leverage score based approaches are built on the importance sampling framework. Here one samples $\{\bm{w}_{i}\}_{i=1}^{s}$ from a distribution $q(\bm{w})$ that needs to be designed, and then uses the following feature mapping in Eq. (5):

To design the distribution $q$, one makes use of the ridge leverage function in KRR:

where $\lambda$ is the KRR regularization parameter. Define

The quantity $d_{\bm{K}}^{\lambda}\ll n$ determines the number of independent parameters in a learning problem and hence is referred to as the number of effective degrees of freedom . With the above notation, the distribution $q$ designed in is given by

Compared to standard Monte Carlo sampling for RFF, leverage score sampling requires fewer Fourier features and enjoys nice theoretical guarantees (see the next section for details). Note that $q(\bm{\omega})$ can be also defined by the integral operator rather than the Gram matrix used above, but we do not strictly distinguish these two cases. The typical leverage score based sampling algorithm for RFF is illustrated in as below.

LS-RFF (Leverage Score-RFF) : It uses a subset of data to approximate the matrix $\bm{K}$ in Eq. (21) so as to compute $d_{\bm{K}}^{\lambda}$. LS-RFF needs $\mathcal{O}(ns^{2}+s^{3})$ time to generate refined random features, which can be used in KRR and SVM for prediction.

SLS-RFF (Surrogate Leverage Score-RFF) : To avoid inverting an $s\times s$ matrix in LS-RFF, SLS-RFF designs a simple but effective surrogate leverage function

where the additional term $n\bm{I}$ and the coefficient $1/(n^{2}\lambda)$ in Eq. (23) ensure that $L_{\lambda}$ is a surrogate function that upper bounds the function $l_{\lambda}$ in Eq. (20). One then samples random features from the surrogate distribution $Q(\bm{w})=\frac{L_{\lambda}({\bm{\omega}})}{\int L_{\lambda}({\bm{\omega}})\text{d}{\bm{\omega}}}$, which has the same time complexity $\mathcal{O}(ns^{2})$ as RFF. SLS-RFF and can be applied to KRR and Canonical Correlation Analysis .

Note that leverage scores sampling is a powerful tool used in sub-sampling algorithms for approximating large kernel matrices with theoretical guarantees, in particular in Nyström approximation. Research on this topic mainly focuses on obtaining fast leverage score approximation due to inversion of an $n$-by-$n$ kernel matrix, e.g., two-pass sampling (LS-RFF belongs to this), online setting , path-following algorithm , or developing various surrogate leverage score sampling based algorithms .

2 Re-weighted random features

Here we briefly review three re-weighted methods: KA-RFF by kernel alignment, KP-RFF by kernel polarization, and CLR-RFF by compressed low-rank approximation.

KA-RFF (Kernel Alignment-RFF) : It pre-computes a large number of random features that are generated by RFF, and then select a subset of them by solving a simple optimization problem based on kernel alignment . In particular, the optimization problem is

KP-RFF (Kernel Polarization-RFF) : It first generates a large number of random features by RFF and then selects a subset from them using an energy-based scheme

CLR-RFF (Compression Low Rank-RFF) : It first generates a large number of random features and then selects a subset from them by approximately solving the optimization problem

3 Kernel learning by random features

This class of approaches construct random features using sophisticated learning techniques, e.g., by learning the spectral distribution of kernel from the data.

Representative approaches in this class often involve a one-stage or two-stage process. The two-stage scheme is common when using random features. It first learns the random features, and then incorporates them into kernel methods for prediction. Actually, the above-mentioned leverage sampling and random features selection based algorithms employ this scheme. The algorithm proposed in is a typical method for kernel learning by random features. This method first learns a spectral distribution of a kernel via an implicit generative model, and then trains a linear model by these learned features.

One-stage algorithms aim to simultaneously learn the spectral distribution of a kernel and the prediction model by solving a single joint optimization problem or using a spectral inference scheme. For example, Yu et al. propose to jointly optimize the nonlinear feature mapping matrix $\bm{W}$ and the linear model with the hinge loss. The associated optimization problem can be solved in an alternating fashion with SGD. In , the kernel alignment approach in the Fourier domain and SVM are combined into a unified framework, which can be also solved using an alternating scheme by Langevin dynamics and projection gradient descent. Wilson and Adams construct stationary kernels as the Fourier transform of a Gaussian mixture based on Gaussian process frequency functions. This approach can be extended to learning with Fastfood , non-stationary spectral kernel generalization , and the harmonizable mixture kernel . Moreover, Oliva et al. propose a nonparametric Bayesian model, in which $p(\bm{\omega})$ is modeled as a mixture of Gaussians with a Dirichlet process prior. The parameters of the Gaussian mixture and the classifier/regressor model are inferred using MCMC.

Theoretical Analysis

In this section, we review a range of theoretical results that center around the two questions mentioned in the introduction and restated below:

Approximation: how many random features are needed to ensure a high quality estimator in kernel approximation?

Generalization: how many random features are needed to incur no loss of empirical risk and expected risk in a learning estimator?

Figure 3 provides a taxonomy of representative work on these two questions.

More specifically, Rahimi and Recht provide the earliest result on learning with RFF with Lipschitz continuous loss functions. Their results imply that $\Omega(n)$ random features are sufficient to incur no loss of learning accuracy. This result is improved in , which shows that $\Omega(\sqrt{n}\log n)$ random features or even less suffice for the Gaussian kernel. When using the data-dependent sampling $\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim q({\bm{\omega}})$, the above results are further improved in under various settings. Note that some results above do not directly apply to the squared loss in KRR, whose Lipschitz parameter is unbounded. For squared losses, Rudi et al. show that $\Omega(\sqrt{n}\log n)$ random features by RFF suffice to achieve a minimax optimal learning rate $\mathcal{O}(1/\sqrt{n})$. A more refined analysis is given in under the $p(\bm{\omega})$-sampling and $q(\bm{\omega})$-sampling settings.

Below we discuss the above theoretical work in more details.

Table III summarizes representative theoretical results on the convergence rates, the upper bound of the growing diameter, and the resulting sample complexity under different metrics. Here sample complexity means the number of random features sufficient for achieving a maximum approximation error at most $\epsilon$.

According to the above theorem by covering number, with $s:=\Omega(\epsilon^{-2}d\log({1}/{\epsilon\delta}))$ random features, one can ensure an $\epsilon$ uniform approximation error with probability greater than $1-\delta$. This result also applies to dot-product kernels by random Maclaurin feature maps (see [34, Theorem 8]). The quadrature based algorithm follows this proof framework, and achieves the same error bound with a smaller constant than RFF in Theorem 4 by an extra boundedness assumption. Instead, Fastfood on Gaussian kernels achieves $\mathcal{O}(\sqrt{\log(d/\delta)})$ times approximation error than RFF due to estimates for $\bm{\Gamma}\bm{H}\bm{B}_{2}$ in Eq. (9), which is based on concentration inequalities for Lipschitz continuous functions under the Gaussian distribution.

Different from the above results using Hoeffding’s inequality for the covering number bound in their proof, Sriperumbudur and Szabó revisit the above bound by refined technique of McDiarmid’s inequality, symmetrization and bound the expectation of Rademacher average by Dudley entropy bound.

Under the same assumption of Theorem 4, we have

where $h(d,|\mathcal{S}|,\sigma_{p})$ is an appropriately defined function of $d$, $|\mathcal{S}|$, and $\sigma_{p}$. For better comparison, the above inequality can be rewritten as

For the Gaussian kernel, the approximation guarantee can be further improved. In particular, the following theorem gives a probability bound independent of $d$.

[$\Delta$-spectral approximation ] For $0\leq\Delta<1$, a symmetric matrix $\bm{A}$ is a $\Delta$-spectral approximation of another symmetric matrix $\bm{B}$, if $(1-\Delta)\bm{B}\preceq\bm{A}\preceq(1+\Delta)\bm{B}$, where $A\preceq B$ indicates that $B-A$ is a semi-positive definite matrix.

According to this definition, $\bm{Z}\bm{Z}^{\!\top}+\lambda\bm{I}_{n}$ is $\Delta$-spectral approximation of $\bm{K}+\lambda\bm{I}_{n}$ if

The follow theorem gives the number of random features $s$ that are sufficient to guarantee $\Delta$-spectral approximation.

Let $k$ be a shift-invariant kernel and its associated probability distribution $p({\bm{\omega}})$ (i.e., the Fourier transform of $k$), $\Delta\leq 1/2$, $\delta\in(0,1)$, and $n_{\lambda}:=n/\lambda$. Assume that $\|\bm{K}\|_{2}\geq\lambda$ and $\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim p(\bm{\omega})$. If the total number of random features satisfies

Theorem 7 states that $\Omega(n_{\lambda}\log d_{\bm{K}}^{\lambda})$ random features are sufficient to guarantee $\Delta$-spectral approximation by the matrix Bernstein concentration inequality and effective degree of freedom, where $n_{\lambda}:=n/\lambda$. Under this framework, Choromanski et al. [97, Theorem 5.4] present a non-asymptotic comparison result between RFF and ORF for spectral approximation by virtue of the smallest singular value of $\bm{K}+n\lambda\bm{I}$.

For the Gaussian kernel, let $\widetilde{\Delta}$ be the smallest positive number such that $\widetilde{\bm{K}}+\lambda n\bm{I}_{n}$ is a $\widetilde{\Delta}$-spectral approximation of $\bm{K}+\lambda n\bm{I}_{n}$, where $\widetilde{\bm{K}}$ is an approximate kernel matrix obtained by RFF or ORF. Then, for any $\epsilon>0$ we have

The results in Theorem 7 can be improved if we consider data-dependent sampling, i.e., $\{{\bm{\omega}}_{i}\}_{i=1}^{s}$ are sampled from the empirical ridge leverage score distribution $q({\bm{\omega}})={l_{\lambda}({\bm{\omega}})}/{d^{\lambda}_{\bm{K}}}$ in Eq. (22) instead of the standard $p({\bm{\omega}})$.

Let $k$ be a shift-invariant kernel associated with the empirical ridge leverage score distribution $q(\bm{\omega})$ in Eq. (22), $\Delta\leq 1/2$ and $\delta\in(0,1)$. Assume that $\|\bm{K}\|_{2}\geq\lambda$ and $\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim q(\bm{\omega})$. If the total number of random features satisfies

Theorem 9 shows that if we sample using the ridge leverage function, then $\Omega(d_{\bm{K}}^{\lambda}\log d_{\bm{K}}^{\lambda})$ random features, which is less than $\Omega(n_{\lambda}\log d_{\bm{K}}^{\lambda})$, suffice for spectral approximation of $\bm{K}$.

The authors of generalize the notion of $\Delta$-spectral approximation to ($\Delta_{1},\Delta_{2})$-spectral approximation.

For $\Delta_{1}$,$\Delta_{2}\geq 0$, a symmetric matrix $\bm{A}$ is a $(\Delta_{1},\Delta_{2})$-spectral approximation of another symmetric matrix $\bm{B}$, if $(1-\Delta_{1})\bm{B}\preceq\bm{A}\preceq(1+\Delta_{2})\bm{B}$.

This definition is motivation by the argument that the quantities $\Delta_{1}$ and $\Delta_{2}$ in the upper and lower bounds may have different impact on the generalization performance. Using this definition, Zhang et al. derive the following approximation guarantees when one quantizes each random Fourier feature $\bm{\omega}_{i}$ to a low-precision $b$-bit representation, which allows more features to be stored in the same amount of space.

Let $\widetilde{\bm{K}}$ be an $s$-features $b$-bit LP-RFF approximation of a kernel matrix $\bm{K}$ and $\delta\in(0,1)$. Assume that $\|\bm{K}\|_{2}\geq\lambda\geq\delta^{2}_{b}=2/(2^{b}-1)^{2}$ and define $a:=8\operatorname{Tr}(\bm{K}+\lambda\bm{I}_{n})^{-1}(\bm{K}+\delta^{2}_{b}\bm{I}_{n})$. For $\Delta_{1}\leq 3/2$ and $\Delta_{2}\in[\frac{\delta_{b}^{2}}{\lambda},\frac{3}{2}]$, if the total number of random features satisfies

Theorem 10 shows that when the quantization noise is small relative to the regularization parameter, using low precision has minimal impact on the number of features required for the ($\Delta_{1},\Delta_{2})$-spectral approximation. In particular, as $s\rightarrow\infty$, $\Delta_{1}$ converges to zero for any precision $b$, whereas $\Delta_{2}$ converges to a value upper bounded by $\delta_{b}^{2}/\lambda$. If $\delta_{b}^{2}/\lambda\ll\Delta_{2}$, using $b$-bit precision has negligible effect on the number of features required to attain this $\Delta_{2}$ see Table III for a summary.

2 Risk and generalization property

The above results on approximation error are a means to an end. More directly related to the learning performance is understanding generalization properties of random features based algorithms. To this end, a series of work study the generalization properties of algorithms based on $p(\bm{\omega})$-sampling and $q(\bm{\omega})$-sampling. Under different assumptions, theoretical results have been obtained for loss functions with/without Lipschitz continuity and for learning tasks including KRR and SVM . Apart from supervised learning with random features, results on randomized nonlinear component analysis refer to , random features with matrix sketching , doubly stochastic gradients scheme , statistical consistency .

Before we detail these theoretical results, we summarize the standard assumptions imposed in existing work. Some assumptions are technical, and thus familiarity with statistical learning theory (see Section 2.1) would be helpful. We organize these assumptions in four categories as shown in Figure 4, including i) the existence of $f_{\rho}$ (Assumption 1) and its stronger version (Assumption 8); ii) quality of random features (Assumptions 2, 6, 7); iii) noise conditions (Assumptions 3, 9, 10); iv) eigenvalue decay (Assumptions 4, 5).

We first state three basic assumptions, which are needed in all of the (regression) results to be presented.

In regression task, we assume $f_{\rho}\in\mathcal{H}$.

Note that since we consider a potentially infinite dimensional RKHS $\mathcal{H}$, possibly universal , the existence of the target function $f_{\rho}$ is not automatic. However, if we restrict to a bounded subspace of $\mathcal{H}$, i.e., $\mathcal{H}_{R}=\{f\in\mathcal{H}:\|f\|\leq R\}$ with $R<\infty$ fixed a prior, then a minimizer of the risk $\mathcal{E}(f)$ always exists as long as $\mathcal{H}_{R}$ is not universal. If $f_{\rho}$ exists, then it must lie in a ball of some radius $R_{\rho,\mathcal{H}}$. The results in this section do not require prior knowledge of $R_{\rho,\mathcal{H}}$ and they hold for any finite radius.

This noise condition is weaker than the boundedness on $y$. It is satisfied when $y$ is bounded, sub-Gaussian, or sub-exponential. In particular, if $y\in[-\frac{b}{2},\frac{b}{2}]$ almost surely with $b>0$, then Assumption 3 is satisfied with $\varsigma=B=b$.

The above three assumptions are needed in all theoretical results for regression presented in this section, so we omit them when stating these results. We next introduce several additional assumptions, which are needed in some of the theoretical results.

Eigenvalue Decay Assumptions: The following assumption, which characterizes the “size” of the RKHS $\mathcal{H}$ of interest, is often discussed in learning theory.

A kernel matrix $\bm{K}$ admit the following three types of eigenvalue decays: 1) Geometric/exponential decay: $\lambda_{i}(\bm{K})\propto n\exp({-i^{1/c}})$, which leads to $d^{\lambda}_{\bm{K}}\lesssim\log(R_{0}/\lambda)$; 2) Polynomial decay: $\lambda_{i}(\bm{K})\propto ni^{-2a}$, which implies $d^{\lambda}_{\bm{K}}\lesssim(1/\lambda)^{1/2a}$; 3) Harmonic decay: $\lambda_{i}(\bm{K})\propto n/i$, which results in $d^{\lambda}_{\bm{K}}\lesssim(1/\lambda)$.

We give some remarks on the above assumption. For shift-invariant kernels, if the RKHS is small, the eigenvalues of the kernel matrix $\bm{K}$ often admit a fast decay. Consequently, functions in the RKHS are smooth enough that a good prediction performance can be achieved. On the other hand, if the RKHS is large and the eigenvalues decay slowly, then functions in the RKHS are not smooth, which would lead to a sub-optimal error rate for prediction. It can be linked to the integral operator characterizing the hypothesis space, defined as $\Sigma:L^{2}_{\rho_{\mathcal{X}}}\rightarrow L^{2}_{\rho_{\mathcal{X}}}$ such that

The integral operator $\Sigma$ plays a significant role in characterizing the hypothesis space. In particular, the decay rate of the spectrum of $\Sigma$ quantifies the capacity of the hypothesis space in which we search for the solution. This capacity in turn determines the number of random features required for accurate learning. Rudi and Rosasco consider the following assumption on $\Sigma$.

There exist $Q>0$ and $\gamma\in$ such that for any $\lambda>0$, we have

The effective dimension $\mathcal{N}(\lambda)$ measures the “size” of the RKHS, and is in fact the operator form of $d^{\lambda}_{\bm{K}}$ in Eq. (21). Assumption 5 holds if the eigenvalues $\lambda_{i}$ of $\Sigma$ decay as $i^{-1/\gamma}$, which corresponds to the eigenvalue decay of $\bm{K}$ in Assumption 4 with $\gamma:=1/(2a)$ . The case $\gamma=0$ is the more benign situation, whereas $\gamma=1$ is the worst case.

Quality of Random Features: Here we introduce several technical assumptions on the quality of random features. The leverage score in Eq. (20) admits the operator form

which is also called as the maximum random features dimension . By defintion we always have $\mathcal{N}(\lambda)\leq\mathcal{F}_{\infty}(\lambda)$. Roughly speaking, when the random features are “good”, it is easy to control their leverage scores in terms of the decay of the spectrum of $\Sigma$. Further, fast learning rates using fewer random features can be achieved if the features are compatible with the data distribution in the following sense.

With the above definition of $\mathcal{F}_{\infty}(\lambda)$, assume that there exist $\varrho\in$, and $F>0$ such that $\mathcal{F}_{\infty}(\lambda)\leq F\lambda^{-\varrho},\forall\lambda>0$.

It always holds that $\mathcal{F}_{\infty}(\lambda)\leq\kappa^{2}\lambda^{-1}$ when $z$ is uniformly bounded by $\kappa$. So the worst case is $\varrho=1$, which means that the random features are sampled in a problem independent way. The favorable case is $\varrho=\gamma$, which means that $\mathcal{N}(\lambda)\leq\mathcal{F}_{\infty}(\lambda)\leq\mathcal{O}(n^{-\alpha\gamma})$. In , the authors consider the following assumption.

The feature mapping $z(\bm{\omega},\bm{x})$ is called optimized if there is a small constant $\lambda_{0}$ such that for any $\lambda\leq\lambda_{0}$, $\mathcal{F}_{\infty}(\lambda)\leq\mathcal{N}(\lambda)=\sum_{i=1}^{\infty}\frac{\lambda_{i}(\Sigma)}{\lambda_{i}(\Sigma)+\lambda}$.

Under the previous definitions, Assumption 7 holds only when $\mathcal{F}_{\infty}(\lambda)=\mathcal{N}(\lambda)$. This assumption is stronger than the compatibility condition in Assumption 6. Note that Assumption 7 is satisfied when sampling from $q(\bm{\omega})$.

Source condition on $f_{\rho}$: The following assumption states that $f_{\rho}$ has some desirable regularity properties.

There exist $1/2\leq r\leq 1$ and $g\in{L_{\rho_{\mathcal{X}}}^{2}}$ such that $f_{\rho}(\bm{x})=(\Sigma^{r}g)(\bm{x})$ almost surely.

Since $\Sigma$ is a compact positive operator on $L_{\rho_{\mathcal{X}}}^{2}$, its $r$-th power $\Sigma^{r}$ is well defined for any $r>0$.A more general condition ($r>0$) is often considered in approximation theory; see . Assumption 8 imposes a form of regularity/sparsity of $f_{\rho}$, which requires the expansion of $f_{\rho}$ on the basis given by the integral operator $\Sigma$. Note that this assumption is more stringent than the existence of $f_{\rho}$ in $\mathcal{H}$. The latter is equivalent to Assumption 8 with $r=\frac{1}{2}$ (the worst case), in which case $f_{\rho}\in\mathcal{H}$ need not have much regularity/sparsity.

Noise Condition: The following two assumptions on noise are considered in random features for classification.

The points in $\mathcal{X}$ can be collected into two sets according to their labels as follows

For $i\in\{\pm 1\}$, the distance of a point $\bm{x}\in X_{i}$ to the set $X_{-i}$ is denoted by $\Delta(\bm{x})$. We say that the data distribution satisfies a separation condition if there exists $\Delta>0$ such that $\rho_{X}(\Delta(\bm{x})<c)=0$.

for some $q>0$. The separation condition in Assumption 10 is an extreme case of the Tsybakov’s noise assumption with $q=\infty$. It is clear that noise-free distributions satisfy this separation assumption, since the conditional probability $\eta$ is bounded away from $1/2$.

2.2 Squared loss in KRR

In this section, we review theoretical results on the generalization properties of KRR with squared loss and random features, for both the $p(\bm{\omega})$-sampling (data-independent) and $q(\bm{\omega})$-sampling (data-dependent) settings. Table IV summarizes these results for the excess risk in terms of the key assumptions imposed, the learning rates, and the required number of random features.

We begin with the remarkable result by Rudi and Rosasco . They are among the first to show that under some mild assumptions and appropriately chosen parameters, $\Omega(\sqrt{n}\log n)$ random features suffice for KRR to achieve minimax optimal rates.

Suppose that Assumption 8 (source condition) holds with $r\in[\frac{1}{2},1]$, Assumption 6 (compatibility) holds with $\varrho\in$, and Assumption 5 (capacity) holds with $\gamma\in$. Assume that $n\geq n_{0}$ and choose $\lambda:=n^{\frac{1}{2r+\gamma}}$. If the number of random features satisfies

where $c_{0}$, $c_{1}$ are constants independent of $(n$, $\lambda$, $\delta)$, and $n_{0}$ does not depends on $n$, $\lambda$, $f_{\rho}$, or $\rho$.

Theorem 11 unifies several results in that impose different assumptions. The simplest result is Theorem 1 in , which only requires the three basic Assumptions 1–3 on existence, boundedness and continuity, corresponding to the the worst case of Theorem 11 with $\varrho=\gamma=1$ and $r=1/2$. In this case, by choosing $\lambda=n^{-1/2}$, we require $\Omega(\sqrt{n}\log n)$ random features to achieve the minimax convergence rate $\mathcal{O}(n^{-1/2})$; also see Table IV.

A more refined result is given in Theorem 2 in , which accounts for the capacity of the RKHS and the regularity of $f_{\rho}$, as quantified by the parameters $\gamma\in$ (Assumption 5) and $r\in[\frac{1}{2},1]$ (Assumption 8), respectively. Under these conditions and choosing $\lambda:=n^{-\frac{1}{2r+\gamma}}$, we require \Omega\big{(}n^{\frac{1+\gamma(2r-1)}{2r+\gamma}}\log n\big{)} random features to achieve the convergence rate \mathcal{O}\big{(}n^{-\frac{2r}{2r+\gamma}}\big{)}. Note that $\gamma=1$ is the worst case, where the eigenvalues of $\bm{K}$ have the slowest decay, and $\gamma=1/(2a)\in(0,1)$ means that the eigenvalues follow a polynomial decay $\lambda_{i}\propto ni^{-2a}$. Table IV presents this result with $\gamma:=1/(2a)$ for better comparison with the other results.

The above two results apply to the standard RFF setting with data-independent sampling. When $\{\bm{\omega}_{i}\}_{i=1}^{s}$ are sampled from a data-dependent distribution satisfying the compatibility condition in Assumption 6 with $\varrho\in$, then Theorem 3 in provide an improved result. In this case, by choosing $\lambda:=n^{-\frac{1}{2r+\gamma}}$, we require \Omega\big{(}n^{\frac{\varrho+(1+\gamma-\varrho)(2r-1)}{2r+\gamma}}\log n\big{)} random features to achieve the convergence rate \mathcal{O}\big{(}n^{-\frac{2r}{2r+\gamma}}\big{)}.

If the compatibility condition is replaced by the stronger Assumption 7 (optimized distribution), satisfied by $q(\bm{\omega})$-sampling, the work derives an improved bound that is the sharpest to date. Below we state a general result from that covers both $p(\bm{\omega})$- and $q(\bm{\omega})$-sampling.

Suppose that the regularization parameter $\lambda$ satisfies $0\leq n\lambda\leq\lambda_{1}$. We consider two sampling schemes.

$\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim p(\bm{\omega})$: if $s\geq(5z_{0}^{2}/\lambda)\log({16d_{\bm{K}}^{\lambda}}/{\delta})$ and $|z(\bm{\omega},\bm{x})|\leq z_{0}$,

$\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim q(\bm{\omega})$: if $s\geq 5d^{\lambda}_{\bm{K}}\log\left(16d_{\bm{K}}^{\lambda}/\delta\right)$,

where we recall that \mathcal{E}\big{(}f_{\bm{z},\lambda}\big{)}-\mathcal{E}\left(f_{\rho}\right) is the excess risk of standard KRR with an exact kernel (see Section 2).

Remark: A sharper convergence rate can be achieved if the Rademacher complexity used in is substituted by the local Rademacher complexity , see for details.

For $p(\bm{\omega})$-sampling, Theorem 12 improves on the results of under the exponential and polynomial decays. Specifically, if $\{{\bm{\omega}}_{i}\}_{i=1}^{s}\sim p(\bm{\omega})$, Theorem 12 requires $s\propto 1/\lambda\log d^{\lambda}_{\bm{K}}$. Specialized to the exponential decay case, this result requires $\Omega(\sqrt{n}\log\log n)$ random features to achieve an $\mathcal{O}(n^{-1/2})$ learning rate, which is an improvement compared to with $\Omega(\sqrt{n}\log n)$ random features.

For $q(\bm{\omega})$-sampling, Theorem 12 shows that if $\lambda=n^{-1/2}$, then $s\propto d^{\lambda}_{\bm{K}}\log d^{\lambda}_{\bm{K}}$ random features is sufficient to incurs no loss in the expected risk if KRR, with a minimax learning rate $\mathcal{O}(n^{-1/2})$. Corollaries of this result under three different regimes of eigenvalue decay are summarized in Table IV.

Carratino et al. extend the result of to the setting where KRR is solved by stochastic gradient descent (SGD). They show that under the basic Assumptions 1–3 and some mild conditions for SGD, $\Omega(\sqrt{n})$ random features suffice to achieve the minimax learning rate $\mathcal{O}(n^{-1/2})$. This result matches those for standard KRR with an exact kernel . The above results can be improved if in addition the source condition in Assumption 8 holds, in which case $\Omega(n^{\frac{1+\alpha(2r-1)}{2r+\alpha}})$ random features suffice to achieve an $\mathcal{O}(n^{-\frac{2r}{2r+\alpha}})$ learning rate.

The work in shows that if the randomized feature map is bounded (which is weaker than Assumption 2), then we have the following out-of-sample bound

If we choose $\lambda:=n^{-1/2}$, then $\Omega(n)$ random features are sufficient to ensure an $\mathcal{O}(n^{-1/2})$ rate in the out-of-sample bound.

In this section, we empirically evaluate the kernel approximation and classification performance of representative random features algorithms on several benchmark datasets. All experiments are implemented in MATLAB and carried out on a PC with Intel® i7-8700K CPU (3.70 GHz) and 64 GB RAM. The source code of our implementation can be found in http://www.lfhsgre.org.

Kernel: We choose the popular Gaussian kernel, zero/first-order arc-cosine kernels, and polynomial kernels for experiments.

where the kernel width parameter $\varsigma$ is tuned via 5-fold inner cross validation over a grid of $\{0.01,0.1,1,10,100\}$.

To evaluate the Gaussian kernel, we conduct the following representative algorithms for comparison: RFF , ORF , SORF , ROM , Fastfood , QMC , SSF , GQ , and LS-RFF . These algorithms include both data-independent and data-dependent approaches and involve a variety of techniques including Monte Carlo and quasi-Monte Carlo sampling, quadrature rules, variance reduction, and computational speedup using structural/circulant matrices.

ii) arc-cosine kernels: Different from Gaussian kernels and polynomial kernels, the designed arc-cosine kernels can be closely connected to neural networks, which include feature spaces that mimic the sparse, nonnegative, distributed representations of single-layer threshold networks. The used zeroth order kernel is given explicitly by

which corresponds to the Heaviside step function $\sigma({\bm{\omega}}^{\top}\bm{x})=\frac{1}{2}(1+\operatorname*{sign}({\bm{\omega}}^{\top}\bm{x}))$ in Eq. (6). The first order kernel is

which corresponds to the ReLU activation function $\sigma({\bm{\omega}}^{\top}\bm{x})=\max\{0,{\bm{\omega}}^{\!\top}\!\bm{x}\}$ in Eq. (6).

Here we consider the zero/first-order arc-cosine kernel and compare these ten algorithms (used for Gaussian kernel approximation) as well. Note that, the theoretical foundation behind random features, Bocher’s theorem, is invalid to arc-cosine kernels. Thankfully, according to the formulation of arc-cosine kernels admitting in Eq. (6), the Monte Carlo sampling (e.g., RFF) is able to used for arc-cosine kernel approximation. In this case, the remaining algorithms, e.g., ORF, QMC, and Fastfood, on various sampling strategies, can be still applicable to arc-cosine kernels, at least in the algorithmic aspect.

iii) Polynomial kernel: This is a widely used family of dot product kernels given by

where $b$ is the order. In our experiments, the order is set to $b=2$. Note that, different from Gaussian kernels and arc-cosine kernels, polynomial kernels admit neither the Bochner’s theorem nor the sampling formulation in Eq. (6), so classical random features based algorithms are applicable to arc-cosine kernels but still invalid to polynomial kernels even though both of them are dot-product. As a result, algorithms for polynomial kernel approximation are often totally different. In this survey, we include three representative approaches for evaluation, including Random Maclaurin (RM) , Tensor Sketch (TS) , and Tensorized Random Projection (TRP) .

Datasets: We consider eight non-image benchmark datasets, two representative image datasets, and a ultra-large scale dataset for evaluation. Table VI gives an overview of these datasets including the number of feature dimension, training samples, test data, training/test split, and the normalization scheme. These eight non-image benchmark datasets can be downloaded from https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ or the UCI Machine Learning Repositoryhttps://archive.ics.uci.edu/ml/datasets.html.. Some datasets include a training/test partition, denoted as “no” in the random split column. For the other datasets, we randomly pick half of the data for training and the rest for testing, denoted as “yes” in the random split column. There are two typical normalization schemes used in these datasets: “mapstd” and “minmax”. The “mapstd” scheme sets each sample’s mean to 0 and deviation to 1, while the “minmax” scheme is a standard min-max scaling operation mapping the samples to the bounded set $^{d}$. Two representative image datasets are the MNIST handwritten digits dataset and the CIFAR10 natural image classification dataset , summarized in the last two rows in Table VI. The MNIST dataset contains 60,000 training samples and 10,000 test samples, each of which is a $28\times 28$ gray-scale image of a handwritten digit from 0 to 9. Here the “minmax” normalization scheme means that each pixel value is divided by 255. The CIFAR10 dataset consists of 60,000 color images of size $32\times 32\times 3$ in 10 categories, with 50,000 for training and 10,000 for test. Besides, apart from medium/large scale datasets in our experiments, we also evaluate the compared approaches on a ultra-large scale dataset MNIST 8M , which is derived from the MNIST dataset by random deformations and translations. It shares the same number of feature dimension and test data with the MNIST dataset, but has 8,100,000 training data.

2 Results for the Gaussian Kernel

Here we test various random features based algorithms, including RFF , ORF , SORF , ROM , Fastfood , QMC , SSF , GQ , LS-RFF for kernel approximation and then combine these algorithms with lr/liblinear for classification on eight non-image benchmark datasets, refer to Appendix B.1 for details. Here we summarize the best performing algorithm on each dataset in terms of the approximation quality and classification accuracy in Table VII, where we distinguish the small $s$ case (i.e., $s=2d$ or $s=4d$) and the large $s$ case (i.e., $s=16d$ or $s=32d$). The notation “-” therein means that there is no statistically significant difference in the performance of most algorithms.

In terms of approximation error, we find that SSF, ORF, and QMC achieve promising approximation performance in most cases. Recall that the goal of using random features is to find a finite-dimensional (embedding) Hilbert space to approximate the original infinite-dimensional RKHS so as to preserve the inner product. To achieve this goal, SSF, QMC, and ORF are based on a similar principle, namely, generating random features that are as independent/complete as possible to reduce the randomness in sampling. Regarding to SSF, we find that SSF performs well under the small $s$ case, but the significant improvement does not hold for the large $s$ case. This might be because, a few points can be adequate in SSF, additional points (i.e., a larger $s$) may have a small marginal benefit in variance reduction under the large $s$ setting. Consequently, the approximation error of SSF sometimes stays almost the same with a larger number of random features. QMC and ORF seek for variance reduction on random features. Nevertheless, they often work well in the large $s$ case. As demonstrated by the expression for variance of ORF and convergence rate in QMC , this theoretical result is consistent with the numerical performance of ORF and QMC, which may explain the reason why they work better in a large $s$ setting than a small $s$ case.

Results on arc-cosine kernels and polynomial kernels can be in Appendix B. Besides, apart from the above used medium/large scale datasets in our experiments, we also evaluate the compared approaches on a ultra-large scale dataset MNIST 8M with millions of data. Due to the memory limit, following the doubly stochastic framework , we incorporate these random features based approaches under the data streaming setting for the reduction of time and space complexity.

2.2 Classification results on MNIST and CIFAR10

Here we consider the MNIST and CIFAR10 datasets, on which we test these random features based algorithms for kernel approximation and then combine these algorithms with liblinear for image classification. In our experiment, we use the Gaussian kernelAs indicated by , (convolutional) NTK generally performs better than Gaussian kernel but it is still non-trivial to obtain a efficient random features mapping for (convolutional) NTK without much loss on prediction., whose kernel width $\varsigma$ is tuned by 5-fold cross validation over the grid $\varsigma=[0.01,0.1,1,10,100]$. For the MNIST database, we directly use the original 784-dimensional feature as the data. For better performance on the CIFAR10 dataset, we use VGG16 with batch normalization pre-trained on ImageNet as a feature extractor. We fine-tune this model on the CIFAR10 dataset with 240 epochs and a mini-batch size 64. The learning rate is initialized at 0.1 and then divided by 10 at the 120-th, 160-th, and 200-th epochs. For each color image, a 4096 dimensional feature vector is obtained from the output of the first fully-connected layer in this fine-tuned neural network.

Figure 6(a) shows the approximation error, the time cost (sec.), and the classification accuracy by liblinear across a range of $s=1000$ to $s=10,000$ random features on the MNIST database. We find that ORF and SSF yield the best approximation quality. Despite that most algorithms achieve different approximation errors, there is no significant difference in the test accuracy, which corresponds to the results on non-image datasets. Similar results are observed on the CIFAR10 dataset with $s=5000$ to $s=12,000$ random features; see Figure 6(b). Note that most algorithms take the similar time cost on generating random features except for the data-dependent algorithm LS-RFF. Several structured based approaches (e.g., Fastfood, SORF, ROM) do not achieve significant reduction on time cost due to the relatively inefficient Matlab built-in function to implement the Walsh-Hadamard transform.

In the previous sections, we review random features based algorithms and their theoretical results, that works under a fixed $d$ setting with $s\ll n$. Random features based approaches are simple in formulation but enjoy nice empirical validations and theoretical guarantees in kernel approximation and generalization properties. Recently, analysis of over-parameterized models has attracted a lot of attention in learning theory, partly due to the observation of several intriguing phenomena, including capability of fitting random labels, strong generalization performance of overfitted classifiers and double descent in the test error curve . Moreover, Belkin et al. point out that the above phenomena are not unique to deep networks but also exist in random features and random forests. In Figure 7, we report the empirical training error, the test error, and the kernel approximation error of random features regression as a function of $s/n$ on the sonar dataset and the MNIST dataset . Even with $n$, $d$, $s$ only in the hundreds, we can still observe that as $s$ increases, the training error reduces to zero and the approximation error monotonously decreases. However, the test error exhibits double descent, i.e., a phase transition at the interpolation threshold: moving away from this threshold on both sides trends to reduce the generalization error. This is somewhat striking as it goes against the conventional wisdom on bias-variance trade-off: predictors that generalize well should trade off the model complexity against training data fitting.

The above observations have motivated researchers to build on the elegant theory of random features to provide an analysis of neural networks in the over-parameterized regime. To be specific, RFF can be regarded as a two-layer (large-width) neural network, where the weights in the first layer are chosen randomly/fixed and only the output layer is optimized. This is a typical over-parameterized model if we take $s\gg n$. As such, two-layer neural networks in this regime are more amenable to theoretical analysis as compared to general arbitrary deep networks. This is a potentially fruitful research direction, and one hand, the optimization and generalization of such model have been studied in in deep learning theory. On the other hand, in order to explain the double descent curve of random features in over-parameterized regimes, we often work in a high dimensional setting, which is more subtle than classical results in standard settings, as indicated by recent random matrix theory (RMT) . An intuitive example is, $\|\bm{K}-\bm{Z}\bm{Z}^{\!\top}\|_{\operatorname{F}}\rightarrow 0$ always hold in low/high dimensions as $s\rightarrow\infty$ but $\|\bm{K}-\bm{Z}\bm{Z}^{\!\top}\|_{2}\rightarrow 0$ does not hold for $n,d,s\rightarrow\infty$. Accordingly, in this section, we provide an overview on analysis of (high dimensional) random features in over-parameterized setting, especially on double descent. We remark upfront that the random features model on double descent is not the only way for analyzing DNNs. Many other approaches, with different points of views, have been proposed for deep learning theory, but they are out of scope of this survey.

Here we briefly introduce the problem setting of high dimensional random features in over-parameterized regimes, and then discuss the techniques used in various studies.

In a similar spirit, Mei and Montanari use RMT to study the spectral distribution of the Gram matrix $\bm{Z}=\sigma(\bm{X}\bm{W}^{\!\top}/\sqrt{d})/\sqrt{d}$ by considering the Stieltjes transform of a related random block matrix, and show that, under least squares regression setting in an asymptotic viewpoint, both the bias and variance have a peak at the interpolation threshold $\psi_{2}=1$ and diverge there when $\lambda\rightarrow 0$. Under this framework, according to the randomness stemming from label noise, initialization, and training features, a refined bias-variance decomposition is conducted by and further improved by using the analysis of variance. Apart from refined error decomposition schemes, the authors of consider a general setting on convex loss functions, transformation matrix, and activation functions for regression and classification. Here the techniques used for analysis are not limited to RMT. Instead, replica method (a non-rigorous heuristic method from statistical physics) used in and the convex Gaussian min-max (CGMM) theorem used in are two alternative way to derive the desired results. Note that, CGMM requires the data to be Gaussian, which might restrict the application scope of their results but is still a common-used technical tool for max-margin linear classifier , boosting classifiers , and adversarial training for linear regression in over-parameterized regimes. Admittedly, the applied replica method in statistical physics is quite different from for tackling inverse random matrices in RMT. However, most of the above methods admit the equivalence between the considered data model and the Gaussian covariate model. That means, problem (3) with Gaussian data can be asymptotically equivalent to

2 Discussion on Random Features and DNNs

Admittedly, the above results may appear pessimistic due to the simple architecture. Nevertheless, random features is still an effective tool, at least the first step, for analyzing and understanding DNNs in certain regimes, and we believe its potential has yet to be fully exploited. Note that the random features model is still a strong and universal approximator in the sense that the RKHSs induced by a broad class of random features are dense in the space of continuous functions. While the aforementioned results show that the number of required features may be exponential in the worst case, a more refined analysis can still provide useful insights for DNNs. One potential way forward in deep learning theory is to use the random features model to analyze DNNs with pruning. For example, the best paper in the Seventh International Conference on Learning Representations (ICLR2019) put forward the following Lottery Ticket Hypothesis: a deep neural network with random initialization contains a small sub-network which, when trained in isolation, can compete with the performance of the original one. Malach et al. provide a stronger claim that a randomly-initialized and sufficiently over-parameterized neural network contains a sub-network with nearly the same accuracy as the original one, without any further training. Their analysis points to the equivalence between random features and the sub-network model. As such, the random features model is potentially useful for network pruning in terms of, e.g., guiding the design of neurons pruning for accelerating computations, and understanding network pruning and the full DNNs.

In this survey, we systematically review random features based algorithms and their associated theoretical results. We also give an overview on generalization properties of high dimensional random features in over-parameterized regimes on double descent, and discuss the limitations and potential of random features in the theory development for neural networks. Below we provide additional remarks and discuss several open problems that are of interest for future research.

As a typical data independent method, random features are simpler to implement, easy to parallelize, and naturally apply to streaming or dynamic data. Current efforts on Nyström approximation by a preconditioned gradient solver parallelized with multiple GPUs and quantum algorithms can guide us to design powerful implementation for random features to handling millions/billions data.

Experimental comparisons show that better kernel approximation does not directly translate to lower generalization errors. We partly answer this question in the current survey but it may be not sufficient to explain this phenomenon. We believe this issue deserves further in-depth study.

Kernel learning via the spectral density is a popular direction , which can be naturally combined with Generative Adversarial Networks (GANs); see for details. In this setting, one may associate the learned model with an implicit probability density that is flexible to characterize the relationships and similarities in the data. This is an interesting area for further research.

The double descent phenomenon has been observed and studied in random features model by various technical tools under different settings. Current theoretical results, such as those in , may be extended to a more general setting with less restricted assumptions on data generation, model formulation, and the target function. Besides, more refined analysis and delicate phenomena beyond double descent have been investigated on the linear model, e.g., multiple descent phenomena and optimal (negative) regularization . Understanding these more delicate phenomena for random features requires further investigation and refined analysis.

There exist significant gaps between the random features model and practical neural networks, both in theory and empirically. Even for fitting simple quadratic or mixture models, the random features model cannot achieve a zero error with a finite number of neurons in general, while NTK and fully trained networks can . Numerical experiments indicate that the prediction performance of NTK and CNTK may significantly degrade if the random features are generated from practically sized nets .

Despite the limitations of existing theory, random features models are still useful for understanding and improving DNNs. For example, understanding the equivalence between the random features model and weight pruning in the Lottery Ticket Hypothesis , may be promising future directions.

We hope that this survey will stimulate further research on the above open problems.

The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program / ERC Advanced Grant E-DUALITY (787960). This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information. This work was supported in part by Research Council KU Leuven: Optimization frameworks for deep kernel machines C14/18/068; Flemish Government: FWO projects: GOA4917N (Deep Restricted Kernel Machines: Methods and Foundations), PhD/Postdoc grant. This research received funding from the Flemish Government (AI Research Program). This work was supported in part by Ford KU Leuven Research Alliance Project KUL0076 (Stability analysis and performance improvement of deep reinforcement learning algorithms), EU H2020 ICT-48 Network TAILOR (Foundations of Trustworthy AI - Integrating Reasoning, Learning and Optimization), Leuven.AI Institute; and in part by the National Natural Science Foundation of China 61977046, in part by National Science Foundation grants CCF-1657420 and CCF-1704828, and in part by SJTU Global Strategic Partnership Fund (2020 SJTU-CORNELL) and Shanghai Municipal Science and Technology Major Project (2021SHZDZX0102).

where the first term in the right hand is the expected error difference between the original KRR and its random features approximation version. The second term in the right hand is the excess risk of KRR, which is independent of the quality of kernel approximation. Specifically, the first term can be further expressed by the representer theorem

I. the residual matrix $\bm{E}$ is semi-positive definite and $\widetilde{\bm{K}}_{1}$, $\widetilde{\bm{K}}_{2}$ are non-singular.

II. $n\lambda\leq\lambda_{1}(\widetilde{\bm{K}}_{1})\leq\lambda_{1}(\widetilde{\bm{K}}_{2})$, and $\widetilde{\bm{K}}_{2}$ admits (at least) polynomial decay.

which can be achieved by a geometry explanation in Figure 8. By virtue of Eq. (33) and Assumption II, we have

The left-hand of the above inequality can be further improved as

where $d_{\widetilde{\bm{K}}_{2}}^{\lambda}$ is the “effective dimension” of $\widetilde{\bm{K}}_{2}$ defined in Eq. (21) and the last inequality follows from Assumption II.

Appendix B Experiments

In this section, we detail the experimental settings and present the comparison results on the compared approaches on several benchmark datasets across various kernels. This part is organized as follows.

In Section B.1, we present experimental results across the Gaussian kernel on eight non-image datasets in terms of approximation error, the time cost (sec.) for generating random features mappings, classification accuracy by linear regression and liblinear.

Results on approximation error and test accuracy (by linear regression) across arc-cosine kernels and polynomial kernels are presented in Sections B.2 and B.3, respectively.

In Section B.4, a ultra-large scale dataset is applied to further validate the related algorithms.

Figures 9, 10 show the approximation error for the Gaussian kernel, the time cost (sec.) of generating randomized feature mappings, and the test accuracy yielded by linear regression and liblinear on the eight datasets, respectively. We see that as the number of random features increases, these algorithms achieve a smaller approximation error and a higher classification accuracy for both classifiers. We notice some interesting phenomena in terms of the relation between approximation quality and prediction performance, depending on whether the feature dimension is low (i.e., $s=2d$ or $s=4d$) or high (i.e., $s=16d$ or $s=32d$). In particular, the algorithms with the best kernel approximation performance are often different in the low-dimensional case and the high dimensional case. Therefore, no algorithm always dominate the others. On the other hand, while the approximation quality of these algorithms varies, their prediction performance are often similar. Further, to better understand the above observations, we summarize the best performing algorithm on each dataset in terms of the approximation quality and classification accuracy in Table VII, as illustrated in our main text (refer to Section 6.2.1).

Regarding to computational efficiency, most algorithms achieve the similar time cost on generating random features except SSF and LS-RFF. SSF requires constructing the transformation matrix by minimizing the discrete Riesz 0-energy in advance; LS-RFF is a data-dependent algorithm that needs to calculate the approximated ridge leverage score. Nevertheless, Fastfood/SORF/ROM does not achieve the reduction on time cost, which appears contradictory to the underlying theoretical result on time complexity. This might be because, one hand, the feature dimension of the used datasets in our experiments often ranges from 10 to 100, except for the image datasets. In this case, it appears difficult to observe the computational saving from $\mathcal{O}(sd)$ to $\mathcal{O}(s\log d)$ or $\mathcal{O}(d\log d)$. On the other hand, in our experiments, due to the relatively inefficient Matlab implementation of Fast Discrete Walsh-Hadamard Transform, typical algorithms (e.g., Fastfood/SORF/ROM) do not show a significant reduction on computational efficiency than RFF.

B.2 Results on Arc-cosine kernels

As mentioned before, according to Eq. (6), various algorithms based on different sampling strategies can be still applicable to arc-cosine kernels, e.g., ORF, QMC, and Fastfood. Accordingly, eight representative algorithms are taken into comparison on arc-cosine kernels, including RFF, ORF, SORF, ROM, Fastfood, QMC, SSF, and GQ.

Figures 5.2.2, 5.2.2 show the approximation error and test accuracy across the zero/first-order arc-cosine kernels, respectively. It can be observed that in most cases SSF and QMC achieve a lower approximation error than the other approaches, which corresponds to the theoretical findings. However, there is no distinct difference on approximation between RFF and ORF/SORF. In fact, the current theoretical results on ORF/SORF for variance reduction are only valid to the Gaussian kernel. Whether such results can be transferred to arc-cosine kernels are still unclear. In general, the approximation performance and time cost (see Figure 14(a) and 14(b)) of these algorithms on arc-cosine kernels are similar to that on the Gaussian kernel, though the approximation error value is often larger than that for the Gaussian kernel. This is because, according to Eq. (6), we actually conduct a $d$-dimensional integration approximation, the smoothness of the integrand $\sigma(\bm{\omega}^{\top}\bm{x})\sigma(\bm{\omega}^{\!\top}\bm{x}^{\prime})$ would significantly effect the approximation performance as indicated by sampling theory. In the term of classification performance, the difference in test accuracy of most algorithms is relatively small, which shows the similar tendency with that of the Gaussian kernel.

B.3 Results on Polynomial kernels

For polynomial kernel approximation, we include three representative approaches, tensorized random projections (TRP) , TensorSketch (TS) , and random Maclaurin (RM) sketch evaluated on eight datasets for approximation and prediction. Since the polynomial kernel can be written as a special type of tensor product, TS and TRP work in this setting by sketching a tensor product of arbitrary vectors, which is different from RM using Maclaurin expansion. Figure 5.2.2 shows that, TS and TRP have the similar test accuracy, but significantly perform better than RM, as RM’s generality is not required for the polynomial kernel. Besides, Figure 14(c) shows that RM is quite computational efficient due to its Maclaurin expansion scheme; while TS takes much time on generating random features since it utilizes a fixed sampling probability to compute the tensor sketch; while TRP works in a flexible sampling strategy proportional to its Maclaurin coefficient.

B.4 Results on the MNIST-8M dataset

Here we evaluate the compared ten algorithms across the Gaussian kernel and arc-cosine kernels on the MNIST-8M dataset . Due to the memory limit, following the doubly stochastic framework , we incorporate these random features based approaches under the data streaming setting for the reduction of time and space complexity. The experimental setting on this dataset follows with : the feature dimension $d=784$ is reduced to $100$ by PCA; the number of random features $s$ is set to $4096$; the used Gaussian RBF kernel with kernel bandwidth $\varsigma$ equaling to four times the median pairwise distance; logistic regression with the regularization parameter $\lambda=0.0005$ for this multi-class classification task; the batch size is set to be $2^{20}$ and feature block to be $2^{15}$. Besides, we report the total time cost of each algorithm on generating feature mapping, training process and test process for evaluation.

Table IX reports the approximation error, training error, test error, and the total time cost of each algorithm across the Gaussian kernel and the zero/first-order arc-cosine kernels under $s=4096$. It can be found that, ORF/SORF and SSF achieve the best approximation performance on the Gaussian kernel, but ORF fails to significantly improve the approximation ability on arc-cosine kernels. This is consistent with previous discussion on medium datasets in Section B.2.