We study the theoretical properties of random Fourier features classification with Lipschitz continuous loss functions such as support vector machine and logistic regression. Utilizing the regularity condition, we show for the first time that random Fourier features classification can achieve $O(1/\sqrt{n})$ learning rate with only $\Omega(\sqrt{n} \log n)$ features, as opposed to $\Omega(n)$ features suggested by previous results. Our study covers the standard feature sampling method for which we reduce the number of features required, as well as a problem-dependent sampling method which further reduces the number of features while still keeping the optimal generalization property. Moreover, we prove that the random Fourier features classification can obtain a fast $O(1/n)$ learning rate for both sampling schemes under Massart's low noise assumption. Our results demonstrate the potential effectiveness of random Fourier features approximation in reducing the computational complexity (roughly from $O(n^3)$ in time and $O(n^2)$ in space to $O(n^2)$ and $O(n\sqrt{n})$ respectively) without having to trade-off the statistical prediction accuracy. In addition, the achieved trade-off in our analysis is at least the same as the optimal results in the literature under the worst case scenario and significantly improves the optimal results under benign regularity conditions.

§ INTRODUCTION

Kernel methods have been widely used in many machine learning tasks such as regression and classification [Schölkopf and Smola, 2001, Schölkopf et al., 2004], as they provide a simple framework to model highly complicated functional relationships and well-established theoretical guarantees [Caponnetto and De Vito, 2007, Steinwart and Christmann, 2008]. The power of kernel methods comes from the so-called "kernel trick", where it utilizes a feature function $\phi(\cdot)$ to implicitly map the data into a high or possibly infinite dimensional feature space and thus allows non-linear functional learning. However, kernel methods are notorious for being time-consuming, since a typical kernel learning algorithm requires $O(n^3)$ computation and $O(n^2)$ memory, where $n$ is the number of training samples. Due to the prohibitive computational requirements, a flurry of research has been devoted to developing algorithms that efficiently approximate kernel functions [Smola and Schökopf, 2000, Williams and Seeger, 2001, Rahimi and Recht, 2007, Mahoney and Drineas, 2009, Alaoui and Mahoney, 2015, Rudi et al., 2017, Zhang et al., 2015].

Among many approximation frameworks, the random Fourier features (RFFs) method proposed by Rahimi and Recht, 2007 has received great attention recently <cit.>. The key idea of RFFs is to approximate the infinite dimensional feature map $\phi(\cdot)$ with an explicit $s$-dimensional random feature map $\phi_s(\cdot)$ through Bochner's theorem [Bochner, 1932, Rudin, 2017], which states that $\phi_s(\cdot)$ can be constructed through sampling from some spectral measure. Kernel methods are now reduced to linear learning in the feature space, which can be computed via fast linear solver [Shalev-Shwartz et al., 2011]. The computational cost decreases from roughly $O(n^3)$ in time and $O(n^2)$ in space to $O(ns^2)$ and $O(ns)$ respectively. As a result, significant computational savings can be achieved as long as $s \ll n$.

Despite their empirical success [Rahimi and Recht, 2007, Huang et al., 2014, Dai et al., 2014], theoretical understanding of the RFFs is incomplete. In particular, the question of how to choose $s$ in order to obtain the RFFs estimators with performance provably comparable to original kernel methods remains unclear. To this end, several authors study the properties of the RFFs to approximate the kernel function and the kernel Gram-matrix <cit.>. However, all of these works require $s = \Omega(n)$ features to guarantee no loss of prediction accuracy, which translates to no computational savings at all. A highly refined analysis in the context of ridge regression is proposed recently <cit.>. When the spectral measure is used for sampling, they first show that $O(\sqrt{n}\log n)$ features are adequate to guarantee the minimax optimal learning rate $O(1/\sqrt{n})$, the same learning rate obtained with full kernel ridge regression. Furthermore, they prove that the RFFs regression can obtain a fast learning rate at the expense of increasing the number of features. Finally, they demonstrate that using a problem-dependent sampling distribution can significantly reduce the number of features to $s = \Omega(1)$.

A question motivating our study is whether similar theoretical results hold in the classification setting where a key difference is the loss functions employed (Lipschitz continuous loss such as support vector machine and logistic regression). Bach, 2017, Sun et al., 2018, and Li et al., 2021 study the generalization properties of RFFs approximations in the classification setting. They show that RFFs estimators can provide computational gains while still preserving the statistical properties of the original kernel method. Nevertheless, a key requirement in these analyses is to employ a certain problem-dependent sampling distribution. Computing such a distribution often requires $O(n^3)$ in time and $O(n^2)$ in space already and hence is itself intractable. Therefore, whether RFFs classification can provide computational savings without using the problem-dependent sampling distribution remains unclear, and a detailed trade-off between the number of features required and the statistical prediction accuracy is still missing.

A key step in obtaining a better trade-off for RFFs regression is to employ the regularity condition (see Assumption A.$3$). However, this property is not used while analyzing the RFFs classification. In this paper, by incorporating the regularity condition, we improve the optimal results in the literature and provide a definitive answer to questions mentioned above by making the following contributions

* Under suitable regularity condition (Assumption A.$3$), Theorem <ref> shows that RFFs classification only requires $\Omega(\sqrt{n}\log n)$ features to guarantee the minimax optimal $O(1/\sqrt{n})$ learning rate, the same prediction accuracy as the original kernel classification methods. Our analysis allows the computational cost to reduce from $O(n^3)$ in time and $O(n^2)$ in space to $O(n^2)$ and $O(n\sqrt{n})$ respectively, and suggest that for a wide range of classification problems, RFFs approximations provide dramatic computational cost savings without loss of prediction accuracy. To the best of our knowledge, this is the first result confirming that such a computational gain is possible in the classification setting when the standard sampling method is used.

* Using Massart's low noise assumption (Assumption A.$4$), Theorem <ref> further provides a more refined analysis on the generalization properties of the RFF classification estimators. We obtain a sharp $O(1/n)$ learning rate for classification at the expense of more random features required.

* We also discuss how problem-dependent sampling distribution further reduces the computational cost in the $O(1/\sqrt{n})$ rate setting and the $O(1/n)$ rate setting. Our analysis expresses the trade-off between the number of features required and the statistical prediction accuracy in terms of the regularization parameter ($\lambda$) and the effective degree of freedom ($d(\lambda)$) and points out how utilizing the optimized feature can lead to a significant reduction in the computational cost.

* Finally, in Table <ref> and <ref>, we provide a comprehensive comparison between achieved results in this paper and the optimal bound in the literature. The analysis demonstrates that under benign conditions, our study obtains the sharpest bound on the number of features required in literature, while under worst case scenario, we match the optimal results in the literature.

§ BACKGROUND

§.§ Supervised Learning with Kernels

Let $P(x,y)=P_xP(y \mid x)$ be a joint probability density function defined on  $\mathcal{X}\times \mathcal{Y}$ where $\mathcal{X}$ is an instance space and $\mathcal{Y}$ a label space. While in regression tasks $\mathcal{Y} \subset \mathbb{R}$, in classification tasks it is typically the case that $\mathcal{Y}=\{-1, 1\}$. Let $\{(x_i,y_i)\}_{i=1}^n$  be a training set sampled independently from $P(x,y)$. The goal of a supervised learning defined with a kernel function $k$ (and the associated reproducing kernel Hilbert space $\mathcal{H}$) is to find a hypothesis $f \colon \mathcal{X} \rightarrow \mathcal{Y}$ such that $f \in \mathcal{H}$ and $f(x)$ is a good estimate of the label $y \in \mathcal{Y}$ corresponding to a previously unseen instance $x \in \mathcal{X}$. In particular, the learning can be formulated as the following optimization problem

f̂^λ         _f ∈ℋ  1/n∑_i=1^n l(y_i,f(x_i))+ λf_ℋ^2 .
where $l:\mathcal{Y}\times\mathcal{Y}\rightarrow \mathbb{R}_{+}$ is a loss function and $\lambda$ is the regularization parameter to prevent overfitting. As a result of the representer theorem [Schölkopf and Smola, 2001], an empirical risk minimization estimator in this setting can be expressed as $\hat{f}^{\lambda} =\sum_{i=1}^n \alpha_ik(x_i,\cdot)$ with $\alpha \in \mathbb{R}^n$ and the optimization problem can be reformulated as

α̂_k^λ         _α∈ℝ^n  1/n∑_i=1^n l(y_i,(𝐊α)_i) + λα^T 𝐊α ,

where $\mathbf{K}$ is the kernel Gram-matrix with $\mathbf{K}_{i,j} = k(x_i,x_j)$.

Learning Risk The hypothesis $\hat{f}^{\lambda}$ is an empirical estimator and we use the learning risk to assess its ability to capture the relationship between instances and labels given by $P$ [Caponnetto and De Vito, 2007]

\[\mathbb{E}_{P}[l_{\hat{f}^{\lambda}}] = \int_{\mathcal{X}\times \mathcal{Y}} l(y,\hat{f}^{\lambda}(x))dP(x,y)\ ,\]

where we use $l_{f}$ to denote $l(y,f(x))$. When the context is clear, we will omit $P$ from the expectation and write $\mathbb{E}[l_{\hat{f}^{\lambda}}]$.

The empirical distribution $P_n(x,y)$ is given by a sample of $n$ examples drawn independently from $P(x,y)$. The empirical risk is used to estimate the learning risk $\mathbb{E}[l_{\hat{f}^{\lambda}}]$ and it is given by
\[\vspace{-.1em} \mathbb{E}_n[l_{\hat{f}^{\lambda}}] = \frac{1}{n} \sum_{i=1}^n l(y_i,\hat{f}^{\lambda}(x_i))\ .\]
Similar to Rudi and Rosasco, 2017 and Caponnetto and De Vito, 2007, we will assume [The existence of $f_{\mathcal{H}}$ depends on the complexity of $\mathcal{H}$ which is related to the data distribution $P(y|x)$. For more details, please see Caponnetto and De Vito, 2007 and Rudi and Rosasco, 2017.] the existence of $f_\mathcal{H} \in \mathcal{H}$ such that $f_{\mathcal{H}} = \argmin_{f \in \mathcal{H}}\ \mathbb{E}[l_f]$.
Note that $\mathbb{E}[l_{f_{\mathcal{H}}}]$ is the lowest learning risk one can achieve in the reproducing kernel Hilbert space $\mathcal{H}$. Hence, theoretical studies of the estimator $\hat{f}^{\lambda}$ often concern how fast its learning risk $\mathbb{E}[l_{\hat{f}^{\lambda}}]$ converges to $\mathbb{E}[l_{f_{\mathcal{H}}}]$, that is, how fast the excess risk $\mathbb{E}[l_{\hat{f}^{\lambda}}] - \mathbb{E}[l_{f_{\mathcal{H}}}]$ converges to zero. In the remainder of the manuscript, we will refer to the rate at which the excess risk converges to zero as the learning rate.

§.§ Random Fourier Features

Despite providing a flexible non-linear approximation framework, kernel methods suffer from the scalability issue. In particular, kernel supervised learning often requires the store or the inverse of the kernel Gram matrix $\mathbf{K}$ ($O(n^2)$ and $O(n^3)$ computations respectively), which is prohibitive. As a result, many low-rank approximation algorithms have been designed to resolve this issue <cit.>.

Among them, RFFs method is a widely used, simple, and effective technique for scaling up kernel methods. The idea is due to Bochner's theorem [Bochner, 1932], which states that any bounded, continuous, and shift-invariant kernel is the Fourier transform of
a bounded positive measure. Assuming the spectral measure $d\tau$ has a density function $p(\cdot)$, we can write the corresponding kernel as
k(x,y)    =     ∫_𝒱e^-2πiv^T(x-y)dτ(v)

=    ∫_𝒱 (e^-2πi v^Tx)(e^-2πi v^Ty)^*p(v)dv  ,
where $c^{*}$ denotes the complex conjugate of $c \in \mathbb{C}$. Typically, the kernel is real valued and we can ignore the imaginary part <cit.>. Bach, 2017 and Rudi and Rosasco, 2017 further generalize the idea by considering the following decomposition of kernel functions
k(x,y) = ∫_𝒱ψ(v,x)ψ(v,y)p(v)dv  ,

where $\psi \colon \mathcal{V}\times \mathcal{X}\rightarrow \mathbb{R}$ is a continuous and bounded function with respect to $v$ and $x$. Hence, we can approximate the kernel function using its Monte-Carlo estimate
k̃(x,y)    =     1/s∑_i=1^sψ(v_i,x)ψ(v_i,y) ,

=    ϕ_s(x)^Tϕ_s(y)  .

where $\{v_i\}_{i=1}^s$ are sampled independently from the spectral measure $p(v)$ and

\[\phi_s(x) = \frac{1}{\sqrt{s}}[\psi(v_1,x),\dots,\psi(v_s,x)]^{\top}.\vspace{-.5em}\]
We denote the reproducing kernel Hilbert space spanned by $\tilde{k}$ as $\tilde{\mathcal{H}}$ (note that in general $\tilde{\mathcal{H}} \nsubseteq \mathcal{H}$). Let $\tilde{\mathbf{K}}$ be Gram-matrices with entries $\tilde{\mathbf{K}}_{ij} = \tilde{k}(x_i,x_j)$. Then the following equalities can be derived easily from Eq. (<ref>)

\begin{align*}
\begin{aligned}
k(x,y) = \mathbb{E}_{v \sim p}\big[\tilde{k}(x,y) \big]\quad \wedge \quad \mathbf{K} = \mathbb{E}_{v\sim p}[\tilde{\mathbf{K}}] \ .
\end{aligned}
\end{align*}

In addition to the kernel Gram-matrix approximation, Bach, 2017 establishes that any $f \in \mathcal{H}$ can be expressed as [It is not necessarily true that for any $g \in L_2(d\tau)$, there exists a corresponding $f \in \mathcal{H}$.]
f(x) = ∫_𝒱g(v)ψ(v,x)p(v)dv     (∀x ∈𝒳)

where $g \in L_2(d\tau)$ is a real-valued function such that $\|g\|_{L_2(d\tau)}^2 < \infty$ and $\|f\|_{\mathcal{H}}= \min_{g} \|g\|_{L_2(d\tau)}$, with the minimum taken over all possible decompositions of $f$. Thus, one can take an independent sample $\{v_i\}_{i=1}^s \sim p(v)$ (we refer to this sampling scheme as plain RFF) and approximate a function $f \in \mathcal{H}$ by an element from $\tilde{\mathcal{H}}$ as

\begin{align*}
\tilde{f}(\cdot) = \sum_{i=1}^s\alpha_i \psi(v_i,\cdot) = \phi_{s}(\cdot)^{\top}\alpha \quad \text{with} \quad \alpha \in \mathbb{R}^s\ .
\end{align*}

As the latter approximation is simply a Monte Carlo estimate, one could also select an importance weighted probability density function $q(\cdot)$ and sample features $\{v_i\}_{i=1}^s$ from $q$ (we refer to this sampling scheme as weighted RFF).
The function $f$ can then be approximated by

\begin{align*}
\tilde{f}_{q}(\cdot) = \sum_{i=1}^s\alpha_i \psi_q(v_i,\cdot) = \phi_{q,s}(\cdot)^{\top}\alpha \ ,
\end{align*}

with $\psi_q(v_i,\cdot)= \sqrt{p(v_i)/q(v_i)}\psi(v_i,\cdot)$ and $\phi_{q,s}(\cdot) = (1/\sqrt{s})[\psi_q(v_1,\cdot),\cdots,\psi_q(v_s,\cdot)]^{\top}$.

For both plain RFF and weighted RFF, the goal is to find $\tilde{f}$ with minimal norm such that the computation error between $\tilde{f}$ and $f$ is minimized. Similar to Bach, 2017, the RFFs sampling can be formulated as the following optimization problem
f̃ - f_L_2(P_x)^2 + λf̃_ℋ̃^2   .

Note that since $\tilde{f} \notin \mathcal{H}$ in general, we use the $L_2(P_x)$ norm to measure the computation error.

§.§ Integral Operator & Leverage Score Sampling

Kernel methods and RFFs are often studied through the integral operator $L:L_2(P_x) \rightarrow L_2(P_x)$, which we define below
(Lf)(·) = ∫_𝒳k(x,·)f(x)dP_x(x) .
Given the kernel decomposition as Eq. (<ref>), the integral operator can be expressed as an expectation [Bach, 2017]
Lf    =    ∫_𝒳k(x,·)f(x)dP_x(x) ,

=    (∫_𝒱ψ(v,·)⊗ψ(v,·)p(v)dv)f   ,

where $f\otimes g$ is the $L_2(P_x)$ outer product operator such that $\left(f\otimes g\right)h = \langle g, h \rangle_{L_2(P_x)}f$. Finally, if $k$ and $\psi$ are both bounded and continuous, then $L$ is positive definite, self-adjoint and trace-class. In particular, if $\|\psi\| \leq \kappa$, we have $\|L\| \leq \kappa^2$.

Similarly, for kernel $\tilde{k}$, we define the integral operator $L_s:L_2(P_x) \rightarrow L_2(P_x)$:
L_sf    =    ∫_𝒳k̃(x,·)f(x)dP_x(x) ,

=    ∫_𝒳1/s∑_i=1^sψ(v_i,·)ψ(v_i,x)f(x)dP_x  ,

=    (1/s∑_i=1^sψ(v_i,·)⊗ψ(v_i,·))f   .

Hence, $L_s$ can be seen as an empirical estimator of $L$.

The study of the integral operator is important because it provides information on how to select the optimal sampling distribution $q(v)$. A large body of literature shows that finding an optimal sampling distribution $q(v)$ often significantly reduces the number of features required [Bach, 2017, Alaoui and Mahoney, 2015, Avron et al., 2017, Rudi and Rosasco, 2017]. The reason is that random features sampled according to $p(v)$ often focus on approximating the leading eigenvalues of the integral operator $L$. In contrast, a reweighted sampling distribution $q(v)$ allows the random features to span the whole eigenspectrum of $L$.

In light of this, a leverage score based weighted distribution function is first introduced in Alaoui and Mahoney, 2015 in the context of the Nyström approximation [Nyström, 1930, Smola and Schölkopf, 2000, Williams and Seeger, 2001]. Utilizing the importance reweighted nature, Alaoui and Mahoney, 2015 establish a sharp convergence rate of the low-rank estimator based on the Nyström method.

The success of the leverage score distribution further motivates the pursuit of a similar notion for RFFs. In particular, Bach, 2017 first proposes the leverage score sampling based on a leverage score function defined below
τ_λ(v) = p(v)⟨ψ(v,·), (L+λI)^-1ψ(v,·)⟩_L_2(P_x)  .

From our assumption, it follows that there exists a constant $\kappa$ such that $| \psi(v,x) | \leq \kappa$ (for all $v$ and $x$). We now have

\begin{align*}
\begin{aligned}
\tau_{\lambda}(v)\leq p(v)\frac{\kappa^2}{\lambda} \ .
\end{aligned}
\end{align*}

An important property of function $\tau_{\lambda}(v)$ is its relation to the effective number of parameters:

\begin{align*}
\begin{aligned}
\int_{\mathcal{V}}\tau_{\lambda}(v)dv = \text{Tr}\big[L(L+\lambda I)^{-1} \big]:= d(\lambda) \ ,
\end{aligned}
\end{align*}

where $d(\lambda)$ implicitly determines the number of parameters in a supervised learning problem and is thus called the number of effective degrees of freedom [Bach, 2013, Hastie, 2017].

We can now sample features according to $q^*(v) = \tau_{\lambda}(v)/d(\lambda)$, since $q^*(v)$ is a probability density function. Bach, 2017 studies the property of $q^*(v)$ and demonstrates that sampling according to $q^*(v)$ requires fewer Fourier features compared to the standard spectral measure sampling. From now on, we refer to $q^*(v)$ as the ridge leverage score distribution and refer to this sampling strategy as leverage weighted RFF.

§ MAIN RESULTS

In this section, we provide our theoretical analysis on the trade-off between the number of random features and the statistical prediction accuracy. We first discuss the worst case scenario where the estimator achieves the $O(1/\sqrt{n})$ learning rate, followed by demonstrating the trade-off in the fast convergence rate setting.

§.§ $O(1/\sqrt{n})$ Learning Rate

We study the scenario where the RFFs estimator obtains the minimax learning rate $O(1/\sqrt{n})$. As discussed before, kernel supervised learning can be formulated as Eq. (<ref>). Since we are investigating the classification setting, we mainly consider the loss function $l$ to be uniformly Lipschitz continuous functions such as support vector machine and logistic regression. A fatal problem for kernel supervised learning is the computational cost since kernel learning problem such as Eq. (<ref>) often requires the store of the kernel Gram matrix $\mathbf{K}$ or even the inversion of $\mathbf{K}$, which are $O(n^2)$ and $O(n^3)$ computations respectively.

In order to overcome the computation issue, the RFFs provide an efficient way to approximate the kernel function. Specifically, we sample $v_1,\dots,v_s$ according to some importance sampling distribution $q(v)$ to form the random feature vector $\phi_s(\cdot)$. For a given data $(x,y)$, we then approximate the label $y$ with the random feature hypothesis $\tilde{f}_q(x) = \phi_{q,s}(x)^{\top}\beta$. The RFFs learning can be cast as the following optimization problem
f̃^λ         _f̃_q ∈ℋ̃  1/n∑_i=1^n l(y_i,f̃_q(x_i))+ λf̃_ℋ̃^2 .
According to Bach, 2017 and Li et al., 2019, we have $\|\tilde{f}\|_{\tilde{\mathcal{H}}}^2 \leq \|\beta\|_2^2$, as a result, the above optimization can be reformulated as
β̃^λ         _β∈ℝ^s  1/n∑_i=1^n l(y_i,ϕ_q,s^⊤β)+ λβ_2^2 .

The RFFs hypothesis with loss function $l$ can be represented as $\tilde{f}^{\lambda} = \phi_{q,s}^{\top}\tilde{\beta}^{\lambda}$. Through the RFFs approximation, we now only need to store the feature matrix $\Phi_q = [\phi_{q,s}(x_1),\dots,\phi_{q,s}(x_n)]^{\top} \in \mathbb{R}^{n\times s}$. The inversion of $\mathbf{K}$ can be approximated as inverting $\Phi_q^{\top}\Phi_q \in \mathbb{R}^{s\times s}$. Hence the computation cost is now $O(ns)$ and $O(ns^2+s^3)$ respectively. We can see that if $s \ll n$, RFFs method enjoys a huge computational savings. However, a key question is how the choice of $s$ affects the prediction accuracy of $\tilde{f}^{\lambda}$.

In this section, we try to address the above issue. We first list our assumptions below

A.$1$ We assume that the kernel has integral expansion as Eq. (<ref>) such that $\psi(v,x)$ is continuous in both $v$ and $x$ and $|\psi(v,x)| \leq \kappa$ for all $x \in \mathcal{X}$ and $v \in \mathcal{V}$;

A.$2$ Assume that the loss function $l$ in Eq. (<ref>) is uniformly Lipschitz continuous with constant $M$, i.e.,
\[\left|l(y,x_1)-l(y,x_2)\right| \leq M\|x_1-x_2\|_2.\]

A.$3$ Recall $f_{\mathcal{H}} = \argmin_{f \in \mathcal{H}}\ \mathbb{E}[l_f]$, we assume that
\[f_{\mathcal{H}} = L^rg, ~~\textnormal{for~some~} r \in [1/2,1]~\&~ g \in L_2(P_x);\]

Assumptions A.$1$ and A.$2$ are standard assumptions made in classification problems. A.$3$ is a regularity condition that is commonly used in approximation theory [Smale and Zhou, 2003]. It describes the decay rate of the coefficients of $f_{\mathcal{H}}$ along the basis given by the integral operator $L$, which further allows controlling the bias of the estimator. While being overlooked in the classification setting, A.$3$ is a key property used in RFFs regression to obtain a better computation and accuracy trade-off. Utlizing A.$3$ enables us to prove the following refined analysis.

Assume A.$1$, A.$2$ and A.$3$ hold. Suppose we have a measurable function $\tilde{\tau}: \mathcal{V} \rightarrow \mathbb{R}$ such that $\tilde{\tau}(v) \geq \tau_{\lambda}(v)$ almost surely. Denote $d_{\tilde{\tau}} = \int_{\mathcal{V}}\tilde{\tau}(v) dv$, and let $q(v) = \frac{\tilde{\tau}(v)}{d_{\tilde{\tau}}}$. We sample $v_1,\dots,v_s \sim q(v)$ and compute the hypothesis $\tilde{f}^{\lambda}$ by solving the optimization problem in Eq. (<ref>). Let $\delta \in (0,1)$, if we have

\[s \geq 12 d_{\tilde{\tau}}\log\frac{d(\lambda)}{\delta},\]
with probability over $1-\delta$,
𝔼( l_f̃^λ) - 𝔼 ( f_ℋ) ≤2MRλ^r + O(1/√(n))  .

Theorem <ref> expresses the trade-off between the computational cost and statistical efficiency through the regularization parameter $\lambda$, the effective dimension of the problem $d(\lambda)$, and the normalization constant $d_{\tilde{\tau}}$ of the sampling distribution. The regularization parameter $\lambda$ is used as a key quantity in the analysis of supervised learning setting [Caponnetto and De Vito, 2007, Rudi and Rosasco, 2017, Li et al., 2019]. In particular, if we set $\lambda \propto 1/n^{2r}$, we observe that the estimator $\tilde{f}^{\lambda}$ attains the $O(1/\sqrt{n})$ learning rate [Bach, 2017]. As a consequence of Theorem <ref>, we have the following bounds on the number of required features for the two strategies: plain RFF (Corollary <ref>) and leverage weighted RFF (Corollary <ref>).

If the probability density function from Theorem <ref> is the spectral measure $p(v)$, then the upper bound on the learning risk from Eq. (<ref>) holds for all $s \geq 5\kappa^2/\lambda\log\frac{16d(\lambda)}{\delta}$.

We set $\tilde{l}(v) = p(v)\kappa^2/\lambda$ and obtain $d_{\tilde{l}} = \int_{\mathcal{V}}p(v)\kappa^2/\lambda dv = \kappa^2/\lambda$.

Theorem <ref> and Corollary <ref> have several implications on the choice of $\lambda$ and $s$ in the classification setting with plain RFF. In particular, the usual generalization bound for kernel estimator $\hat{f}^{\lambda}$ (i.e., minimizer of Eq. (<ref>)) is $O(1/\sqrt{n})$ <cit.>. As such, if we set $\lambda = O(n^{1/2r})$, we can see that the RFFs estimator $\tilde{f}^{\lambda}$ incurs no loss of prediction accuracy while offering computational gains.

Specifically, in the benign case where $r = 1$, $O(\sqrt{n}\log n)$[We use the fact that $d(\lambda) \ll n$] features is able to achieve the $O(1/\sqrt{n})$ learning rate. Comparing with the existing analysis where $O(n\log n)$ features are required [Rahimi and Recht, 2009, Li et al., 2019], our result is a significant improvement. We also achieve remarkable computational savings: from roughly $O(n^3)$ and $O(n^2)$ in time and space for original kernel methods to $O(n^2)$ and $O(n\sqrt{n})$ for the RFFs approximation. Moreover, when $r > 1/2$, we also obtain computational gain as the number of features required now is $\Omega(n^{1/2r})$ with $2r >1$. In the worst scenario where $r = 1/2$ (equivalent to assuming $f_{\mathcal{H}}$ exists), we recover the results from existing analysis [Rahimi and Recht, 2007, Li et al., 2019].

To our knowledge, this is the first result showing that for a large class of classification problems ($r > 1/2$), RFFs classification can dramatically reduce the computational cost while preserving the optimal generalization properties.

If the probability density function from Theorem <ref> is the ridge leverage score distribution $q^*(v)$, the upper bound on the risk from Eq. (<ref>) holds for all $s \geq  5d(\lambda)\log\frac{16d(\lambda)}{\delta}$.

For this corollary, we set $\tilde{\tau}(v) = \tau_{\lambda}(v)$ and deduce $d_{\tilde{l}} = \int_{\mathcal{V}}\tau_{\lambda}(v)dv = d(\lambda)$.

Corollary <ref> details the number of features required in the leverage weighted RFF setting. Similar to the plain RFF setting, the RFFs estimator obtains $O(1/\sqrt{n})$ rate once we set $\lambda = O( n^{1/2r})$. However, the choice of $s$ now is determined by two factors: the regularity condition $r$ and the decay rate of the eigenspectrum of $L$.

We first consider the benign scenario where $r = 1$. Depending on the eigenspectrum decay rate, we have several different cases. Denote $\{\mu_1,\mu_2,\dots,\}$ to be the eigenvalue of $L$, in the best case where $L$ has finite rank, $d(\lambda)$ remains constant as $n$ grows. We therefore conclude that even $\Omega(1)$ features can guarantee the $O(1/\sqrt{n})$ learning rate. Next, if the eigenspectrum displays exponential decay, i.e., $\mu_i \propto C_0r^i$, we have $d(\lambda) \leq \log (C_0^0/\lambda)$. We can see that $s \geq \log n \log \log n$ is enough to achieve the $O(1/\sqrt{n})$ learning rate. As such, significant computational savings is obtained: from $O(n^3)$ and $O(n^2)$ to $O(n\log^4 n)$ and $O(n \log n)$ respectively. In the case of a slower decay with $\mu_i \propto C_0 i^{-2\gamma}$, we have $d(\lambda) \leq (R_0/\lambda)^{1/(2\gamma)}$ and $s\geq n^{1/4\gamma}\log n$. Hence, substantial computational savings can be achieved even in this case. Furthermore, in the worst case with $\mu_i$ close to $C_0 i^{-1}$, our bound implies that $s \geq n^{1/2}\log n$ features are sufficient.

The analysis for the worst case scenario where $r = 1/2$ is similar. The required numbers of features are $\Omega(1), \Omega(\log n \log \log n), \Omega(n^{1/2\gamma})$ and $\Omega(n\log n)$ for the cases where the eigenspectrum has finite rank, decays exponentially, proportional to $C_0 i^{-2\gamma}$ and close to $C_0i^{-1}$, respectively. Our results demonstrate that huge computational savings are possible as long as the eigenspectrum of $L$ displays fast decay (faster than $i^{-1}$).

sampling scheme   spectrum     this work    li et al. (2021)    learning rate

37emplain rff     finite rank     $s \in \Omega (n^{1/2r})$    $s \in \Omega (n)$     33em$O(1/\sqrt{n})$

exponential decay    $ s \in \Omega(n^{1/2r} \cdot\log \log n)$    $ s \in \Omega(n\cdot\log \log n)$

polynomial decay    $s \in \Omega(n^{1/2r}\cdot \log n)$    $s \in \Omega(n\cdot \log n)$

47emweighted rff    finite rank     $s \in \Omega(1)$    $s \in \Omega(1)$    43em$O(1/\sqrt{n})$

Exponential Decay     $s \in \Omega (\log n \cdot \log \log n)$    $s \in \Omega (\log n \cdot \log \log n)$

$\mu_i \propto i^{-2\gamma}, \gamma \geq 1$     $s \in \Omega (n^{1/4\gamma r} \cdot \log n)$    $s \in \Omega (n^{1/2\gamma} \cdot \log n)$

$\mu_i \propto i^{-1}$     $s \in \Omega (n^{1/2r} \cdot \log n)$    $s \in \Omega (n \cdot \log n)$

The comparison of our results to the sharpest learning rates from prior work [Li et al., 2021], where $r \in [1/2,1]$  .

§.§.§ Comparison with Existing Sharpest Results

Under $O(1/\sqrt{n})$ learning rate setting, Rahimi and Recht, 2009, Bach, 2017, and Li et al., 2021 analyze the trade-off between the number of features and the statistical prediction accuracy. Table <ref> provides a detailed comparison between this work and that from Li et al., 2021. Our results show that we obtain at least the same rate as the previous best rate in the literature when $r = 1/2$, while significantly improving the trade-off under benign conditions ($r > 1/2$).

The first block of rows in Table <ref> illustrates the difference between our work and that from Li et al., 2021 when plain sampling is used. A key feature in our results is that when $r > 1/2$, our results state that the RFFs approximation provides computational gain without trading off for the prediction accuracy. In comparison, results from Li et al., 2021 state that there is no computational gain ($s = \Omega(n)$) if we were to achieve the $O(1/\sqrt{n})$ learning rate. In addition, we recover the results from Li et al., 2021 when $r = 1/2$.

We observe a similar pattern when leverage weighted RFF is used. In particular, our results match those from Li et al., 2021, when the eigenspectrum has finite rank or displays exponential decay. However, as soon as the eigenspectrum has polynomial decay, our result is sharper. Specifically, when the eigenvalue decays polynomially with $\mu_i \propto i^{-2\gamma}$ and $r = 1$, our results show that $\Omega(n^{1/4\gamma \log n})$ features are enough to achieve $O(1/\sqrt{n})$ learning rate, comparing with $\Omega(n^{1/2\gamma})$ features required from Li et al., 2021. When the eigenvalue decays close to $i^{-1}$, our results require $\Omega(n^{1/2}\log n)$ features while Li et al., 2021 require $\Omega(n \log n)$ features.

§.§ Refined Learning Rate

In the previous section, we study the trade-off between the number of features and the statistical prediction accuracy in the $O(1/\sqrt{n})$ minimax learning rate setting. In general, it is hard for classification problems to obtain learning rates sharper than $O\left(1/\sqrt{n} \right)$. However, under some benign conditions, it is possible to obtain $O(1/n)$ convergence rate as demonstrated by Bartlett et al., 2006 and Steinwart and Christmann, 2008. As such, with the help of the following assumption, we derive a sharp learning rate for RFFs classification problems in this section.

A.$4$ Recall that $f_{\mathcal{H}}$ is the optimal estimator in A.$3$. We assume that there exists a constant $G$ such that for all $f \in \mathcal{H}$
\[\mathbb{E}[(f - f_{\mathcal{H}})^2] \leq G\mathbb{E}[l_f - l_{f_{\mathcal{H}}}] \ .\]

Assumption A.$4$ is a widely used condition for classification problems to obtain faster learning rates. It typically requires that the loss function $l$ is uniformly convex and the function space $\mathcal{H}$ is convex and uniformly bounded. It can be shown that many loss functions satisfy this assumption, including squared loss [Bartlett et al., 2005] and hinge loss <cit.>. Additional examples of these loss functions are discussed in Bartlett et al., 2006 and Mendelson, 2002. In addition, since $l$ is Lipschitz continuous, we can rewrite A.$4$ as
\[\mathbb{E}[(l_f-l_{f_{\mathcal{H}}})^2] \leq L^2\mathbb{E}[(f - f_{\mathcal{H}})^2] \leq G L^2\mathbb{E}[l_f - l_{f_{\mathcal{H}}}] \ .\]
This is the variance condition described in <cit.>, which is also linked to the Massart's low noise condition or more generally to the Tsybakov condition [Sun et al., 2018]. Intuitively speaking, the condition requires that the bayes classifier $P(Y=1 \mid X=x)$ is not close to $1/2$ <cit.>.

Assume A.$1$-A.$4$ hold. In addition, we assume the condition for $\tilde{\tau}$, $d_{\tilde{\tau}}$ and $q(v)$ hold as that in Theorem <ref>. Let $\{\tilde{\mu}_1,\tilde{\mu}_2,\dots\}$ be the eigenvalues of the normalized Gram-matrix $(1/n) \tilde{\mathbf{K}}$, $c_1, c_2, c_3$ be some universal constant, and $\delta \in (0,1)$, if we have

\[s \geq 12 d_{\tilde{\tau}}\log\frac{d(\lambda)}{\delta},\]
with probability over $1-\delta$,

𝔼( l_f̃^λ) - 𝔼 ( f_ℋ) ≤2MRλ^r + c_1r̂^* + c_2/nlog1/δ   ,

\[\hat{r}^* \leq c_3 \min _{0 \leq h \leq n}\left(\frac{h}{n}+\sqrt{\frac{1}{n} \sum_{i>h} \tilde{\mu}_{i}}\right).\]

Theorem <ref> covers a wide range of cases and can provide sharp risk convergence rates. In particular, $\hat{r}^*$ has an upper bound of $O(1/\sqrt{n})$ in all cases, which happens when $\hat{\mu}_i$ decays polynomially as $O(n^{-\gamma})$ with $\gamma > 1$ and we let $h = 0$. On the other hand, if $\hat{\mu}_i$ decays exponentially, then setting $h = \ceil{\log n}$ implies that $\hat{r}^* \leq O(\log n /n)$. In the best case, when $(1/n) \mathbf{K}$ has only finite rank $d_K$, then $\hat{r}^* \leq O(1/n)$ by letting $h = d_K+1$. These different upper bounds provide insights into various trade-offs between computational complexity and statistical efficiency. We now split the discussion into two cases: plain RFF and leverage weighted RFF.

Under plain RFF strategy, similar to Corollary <ref>, we have $d_{\tilde{\tau}} \leq \kappa^2/\lambda$. If the eigenvalues decay polynomially, i.e., $\mu_i \propto i^{-\gamma}$ with $\gamma > 1$, then the learning rate is upper bounded by $O(1/\sqrt{n})$. In this case, we need $s =\Omega(n^{1/2r}\log n)$. On the other hand, if $\mu_i$ decays exponentially, $\Omega(n^{1/r})$ features are able to guarantee $O(\log n / n)$ learning rate. Finally, if the eigenspectrum has finite rank, $\Omega(n^{1/r})$ features yield $O(1/n)$ fast learning rate. Under the leverage weighted RFF, the required numbers of features and the corresponding learning rates for the three above cases are: i) $s =\Omega(n^{1/4\gamma r} \log n)$ and $O(1/\sqrt{n})$ (polynomial decay $\mu_i \propto i^{-\gamma}$ with $\gamma > 1$); ii) $s=\Omega(\log n \log \log n)$ and $O(\log n/n)$ (exponential decay); and iii) $s = \Omega(1)$ and $O(1/n)$ (finite rank).

sampling scheme   results   spectrum     number of features     learning rate

67emplain rff   36emthis work
finite rank     $s \in \Omega (n^{1/r})$    $O(1/n)$

exponential decay     $ s \in \Omega(n^{1/r})$    $O(\log n / n)$

$\mu_i \propto i^{-\gamma}$     $s \in\Omega(n^{1/2r}\cdot \log n)$     $O(1/\sqrt{n})$

36emli et al. (2021)     finite rank    $s \in\Omega(n^2)$     $O(1/n)$

exponential decay     $s \in\Omega(n^2)$     $O(\log n /n)$

$\mu_i \propto i^{-\gamma}$      $s \in\Omega(n\cdot \log n)$    $O(1/\sqrt{n})$

67emweighted rff   36emthis work    finite rank     $s \in \Omega(1)$    $O(1/n)$

exponential decay     $s \in \Omega (\log n \cdot \log \log n)$     $O(\log n / n)$

$\mu_i \propto i^{-\gamma}$     $s \in \Omega (n^{1/4\gamma r} \cdot \log n)$     $O(1/\sqrt{n})$

36emsun (2018)     finite rank     $s \in \Omega (1)$    $O(1/n)$

exponential decay    $s \in \Omega (\log^d n \cdot \log \log^d n)$      $O(\log^{d+2} n / n)$

$\mu_i \propto i^{-\gamma}$    $s \in \Omega (n^{\frac{2}{2+\gamma}} \cdot \log n)$      $O(1/n^{\frac{\gamma}{2+\gamma}})$

The comparison of our results to the sharpest results in the literature under fast learning rate setting, where $r \in [1/2,1]$.

§.§.§ Comparison with Existing Sharpest Results

For fast learning rate scenario, Li et al., 2021 and Sun et al., 2018 both study the trade-off between the number of features and the prediction accuracy and obtain a similar bound on the number of features required. Table <ref> compare our results with that in Li et al., 2021 under plain RFF sampling, and that in Sun et al., 2018 under leverage weighted RFF sampling. Similar to the analysis in the $O(1/\sqrt{n})$ scenario, our results strictly dominate previous optimal results in Li et al., 2021 when $r > 1/2$ and match the obtained bound in Li et al., 2021 when $r = 1/2$, because of Assumption A.$3$.

Under weighted RFF, our results match that in Sun et al., 2018 when the eigenspectrum has finite rank. However, when the eigenspectrum displays exponential decay, results from Sun et al., 2018 suffer from the curse of dimension, since both the number of features required and the learning rate obtained depend on the data dimension $d$. In contrast, our analysis does not have this dependency. When the eigenspectrum exhibits a polynomial decay, our results achieve the $O(1/\sqrt{n})$ learning rate while Sun et al., 2018 obtain a more flexible rate that depends on $\gamma$. If $\gamma \leq 2$, our results have a better trade-off as both the number of features and the learning rate are sharper than those from Sun et al., 2018. For example, if $\gamma = 2$, analysis in Sun et al., 2018 shows that $\Omega(n^{1/2}\log n)$ features can obtain $O(n^{-1/2})$ learning rate, whereas our results state that $\Omega(n^{1/8r}\log n )$ features yield the same learning rate. When $\gamma > 2$, the learning rate obtained by Sun et al., 2018 is faster than ours at the cost of increasing the number of features, i.e., $ \Omega (n^{\frac{2}{2+t}} \cdot \log n)$ versus $ \Omega(n^{\frac{1}{4\gamma r}} \cdot \log n)$. In particular, setting $\gamma=4$, we can see that Sun et al., 2018 obtain a fast $O(n^{2/3})$ learning rate, but at the cost of requiring $ \Omega(n^{1/3}\cdot \log n)$ random features. For the same setting, on the other hand, we obtain the minimax optimal $O(1/\sqrt{n})$ learning rate with only $\Omega(n^{1/16r}\cdot \log n)$ random features.

§ CONCLUSION

In this paper, we thoroughly study the generalization properties of RFFs classification with Lipschitz continuous loss such as support vector machine and logistic regression. Our main results for the first time demonstrate that RFFs classification can indeed provide computational gains without hurting the prediction accuracy when plain RFF is used. This is in contrast with all previous results that suggest that computational savings come at the expense of prediction accuracy. Furthermore, our analysis shows that a fast $O(1/n)$ learning rate is possible at the cost of increasing the number of features unless the leverage weighted RFF is used. However, a limitation in our work is that in the worst case where $r = 1/2$, we can see that the current analysis on RFFs classification cannot guarantee a computational gain. Therefore, how to obtain a sharper result in the worst case is an interesting future direction. In addition, how to efficiently approximate the leverage score is also an important future direction since it can often leads to significant reduction in the required number of feature as well as computational cost.

[Alaoui and Mahoney, 2015]
Ahmed Alaoui and Michael W Mahoney.
Fast randomized kernel ridge regression with statistical guarantees.
In Advances in Neural Information Processing Systems, pages
775–783, 2015.

[Avron et al., 2017]
Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya
Velingker, and Amir Zandieh.
Random Fourier features for kernel ridge regression: Approximation
bounds and statistical guarantees.
In International Conference on Machine Learning, pages
253–262, 2017.

[Bach, 2013]
Francis Bach.
Sharp analysis of low-rank kernel matrix approximations.
In Conference on Learning Theory, pages 185–209, 2013.

[Bach, 2017]
Francis Bach.
On the equivalence between kernel quadrature rules and random feature
Journal of Machine Learning Research, 180
(21):0 1–38, 2017.

[Bartlett and Mendelson, 2002]
Peter L Bartlett and Shahar Mendelson.
Rademacher and Gaussian complexities: Risk bounds and structural
Journal of Machine Learning Research, 30
(Nov):0 463–482, 2002.

[Bartlett et al., 2005]
Peter L Bartlett, Olivier Bousquet, Shahar Mendelson, et al.
Local Rademacher complexities.
The Annals of Statistics, 330 (4):0
1497–1537, 2005.

[Bartlett et al., 2006]
Peter L Bartlett, Michael I Jordan, and Jon D McAuliffe.
Convexity, classification, and risk bounds.
Journal of the American Statistical Association, 1010
(473):0 138–156, 2006.

[Bochner, 1932]
Salomon Bochner.
Vorlesungen über Fouriersche Integrale.
In Akademische Verlagsgesellschaft, 1932.

[Caponnetto and De Vito, 2007]
Andrea Caponnetto and Ernesto De Vito.
Optimal rates for the regularized least-squares algorithm.
Foundations of Computational Mathematics, 70
(3):0 331–368, 2007.

[Dai et al., 2014]
Bo Dai, Bo Xie, Niao He, Yingyu Liang, Anant Raj, Maria-Florina F Balcan, and
Le Song.
Scalable kernel methods via doubly stochastic gradients.
Advances in Neural Information Processing Systems,
27:0 3041–3049, 2014.

[Hastie, 2017]
Trevor J Hastie.
Generalized additive models.
In Statistical models in S, pages 249–307. Routledge, 2017.

[Huang et al., 2014]
Po-Sen Huang, Haim Avron, Tara N Sainath, Vikas Sindhwani, and Bhuvana
Kernel methods match deep neural networks on timit.
In 2014 IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP), pages 205–209. IEEE, 2014.

[Koltchinskii, 2011]
Vladimir Koltchinskii.
Oracle Inequalities in Empirical Risk Minimization and Sparse
Recovery Problems: Ecole d’Eté de Probabilités de Saint-Flour
XXXVIII-2008, volume 2033.
Springer Science & Business Media, 2011.

[Li et al., 2019]
Zhu Li, Jean-Francois Ton, Dino Oglic, and Dino Sejdinovic.
Towards a unified analysis of random fourier features.
In International Conference on Machine Learning, pages
3905–3914. PMLR, 2019.

[Li et al., 2021]
Zhu Li, Jean-Francois Ton, Dino Oglic, and Dino Sejdinovic.
Towards a unified analysis of random fourier features.
Journal of Machine Learning Research, 220
(108):0 1–51, 2021.

[Liu et al., 2020]
Fanghui Liu, Xiaolin Huang, Yudong Chen, and Johan AK Suykens.
Random features for kernel approximation: A survey on algorithms,
theory, and beyond.
arXiv preprint arXiv:2004.11154, 2020.

[Mahoney and Drineas, 2009]
Michael W Mahoney and Petros Drineas.
CUR matrix decompositions for improved data analysis.
Proceedings of the National Academy of Sciences, 1060
(3):0 697–702, 2009.

[Mendelson, 2002]
Shahar Mendelson.
Improving the sample complexity using global data.
IEEE transactions on Information Theory, 480
(7):0 1977–1991, 2002.

[Nyström, 1930]
Evert J. Nyström.
Über die praktische Auflösung von Integralgleichungen
mit Anwendungen auf Randwertaufgaben.
Acta Mathematica, 1930.

[Rahimi and Recht, 2007]
Ali Rahimi and Benjamin Recht.
Random features for large-scale kernel machines.
In Advances in neural information processing systems, pages
1177–1184, 2007.

[Rahimi and Recht, 2009]
Ali Rahimi and Benjamin Recht.
Weighted sums of random kitchen sinks: Replacing minimization with
randomization in learning.
In Advances in neural information processing systems, pages
1313–1320, 2009.

[Rudi and Rosasco, 2017]
Alessandro Rudi and Lorenzo Rosasco.
Generalization properties of learning with random features.
In Advances in Neural Information Processing Systems, pages
3218–3228, 2017.

[Rudi et al., 2017]
Alessandro Rudi, Luigi Carratino, and Lorenzo Rosasco.
Falkon: An optimal large scale kernel method.
In Advances in Neural Information Processing Systems, pages
3891–3901, 2017.

[Rudin, 2017]
Walter Rudin.
Fourier analysis on groups.
Courier Dover Publications, 2017.

[Schölkopf and Smola, 2001]
Bernhard Schölkopf and Alexander J. Smola.
Learning with kernels: Support vector machines,
regularization, optimization, and beyond.
MIT Press, 2001.

[Schölkopf et al., 2004]
Bernhard Schölkopf, Koji Tsuda, and Jean-Philippe Vert.
Kernel methods in computational biology.
MIT press, 2004.

[Shalev-Shwartz and Ben-David, 2014]
Shai Shalev-Shwartz and Shai Ben-David.
Understanding machine learning: From theory to algorithms.
Cambridge university press, 2014.

[Shalev-Shwartz et al., 2011]
Shai Shalev-Shwartz, Yoram Singer, Nathan Srebro, and Andrew Cotter.
Pegasos: Primal estimated sub-gradient solver for svm.
Mathematical programming, 1270 (1):0 3–30,

[Smale and Zhou, 2003]
Steve Smale and Ding-Xuan Zhou.
Estimating the approximation error in learning theory.
Analysis and Applications, 10 (01):0 17–41,

[Smola and Schökopf, 2000]
Alex J Smola and Bernhard Schökopf.
Sparse greedy matrix approximation for machine learning.
In Proceedings of the Seventeenth International Conference on
Machine Learning, pages 911–918, 2000.

[Smola and Schölkopf, 2000]
Alexander J. Smola and Bernhard Schölkopf.
Sparse greedy matrix approximation for machine learning.
In Proceedings of the 17th International Conference on Machine
Learning, 2000.

[Sriperumbudur and Szabó, 2015]
Bharath Sriperumbudur and Zoltán Szabó.
Optimal rates for random Fourier features.
In Advances in Neural Information Processing Systems, pages
1144–1152, 2015.

[Steinwart and Christmann, 2008]
Ingo Steinwart and Andreas Christmann.
Support vector machines.
Springer Science & Business Media, 2008.

[Sun et al., 2018]
Yitong Sun, Anna Gilbert, and Ambuj Tewari.
But how does it work in theory? linear svm with random features.
In Advances in Neural Information Processing Systems, pages
3379–3388, 2018.

[Sutherland and Schneider, 2015]
Dougal J Sutherland and Jeff Schneider.
On the error of random Fourier features.
In Proceedings of the Thirty-First Conference on Uncertainty in
Artificial Intelligence, pages 862–871. AUAI Press, 2015.

[Tropp, 2015]
Joel A Tropp.
An introduction to matrix concentration inequalities.
Foundations and Trends® in Machine Learning,
80 (1-2):0 1–230, 2015.

[Tsybakov et al., 2004]
Alexander B Tsybakov et al.
Optimal aggregation of classifiers in statistical learning.
The Annals of Statistics, 320 (1):0 135–166,

[Williams and Seeger, 2001]
Christopher K. I. Williams and Matthias Seeger.
Using the Nyström method to speed up kernel machines.
In Advances in Neural Information Processing Systems 13,

[Williams and Seeger, 2001]
Christopher KI Williams and Matthias Seeger.
Using the nyström method to speed up kernel machines.
In Advances in neural information processing systems, pages
682–688, 2001b.

[Zhang et al., 2015]
Yuchen Zhang, John Duchi, and Martin Wainwright.
Divide and conquer kernel ridge regression: A distributed algorithm
with minimax optimal rates.
The Journal of Machine Learning Research, 160
(1):0 3299–3340, 2015.

§ NOTATION & DEFINITION

In the rest of the appendix, we will denote $P(x,y)$ to be the joint probability density function on $\mathcal{X} \times \mathcal{Y}$, where $\mathcal{X}$ is the instance space and $\mathcal{Y}$ is the label space. We will use $\|A\|$ to denote the vector norm if $A$ is a vector and operator norm if $A$ is a matrix or an operator. We will use $L_2(P_x)$ to denote the space of square-integrable functions with respect to $P_x$, the mariginal distribution of $P(x,y)$ with $\|\cdot\|_{L_2(P_x)}$ being the norm function and $\langle \cdot, \cdot \rangle_{L_2(P_x)}$ being the inner product. In addition, we use $\otimes$ to denote the $L_2(P_x)$ outer product. We will also use Tr$(A)$ to denote the trace of matrix $A$ or operator $A$.

Our analysis heavily depends on the notion of Rademacher complexity which we define below.

Suppose $\{x_i\}_{i=1}^n$ are i.i.d samples from $P_x$. Let $\mathcal{H}$ be a class of functions mapping $\mathcal{X}$ to $\mathbb{R}$. Then, the random variable known as the empirical Rademacher complexity is defined as
R̂_n(ℋ) = 𝔼_σ[sup_f∈ℋ|2/n∑_i=1^nσ_if(x_i)||x_1,⋯,x_n   ]  ,
where $\sigma_1,\cdots,\sigma_n$ are independent Rademacher random variables. The corresponding Rademacher complexity is then defined as
$$R_n(\mathcal{H}) = \mathbb{E}\Big[\hat{R}_n(\mathcal{H})\Big] \ ,$$

where the expectation is taken with respect to $P_x$.

Finally, for random Fourier feature vector $\phi_{q,s}(\cdot)$, we define the operator $\Psi_s : \mathbb{R}^s \rightarrow L_2(P_x)$ as
Ψ_sβ= ϕ_q,s^⊤β .

§ RANDOM FOURIER FEATURES APPROXIMATION OF $F_{\MATHCAL{H}}$

In this section, we consider using RFFs to approximate the optimal function $f_{\mathcal{H}}$ in $\mathcal{H}$. In other words, we first sample $v_1,\dots,v_s$ from some importance sampling distribution $q(v)$ and approximate $f_{\mathcal{H}}$ with $\phi_{q,s}(\cdot)^{\top}\beta$ for some $\beta \in \mathbb{R}^s$. To find the best $\beta$, we form the following optimization problem:
β̃_ℋ:= _β∈ℝ^s ϕ_q,s^⊤β- f_ℋ_L_2(P_x)^2 + λβ_2^2 .

We denote the $\tilde{f}_{\mathcal{H}}$ as the RFFs approximation of $f_{\mathcal{H}}$. Our next theorem indicates the approximation property of $\tilde{f}_{\mathcal{H}}$.

Assume A.$1$ and A.$3$ hold. Suppose we have a measurable function $\tilde{\tau}: \mathcal{V} \rightarrow \mathbb{R}$ such that $\tilde{\tau}(v) \geq \tau_{\lambda}(v)$ almost surely, denote $d_{\tilde{\tau}} = \int_{\mathcal{V}}\tilde{\tau}(v) dv$, and let $q(v) = \frac{\tilde{\tau}(v)}{d_{\tilde{\tau}}}$. Let $\delta \in (0,1)$, if we sample $v_1,\dots,v_s \sim q(v)$ such that
\[s \geq 12 d_{\tilde{\tau}}\log\frac{d(\lambda)}{\delta},\]
with probability over $1-\delta$
\[\left\| \tilde{f}_{\mathcal{H}}- f_{\mathcal{H}}\right\|_{L_2(P_x)} \leq 2R\lambda^{r}.\]

We approximate $f_{\mathcal{H}}$ with $\tilde{f}_q(\cdot) = \phi_{q,s}^{\top}\beta$ and formulate the optimization problem in Eq. (<ref>). According to Eq. (<ref>), we can further rewrite Eq. (<ref>) as
Ψ_sβ- f_ℋ_L_2(P_x)^2 + λβ_2^2  .

The solution of the above optimization problem can be computed as
\[\tilde{\beta}_{\mathcal{H}} = \Psi_s^*\left(\Psi_s\Psi_s^* + \lambda I \right)^{-1}f_{\mathcal{H}} \,\]

where $\Psi^*$ denotes the adjoint operator of $\Psi$. Hence, we have
f̃_ℋ- f_ℋ    =    Ψ_sβ̃_ℋ - f_ℋ   ,

=      Ψ_sΨ_s^*(Ψ_sΨ_s^* + λI )^-1f_ℋ -f_ℋ  ,

=     L_s(L_s + λI )^-1f_ℋ -f_ℋ  ,

=     λ(L_s + λI )^-1f_ℋ  ,

=    λ(L + λI )^-1/2(I + (L + λI )^-1/2(L_s - L) (L + λI )^-1/2)^-1 (L + λI )^-1/2 f_ℋ  ,

=    λ(L + λI )^-1/2(I + (L + λI )^-1/2(L_s - L) (L + λI )^-1/2)^-1 (L + λI )^-1/2 L^rg  ,

where for last step we used assumption A.$3$.

In the mean time, we have
f̃_ℋ    ≤     β̃_ℋ_2^2 = ⟨f_ℋ,(L_s + λI )^-1L_s(L_s + λI )^-1f_ℋ⟩_L_2(P_x)  ,

≤    ⟨f_ℋ,(L_s + λI )^-1 f_ℋ⟩_L_2(P_x) ,

=    ⟨f_ℋ,(L + λI )^-1/2(I + (L + λI )^-1/2(L_s - L) (L + λI )^-1/2)^-1 (L + λI )^-1/2 f_ℋ⟩_L_2(P_x)

We now utilize Lemma <ref> to lower bound $\left(L + \lambda I \right)^{-\frac{1}{2}}(L_s - L) \left(L + \lambda I \right)^{-\frac{1}{2}}$. Specifically, let $\epsilon = 1/2$ in Lemma <ref>, we obtain that when
\[s \geq 12 d_{\tilde{\tau}}\log\frac{d(\lambda)}{\delta},\]
we have with probability over $1-\delta$ such that
\[\left\|\left(L + \lambda I \right)^{-\frac{1}{2}}(L_s - L) \left(L + \lambda I \right)^{-\frac{1}{2}}\right\| \preceq \frac{1}{2}I.\]

f̃_ℋ- f_ℋ_L_2(P_x)   ≤    (1-1/2)^-1 λ(L + λI )^-1L^rg_L_2(P_x)  ,

=    2λ(L + λI )^-(1-r) (L + λI )^-rL^rg_L_2(P_x)  ,

≤    2 λ^rg_L_2(P_x)  ,

=    2Rλ^r .

Finally, recall $r \in [1/2,1]$, we have
\[\|\tilde{\beta}_{\mathcal{H}}\|_2^2 \leq 2 \langle f_{\mathcal{H}},\left(L + \lambda I \right)^{-1} f_{\mathcal{H}}\rangle_{L_2(P_x)}= 2 \langle g,\left(L + \lambda I \right)^{-1} L^{2r}g\rangle_{L_2(P_x)} \leq 2R^2.\]

§ APPROXIMATION ERROR OF INTEGRAL OPERATOR

The following lemma characterizes the approximation error of $L$ using $L_s$.

Let $L$ and $L_s$ be defined as in Eq. (<ref>) and Eq.(<ref>) respectively and recall we denote the leverage score as $\tau_{\lambda}(v)$. Suppose we have a measurable function $\tilde{\tau}: \mathcal{V} \rightarrow \mathbb{R}$ such that $\tilde{\tau}(v) \geq \tau_{\lambda}(v)$ almost surely, denote $d_{\tilde{\tau}} = \int_{\mathcal{V}}\tilde{\tau}(v) dv$, and let $q(v) = \frac{\tilde{\tau}(v)}{d_{\tilde{\tau}}}$. Assume A.$1$ holds and we draw random features $v_1\dots,v_s \sim q(v)$, for some constant $\epsilon > 0$ and $\delta \in (0,1)$, if the number of features satisfies
\[s \geq \left(\frac{2}{\epsilon^2} + \frac{4}{3\epsilon} \right) d_{\tilde{\tau}}\log\frac{d(\lambda)}{\delta},\]
then with probability over $1-\delta$, we have
\[\left\|\left(L+ \lambda I\right)^{-\frac{1}{2}}(L_s - L)\left(L+ \lambda I\right)^{-\frac{1}{2}} \right\| \leq \epsilon.\]

We employ the Bernstein inequality to prove the lemma. To this end, we first define
\[\mathbf{R}_i = \frac{p(v_i)}{q(v_i)}\left(L+ \lambda I\right)^{-\frac{1}{2}}\psi(v_i,\cdot)\otimes\psi(v_i,\cdot) \left(L+ \lambda I\right)^{-\frac{1}{2}},\]
then by Eq. (<ref>) we have
\[\left(L+ \lambda I\right)^{-\frac{1}{2}} L_s\left(L+ \lambda I\right)^{-\frac{1}{2}} = \frac{1}{s}\sum_{i=1}^s \mathbf{R}_i.\]
We then immediately have $\mathbb{E}\left(\mathbf{R}_i\right) = \left(L+ \lambda I\right)^{-\frac{1}{2}} L\left(L+ \lambda I\right)^{-\frac{1}{2}}$.

Furthermore, we have
𝐑_i    =    p(v_i)/q(v_i)(L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2   ,

=    p(v_i)/q(v_i) ⟨ψ(v_i,·), (L+ λI)^-1ψ(v_i,·) ⟩_L_2(P_x)  ,

=    τ_λ(v_i)/q(v_i) =  τ_λ(v_i)/τ̃(v_i)d_τ̃  ,

≤    d_τ̃ .

In addition, $\mathbf{R}_i$ is a self-adjoint operator, we have
𝐑_i𝐑_i^*    =    𝐑_i^*𝐑_i  ,

=    p(v_i)^2/q(v_i)^2(L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2  ,

=    p(v_i)/q(v_i) ⟨ψ(v_i,·), (L+ λI)^-1ψ(v_i,·) ⟩_L_2(P_x) p(v_i)/q(v_i) (L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2  ,

=    τ_λ(v_i)/q(v_i)p(v_i)/q(v_i)(L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2  ,

≤    d_τ̃p(v_i)/q(v_i)(L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2  .

As a result, we have
𝔼( 𝐑_i 𝐑_i^*)    ≼    d_τ̃ 𝔼( p(v_i)/q(v_i)(L+ λI)^-1/2ψ(v_i,·)⊗ψ(v_i,·) (L+ λI)^-1/2)  ,

=   d_τ̃ (L+ λI)^-1/2 L(L+ λI)^-1/2 := 𝐌 .

It is easy to see $m = \|\mathbf{M}\| \leq d_{\tilde{\tau}}$. We now let $\bar{\mathbf{R}}_s = \frac{1}{s}\sum_{i=1}^s \mathbf{R}_i$ and $d = \frac{\textnormal{Tr}(\mathbf{M})}{d_{\tilde{\tau}}} = d(\lambda)$, where we recall $d(\lambda) = \textnormal{Tr}\left(L(L+\lambda I)^{-1}\right)$. Applying the Bernstein inequality in Lemma <ref>, we have
\[P\left(\left\| \bar{\mathbf{R}}_s - \mathbf{R}\right\| \geq \epsilon \right)\leq 4d(\lambda)\exp\left(\frac{-s\epsilon^2}{d_{\tilde{\tau}}+2d_{\tilde{\tau}}}\epsilon/3 \right).\]

By letting $4d(\lambda)\exp\left(\frac{-s\epsilon^2}{d_{\tilde{\tau}}+2d_{\tilde{\tau}}}\epsilon/3 \right) \leq \delta$ and rearrange the terms, we prove the lemma.

§ PROOF OF THEOREM <REF>

To this end, we now consider the estimator $\tilde{f}^{\lambda} \in \tilde{\mathcal{H}}$ obtained as $\tilde{f}^{\lambda}: = \argmin \mathbb{E}_n(l_{\tilde{f}})$ subject to $\|\beta\|_2^2 \leq 2R^2$ and recall $\tilde{f}_{\mathcal{H}} \in \tilde{\mathcal{H}}$, we have the following decomposition
𝔼(l_f̃^λ)    =    𝔼(l_f̃^λ)- 𝔼_n(l_f̃^λ)+𝔼_n(l_f̃^λ)-𝔼_n(l_f̃_ℋ)+𝔼_n(l_f̃_ℋ)-𝔼(l_f̃_ℋ)+𝔼(l_f̃_ℋ)-𝔼(l_f_ℋ)+𝔼(l_f_ℋ)  ,

≤    𝔼(l_f̃^λ)- 𝔼_n(l_f̃^λ) +𝔼_n(l_f̃_ℋ)-𝔼(l_f̃_ℋ)+𝔼(l_f̃_ℋ)-𝔼(l_f_ℋ)+𝔼(l_f_ℋ)  ,

≤    O(1/√(n))  +𝔼(l_f̃_ℋ)-𝔼(l_f_ℋ)+𝔼(l_f_ℋ)  ,

≤    O(1/√(n))  + M f̃_ℋ-f_ℋ _L_2(P_x)+𝔼(l_f_ℋ)  ,

≤    2MRλ^r + O(1/√(n))  +  𝔼(l_f_ℋ) .

where for the second step, we use the fact that $\tilde{f}^{\lambda}$ is the minimizer of the empirical risk. For the third step, we use standard results for Rademacher complexities of $L_2$-balls <cit.> applied to $\mathbb{E}(l_{\tilde{f}^{\lambda}})- \mathbb{E}_n(l_{\tilde{f}^{\lambda}})$ and $\mathbb{E}_n(l_{\tilde{f}_{\mathcal{H}}})-\mathbb{E}(l_{\tilde{f}_{\mathcal{H}}})$. For the fourth step, we utilize the Lipschitz continuity of the loss function and the last step is due to Theorem <ref>.

§ PROOF OF THEOREM <REF>

To prove Theorem <ref>, we rely on the notion of local Rademacher complexity introduced in Lemma <ref>. In order to do that, we need two steps. The first step is to find a proper sub-root function $\hat{\varphi}_n(r)$. The second step is to find the fixed point of $\hat{\varphi}_n(r)$. Hence, the following is devoted to solving these two problems.

We first recall that $\tilde{\mathcal{H}}$ is an RKHS and hence is convex. Therefore, it is a star hull around every of its interior point (i.e., $\operatorname{star}\left(\tilde{\mathcal{H}}, f_{0}\right)=\left\{f_{0}+\alpha\left(f-f_{0}\right) \mid f \in \tilde{\mathcal{H}} \wedge \alpha \in[0,1]\right\} = \tilde{\mathcal{H}} , \forall f_0 \in \tilde{\mathcal{H}}$). In addition, let $f_{\tilde{\mathcal{H}}}$ be the estimator minimizing the expected risk in RKHS $\tilde{\mathcal{H}}$, subject to $\|f\|_{\tilde{\mathcal{H}}} \leq 2R^2$, Assumption A.$4$ then implies that $\mathbb{E}\left(f-f_{\tilde{\mathcal{H}}}\right)^2 \leq B \mathbb{E}\left(l_f - l_{f_{\tilde{\mathcal{H}}}}\right)$. We also recall $\tilde{f}^{\lambda} := \argmin_{f \in \tilde{\mathcal{H}}} \mathbb{E}_n(l_f)$ subject to $\|f\|_{\tilde{\mathcal{H}}} \leq 2R^2$.

Let $\psi(r)$ be a sub-root function and satisfy

\[\psi(r) \geq B L \mathbb{E} R_{n}\left\{f \in \tilde{\mathcal{H}}: L^{2} \mathbb{E}\left(f-f_{\tilde{\mathcal{H}}}\right)^{2} \leq r\right\}.\]

We apply Lemma <ref> to the class $l_{f}-l_{f_{\tilde{\mathcal{H}}}}$ with $T(f)=L^{2} \mathbb{E}\left(f-f_{\tilde{\mathcal{H}}}\right)^{2}$. By the contraction property of Rademacher complexity <cit.> and the symmetry of the Rademacher variables, $L \mathbb{E} R_{n}\left\{f: L^{2} \mathbb{E}\left(f-f_{\tilde{\mathcal{H}}}\right)^{2} \leq r\right\} \geq \mathbb{E} R_{n}\left\{l_{f}-l_{f_{\tilde{\mathcal{H}}}}: L^{2} \mathbb{E}\left(f-f_{\tilde{\mathcal{H}}}\right)^{2} \leq r\right\}$. Now by Lemma <ref> and noticing that $\mathbb{E}_{n}\left(l_{\tilde{f}^{\lambda} }-l_{f_{\tilde{\mathcal{H}}}}\right) \leq 0$, we have for any $\delta \in (0,1)$ and any $r \geq \psi(r)$, with probability at least $1-\delta$

\[\mathbb{E}\left(l_{\tilde{f}^{\lambda} }-l_{f_{\tilde{\mathcal{H}}}}\right) \leq  c_1 r+\frac{c_2}{n}\log \frac{1}{\delta}.\]

We now define
\[\hat{\psi}(r)=c_{3} \hat{R}_{n}\left\{f \in \tilde{\mathcal{H}}: \mathbb{E}_{n}(f-\hat{f})^{2} \leq c_{4} r\right\}+\frac{c_{5}}{n}.\]

By <cit.>, we have $\psi(r) \leq \hat{\psi}(r)$ and $r^* \leq \hat{r}^*$ where $r^*$ and $\hat{r}^*$ are the fixed points of $\psi(r)$ and $\hat{\psi}$ respectively. Therefore, with probability at least $1-2\delta$,
𝔼(l_f̃^λ -l_f_ℋ̃) ≤c_6 r̂^*+c_7/nlog1/δ  .

We are now left to compute the fixed point $\hat{r}^*$. To this end, we notice

\begin{align}
\begin{aligned}
&\hat{R}_n\{ f - f_{\tilde{\mathcal{H}}}, f\in \tilde{\mathcal{H}} \mid \mathbb{E}_n[f-f_{\tilde{\mathcal{H}}}]^2 \leq r \} \leq  &\\
&\hat{R}_n\{ f - g, f, g \in \tilde{\mathcal{H}} \mid \mathbb{E}_n[f-g]^2 \leq r \} =  &\\
&2\hat{R}_n\{ f , f\in \tilde{\mathcal{H}} \mid \mathbb{E}_n[f]^2 \leq r/4 \}  &
\end{aligned}\nonumber%\label{apn:rade_bound}
\end{align}

We now let $\tilde{\mu}_1,\dots,\tilde{\mu}_s$ to be the eigenvalues of the normalized Gram-matrix $1/n\tilde{\mathbf{K}}$, by Lemma <ref>, we have
\[2\hat{R}_n\{ f , f\in \tilde{\mathcal{H}} \mid \mathbb{E}_n[f]^2 \leq r/4 \} \leq 2 \left(\frac{2}{n}\sum_{i=1}^n\min\{r/4,\tilde{\mu}_i\}\right)^{1/2}.\]

For sub-root function $\hat{\psi}(r) = 2 \left(\frac{2}{n}\sum_{i=1}^n\min\{r/4,\tilde{\mu}_i\}\right)^{1/2} + \frac{c_5}{n}$, we first notice that adding some constant $c$ will only increase the fixed point by at most $2c$. Therefore, it suffice to compute the fixed point of

\[\hat{r} \leq  2M \left(\frac{2}{n}\sum_{i=1}^n\min\{\hat{r}/4,\tilde{\mu}_i\}\right)^{1/2}.\]

To this end, we have

(r̂/2M)^2     ≤    2/n ∑_i=1^n min{r̂/4, μ̃_i}  ,

=   2/n min_S ⊆{1, …, n}(∑_i ∈S r̂/4+∑_i ∉S μ̃_i) ,

=   2/n min_0 ≤h ≤n( h r̂/4+∑_i>h μ̃_i)  .

Finally, we solve the quadratic inequality for each value of $h$ implies

\[\hat{r} \leq c \min _{0 \leq h \leq n}\left(\frac{h}{n}+\sqrt{\frac{1}{n} \sum_{i>h} \tilde{\mu}_{i}}\right).\]

Recall the estimator $\tilde{f}^{\lambda} \in \tilde{\mathcal{H}}$ is obtained as $\tilde{f}^{\lambda}: = \argmin \mathbb{E}_n(l_{\tilde{f}})$ subject to $\|\beta\|_2^2 \leq 2R^2$ and recall $f_{\tilde{\mathcal{H}}} \in \tilde{\mathcal{H}}$.
By Eq. <ref>, we have
𝔼(l_f̃^λ-l_f_ℋ̃)    ≤    c_6r̂^* + c_7/nlog1/δ  ,

Equivalently, we can rewrite the above equation as
𝔼(l_f̃^λ)-𝔼(l_f_ℋ)   ≤    𝔼(l_f_ℋ̃) - 𝔼(l_f_ℋ) + c_1r̂^* + c_2/nlog1/δ  ,

≤    𝔼(l_f̃_ℋ) - 𝔼(l_f_ℋ) + c_1r̂^* + c_2/nlog1/δ  ,

≤    M f̃_ℋ - f_ℋ _L_2(P_x)+ c_1r̂^* + c_2/nlog1/δ  ,

≤    2MRλ^r + c_1r̂^* + c_2/nlog1/δ  ,

where the second step is because $\mathbb{E}\left( l_{f_{\tilde{\mathcal{H}}}}\right) \leq \mathbb{E}\left(l_{\tilde{f}_{\mathcal{H}}}\right)$.

§ MATRIX BERNSTEIN INEQUALITY

The next lemma is the matrix Bernstein inequality, cited from <cit.> which is a restatement of Corollary 7.3.3 in Tropp, 2015 with some fix in the typos.
Let $\mathbf{R}$ be a fixed $d_1 \times d_2$ matrix over the set of complex/real numbers. Suppose that $\{\mathbf{R}_1,\cdots,\mathbf{R}_n\}$ is an independent and identically distributed sample of $d_1 \times d_2$ matrices such that
\[\mathbb{E}[\mathbf{R}_i] = \mathbf{R} \qquad \text{and} \qquad \|\mathbf{R}_i\|_2 \leq L\ ,\]

where $L>0$ is a constant independent of the sample.
Furthermore, let $\mathbf{M}_1, \mathbf{M}_2$ be semidefinite upper bounds for the matrix-valued variances

\begin{align*}
\begin{aligned}
& \mathrm{Var}_1[\mathbf{R}_i] \preceq \mathbb{E}[\mathbf{R}_i\mathbf{R}_i^{T}] \preceq \mathbf{M}_1 & \\
& \mathrm{Var}_2[\mathbf{R}_i] \preceq \mathbb{E}[\mathbf{R}_i^{T}\mathbf{R}_i]\preceq \mathbf{M}_2\ . &
\end{aligned}
\end{align*}

Let $m = \max(\|\mathbf{M}_1\|_2,\|\mathbf{M}_2\|_2)$ and $d =\frac{\text{Tr}(\mathbf{M}_1)+ \text{Tr}(\mathbf{M}_2)}{m}.$
Then, for $\epsilon \geq \sqrt{m/n}+2L/3n$, we can bound
\[\bar{\mathbf{R}}_n = \frac{1}{n}\sum_{i=1}^{n}\mathbf{R}_i\]
around its mean using the concentration inequality
\[P(\|\bar{\mathbf{R}}_n - \mathbf{R}\|_2 \geq \epsilon) \leq 4d\exp\Bigg(\frac{-n\epsilon^2/2}{m+2L\epsilon/3}\Bigg)\ .\]

§ LOCAL RADEMACHER COMPLEXITIES

In the refined case, our analysis relies on the local Rademacher complexities technique developed by Bartlett et al., 2005. As shall be seen later, local Rademacher complexity is often linked with the so-called sub-root function which we define below.

Let $\varphi: [0,\infty) \rightarrow [0,\infty)$ be a function. Then, $\varphi(r)$ is called a sub-root function if, for all $r>0$, $\varphi(r)$ is non-decreasing and $\frac{\varphi(r)}{r}$ is non-increasing.

A sub-root function has the following property.

If $\varphi(r)$ is a sub-root function, then $\varphi(r) = r$ has a unique positive solution $r^*$. In addition, we have that $r \geq \varphi(r) $ if and only if $r \geq r^*$.

After we introduce the Rademacher complexity and the sub-root function, we are now ready to state the theorem on local Rademacher complexity we need in the refined analysis.

Let $\mathcal{F}$ be a class of functions with bounded ranges and assume that there are some functional $T: \mathcal{F} \rightarrow \mathbb{R}^{+}$ and some constant $B$ such that for every $f \in \mathcal{F}, \textnormal{Var}[f] \leq T(f) \leq B P f .$ Let $\psi$ be a sub-root function and $r^{*}$ be the fixed point.
Assume that $\psi$ satisfies, for any $r \geq r^{*}$

\[\psi(r) \geq B R_{n}\{f \in \mathcal{F}: T(f) \leq r\}.\]

Then, for any $f \in \mathcal{F}$, some constant $D>1$ and every $\delta \in (0,1)$, with probability at least $1-\delta$

\[\mathbb{E}(f) \leq \frac{D}{D-1} \mathbb{E}_n(f) +\frac{c_{1} D}{B} r^{*}+\frac{c_2}{n} \log \frac{1}{\delta},\]
where $c_1,c_2$ are some universal constants.

In addition to the above theorem, one also need to quantitatively characterize the local Rademacher complexity in order to obtain a tight upper bound. Fortunately in the kernel case, we are able to do so as illustrate in the following lemma.
Let $k$ be a positive definite kernel function with reproducing kernel Hilbert space $\mathcal{H}$ and let $\hat{\lambda}_1 \geq \cdots \geq \hat{\lambda}_n$ be the eigenvalues of the normalized Gram-matrix $(1/n)\mathbf{K}$. Then, for all $r >0$ and $f \in \mathcal{H}$,
\[\hat{R}_n\{f\in \mathcal{H} \mid \mathbb{E}_n[f^2] \leq r\} \leq \left(\frac{2}{n}\sum_{i=1}^n\min\{r,\hat{\lambda}_i\}\right)^{1/2}.\]