Year: 2022

Abstract:
In this article, we investigate the spectral behavior of random features kernel matrices of the type K = Ew[σ(wxi)σ(wxj)]i,j=1, with nonlinear function σ(·), data x1, . . . ,xn ∈ R, and random projection vector w ∈ R having i.i.d. entries. In a high-dimensional setting where the number of data n and their dimension p are both large and comparable, we show, under a Gaussian mixture model for the data, that the eigenspectrum of K is independent of the distribution of the i.i.d. (zero-mean and unit-variance) entries of w, and only depends on σ(·) via its (generalized) Gaussian moments Ez∼N (0,1)[σ(z)] and Ez∼N (0,1)[σ(z)]. As a result, for any kernel matrix K of the form above, we propose a novel random features technique, called Ternary Random Feature (TRF), that (i) asymptotically yields the same limiting kernel as the original K in a spectral sense and (ii) can be computed and stored much more efficiently, by wisely tuning (in a data-dependent manner) the function σ and the random vector w, both taking values in {−1, 0, 1}. The computation of the proposed random features requires no multiplication, and a factor of b times less bits for storage compared to classical random features such as random Fourier features, with b the number of bits to store full precision values. Besides, it appears in our experiments on real data that the substantial gains in computation and storage are accompanied with somewhat improved performances compared to state-of-the-art random features compression/quantization methods.

heading: 1 INTRODUCTION
text: Kernel methods are among the most powerful machine learning approaches with a wide range of successful applications (Schölkopf & Smola, 2018) which, however, suffer from scalability issues in large-scale problems, due to their high space and time complexities (with respect to the number of data n). To address this key limitation, a myriad of random features based kernel approximation techniques have been proposed (Rahimi & Recht, 2008; Liu et al., 2021a): random features methods randomly project the data to obtain low-dimensional nonlinear representations that approximate the original kernel features. This allows practitioners to apply them, with a large saving in both time and space, to various kernel-based downstream tasks such as kernel spectral clustering (Von Luxburg, 2007), kernel principal component analysis (Schölkopf et al., 1997), kernel canonical correlation analysis (Lai & Fyfe, 2000), kernel ridge regression (Vovk, 2013), to name a few. A wide variety of these kernels can be written, for data points xi,xj ∈ Rp, in the form
κ(xi,xj) = Ew [ σ ( wTxi ) σ ( wTxj )] (1)
with w ∈ Rp having i.i.d. entries, which can be “well approximated” by a sample mean 1 m ∑m t=1 σ ( wTt xi ) σ ( wTt xj ) over m independent random features for m sufficiently large. For instance, taking σ(x) = [cos(x), sin(x)] and w with i.i.d. standard Gaussian entries, one obtains the popular Random Fourier Features (RFFs) that approximate the Gaussian kernel (and the
Laplacian kernel for Cauchy distributed w with the same choice of σ) (Rahimi & Recht, 2008); for σ(x) = max(x, 0), one approximates the first order Arc-cosine kernel; and the zeroth order Arc-cosine kernel (Cho, 2012) with σ(x) = (1 + sign(x))/2, etc.
As shall be seen subsequently, (random) neural networks are, to a large extent, connected to kernel matrices of the form (1). More specifically, the classification or regression performance at the output of random neural networks are functionals of random matrices that fall into the wide class of kernel random matrices. Perhaps more surprisingly, this connection still exists for deep neural networks which are (i) randomly initialized and (ii) trained with gradient descent, as testified by the recent works on neural tangent kernels (Jacot et al., 2018), by considering the “infinitely many neurons” limit, that is, the limit where the network widths of all layers go to infinity simultaneously. This close connection between neural networks and kernels has triggered a renewed interest for the theoretical investigation of deep neural networks from various perspectives, including optimization (Du et al., 2019; Chizat et al., 2019), generalization (Allen-Zhu et al., 2018; Arora et al., 2019; Bietti & Mairal, 2019), and learning dynamics (Lee et al., 2019; Advani et al., 2020; Liao & Couillet, 2018a). These works shed new light on the theoretical understanding of deep neural network models and specifically demonstrate the significance of studying networks with random weights and their associated kernels to assess the mechanisms underlying more elaborate deep networks.
In this article, we consider the random feature kernel of the type (1), which can also be seen as the limiting kernel of a single-hidden-layer neural network with a random first layer. By assuming a high-dimensional Gaussian Mixture Model (GMM) for the data {xi}ni=1 with xi ∈ Rp, we show that the centered kernel matrix1
K , P{κ(xi,xj)}ni,j=1P, P , In − 1
n 1n1
T n, (2)
is asymptotically (as n, p→∞ with p/n→ c ∈ (0,∞)) equivalent, in a spectral sense, to another random kernel matrix K̃ which depends on the GMM data statistics and the generalized Gaussian moments E[σ′(z)],E[σ′′(z)] of the activation function σ(·), but is independent of the specific law of the i.i.d. entries of the random vector w, as long as they are normalized to have zero mean and unit variance. As such, one can design novel random features schemes with limiting kernels asymptotically equivalent to the original K. For instance, define
κter(xi,xj) , Ewter [ σter ( xTi w ter ) σter ( xTj w ter )]
(3)
with wter ∈ Rp having i.i.d. entries taking value wteri = 0 (with probability ) and value wteri ∈ { −(1− )− 12 , (1− )− 12 } (4)
each with probability 1/2− /2, where ∈ [0, 1) represents the level of sparsity of w, and
σter(t) = −1 · δt<s− + 1 · δt>s+ (5)
for some thresholds s− < s+ chosen to match the generalized Gaussian moments E[σ′(z)],E[σ′′(z)] of any σ function (e.g., ReLU, cos, sin) widely used in random features or neural network contexts. The proposed Ternary Random Features (TRFs, with limiting kernel matrices defined in (3) asymptotically “matching” any random features kernel matrices in a spectral sense) have the computational advantage of being sparse and not requiring multiplications but only additions, as well as the storage advantage of being composed of only a finite set of words, e.g., {−1, 0, 1} for = 0. Given the urgent need for environmentally-friendly, but still efficient, neural networks such as binary neural networks (Hubara et al., 2016; Lin et al., 2015; Zhu et al., 2016; Qin et al., 2020; Hubara et al., 2016), pruned neural networks (Liu et al., 2015; Han et al., 2015a;b), weights-quantized neural networks (Gupta et al., 2015; Gong et al., 2014), we hope that our analysis will open a new door to a random matrix-improved framework of computationally efficient methods for machine learning and neural network models more generally.
1Left- and right-multiplying the kernel matrices by P is equivalent to centering the data in the kernel feature space, which is a common practice in kernel learning and plays a crucial role in multidimensional scaling (Joseph & Myron, 1978) and kernel PCA (Schlkopf et al., 1998). In the remainder of this paper, whenever we use kernel matrices, they are considered to have been centered in this way.

heading: 1.1 CONTRIBUTIONS
text: Our main results are summarized as follows.
1. By considering a high-dimensional Gaussian mixture model for the data, we show (Theorem 1) that for K defined in (2), ‖K− K̃‖ → 0 as n, p→∞, where K̃ is a random matrix independent of the law of w, and depends on the nonlinear σ(·) only via its generalized Gaussian moments E[σ′(z)] and E[σ′′(z)] for z ∼ N (0, 1).
2. We exploit this result to propose a computationally efficient random features approach, called Ternary Random Features (TRFs), with asymptotically the same limiting kernel as any random features-type kernel matrices K of the form (2), while inducing no multiplication and b times less memory storage for its computation, with b the number of bits to store full precision values, e.g., b = 32 in the case of a single-precision floating-point format.
3. We provide empirical evidence on various random-features based algorithms, showing the computational and storage advantages of the proposed TRF method, while achieving competitive performances compared to state-of-the-art random features techniques. As a first telling example, in Figure 1 on the Cov-Type dataset from the UCI ML repository,TRFs achieve similar logistic regression performance as the LP-RFF method proposed by Zhang et al. (2019), with 8 times less memory, and 32 or 16 times less memory than the Nyström approximation of the Gaussian kernel (using full precision of 32 or 16 bits) (Williams & Seeger, 2001); see the shift in the x-axis (memory) of Figure 1.

heading: 1.2 RELATED WORK
text: Random features kernel and random neural networks. Random features methods were first proposed to relieve the computational and storage burden of kernel methods in large-scale problems when the number of training samples n is large (Schölkopf & Smola, 2018; Rahimi & Recht, 2008; Liu et al., 2021a). For instance, Random Fourier Features can be used to approximate the popular Gaussian kernel, when the number of random features m is sufficiently large (Rahimi & Recht, 2008). Since (deep) modern neural networks are routinely trained with random initialization, random features method is also considered a stylish model to analyze neural networks (Neal, 1996; Williams, 1997; Novak et al., 2018; Matthews et al., 2018; Lee et al., 2017; Louart et al., 2018). In particular, by focusing on the large n,m regime, the analysis of such models led to the so-called double descent theory (Advani et al., 2020; Mei & Montanari, 2019; Liao et al., 2020) for neural nets.
Computationally efficient random features methods. In an effort to further reduce the computation and storage costs of random features models, various quantization and binarization methods were proposed (Goemans & Williamson, 1995; Charikar, 2002; Li & Slawski, 2017; Li & Li, 2019; 2021;
Agrawal et al., 2019; Zhang et al., 2019; Liao et al., 2021; Couillet et al., 2021). More precisely, Agrawal et al. (2019) combined RFFs with a data-dependent feature selection approach to reduce the computational cost, while preserving the statistical guarantees of (using the original set of) RFFs. Zhang et al. (2019) proposed a low-precision approximation of RFFs to significantly reduce the storage while generalizing as well as full-precision RFFs. Li & Li (2021) designed quantization schemes of RFFs for arbitrary choice of the Gaussian kernel parameter. Our work improves these previous efforts by (i) considering a broader family of random features kernels beyond RFFs and by (ii) proposing the TRF approach that is both sparse and quantized, while asymptotically yielding the same limiting kernel spectral structure, and thus algorithmic performances (Cortes et al.).
Random matrix theory and neural networks. Random matrix theory (RMT), as a powerful and flexible tool to investigate the (asymptotic) behavior of large-scale systems, is recently gaining popularity in the analysis of (deep) neural networks (Martin & Mahoney, 2018; 2019). In this respect, Pennington & Worah (2017) derived the eigenvalue distribution of the Conjugate Kernel (CK) in a single-hidden-layer random neural network model. This result was then generalized to a broader class of data distributions (Louart et al., 2018) and to a multi-layer scenario (Benigni & Péché, 2019; Pastur, 2020). Fan & Wang (2020) went beyond the general i.i.d. assumption (on the entries of data vectors) and studied the spectral properties of the CK and neural tangent kernel for data that are approximately “pairwise orthogonal.” Our work improves (Fan & Wang, 2020) by studying the random features kernel for more structured GMM data, and is thus more adapted to machine learning applications such as classification. As far as the study of random features kernels under GMM data is concerned, the closest work to ours is (Liao & Couillet, 2018b) where the kernel matrix K defined in (2) is studied for GMM data, but only for a few specific activation functions and Gaussian w, see Footnote 5 below for a detailed account of the technical differences between this work and (Liao & Couillet, 2018b). Here, we provide a universal result with respect to the (much broader class of) activation functions and random w, and propose a computation and storage efficient random features technique well tuned to match the performances of any commonly used random features kernels.

heading: 1.3 NOTATIONS AND ORGANIZATION OF THE ARTICLE
text: In this article, we denote scalars by lowercase letters, vectors by bold lowercase, and matrices by bold uppercase. We denote the transpose operator by (·)T, we use ‖ · ‖ to denote the Euclidean norm for vectors and spectral/operator norm for matrices. For a random variable z, E[z] denotes the expectation of z. The notation δx∈A is the Kronecker delta taking value 1 when x ∈ A and 0 otherwise. Finally, 1p and Ip are respectively the vector of all one’s of dimension p and the identity matrix of dimension p× p. The remainder of the article is structured as follows. In Section 2, we describe the random features model under study along with our working assumptions. We then present our main technical results in Section 3 on the spectral characterization of random features kernel matrices K and its practical consequences, in particular the design of cost-efficient ternary random features (TRFs) leading asymptotically to the same kernel as any generic random features. We provide empirical evidence in Section 4 showing the computational and storage advantage along with competitive performance of TRFs compared to state-of-the-art random features approaches. Conclusion and perspective are placed in Section 5.

heading: 2 SYSTEM SETTINGS
text: Let W ∈ Rm×p be a random matrix having i.i.d. entries with zero mean and unit variance. The random features matrix Σ ∈ Rm×n of some data X ∈ Rp×n is defined as Σ , σ(WX) for some nonlinear activation function σ : R→ R applied entry-wise on WX. We denote G the associated random features Gram matrix
G , 1
m ΣTΣ =
1 m σ(WX)Tσ(WX) = 1 m m∑ t=1 σ ( XTwt ) σ ( wTt X ) ∈ Rn×n (6)
which is a sample mean of the expected kernel defined in (1). With P , In − 1n1n1 T n, we consider the expected and centered random features kernel matrix K defined in (2), which plays a fundamental
role in various random features kernel-based learning methods such as kernel ridge regression, logistic regression, support vector machines, principal component analysis, or spectral clustering.
Let x1, · · · ,xn ∈ Rp be n independent data vectors belonging to one of K distributional classes C1, · · · , CK , with class Ca having cardinality na. We assume that xi follows a Gaussian Mixture Model (GMM), that is, for xi ∈ Ca,
xi = µa/ √ p+ zi (7)
with zi ∼ N (0p,Ca/p) for some mean µa ∈ Rp and covariance Ca ∈ Rp×p associated to class Ca.
In high dimensions, under “non triviality assumptions”,2 the data vectors drawn from (7) can be shown to be neither very “close” nor very “far” from each other, irrespective of the class they belong to, see (Couillet et al., 2016). We place ourselves under such non-trivial conditions, by imposing, as in (Couillet et al., 2016; 2018), the following growth rate conditions.
Assumption 1 (High-dimensional asymptotics) As n → ∞, we have (i) p/n → c ∈ (0,∞) and na/n → ca ∈ (0, 1); (ii) ‖µa‖ = O(1); (iii) for C◦ = ∑K a=1 na n Ca and C ◦ a = Ca − C◦,
‖Ca‖ = O(1), tr(C◦a) = O( √ p) and tr(CaCb) = O(p) for a, b ∈ {1, · · · ,K}. We denote τ , tr(C◦)/p that is assumed to converge in (0,∞).
Remark 1 (Beyond Gaussian mixture data) While the theoretical results in this paper are derived for Gaussian mixture data in (7), under the non-trivial setting of Assumption 1, we conjecture that they can be extended to a much broader family of data distributions beyond the Gaussian setting, e.g., to the so-called concentrated random vector family (Seddik et al., 2020; Louart & Couillet, 2018), under similar non triviality assumptions. See Section A.1 in the appendix for more discussions.
To cover a large family of random features, we assume that the random projector matrix W has i.i.d. entries of zero mean, unit variance and bounded fourth-order moment, with no restriction on their particular distribution.
Assumption 2 (On random projection matrix) The random matrix W has i.i.d. entries such that E[Wij ] = 0, E[W 2ij ] = 1 and E[W 4 ij ] < λW for some constant λW <∞.
We consider the family of activation functions σ(·) satisfying the following assumption.
Assumption 3 (On activation function) The function σ is at least twice differentiable (in the sense of distributions when applied to a random variable having a non-degenerate distribution function, see Remark 3 in Section A.1 of the appendix), with max{E|σ(z)|,E|σ2(z)|,E|σ′(z)|,E|σ′′(z)|} < λσ for some constant λσ <∞ and z ∼ N (0, 1).

heading: 3 MAIN RESULT
text: Our objective is to characterize the high-dimensional spectral behavior of the centered and expected random features kernel matrix K defined in (2). It turns out, somewhat surprisingly, that under the non-trivial setting of Assumption 1, one may mentally picture the high-dimensional data vectors as (i) being asymptotically pairwise orthogonal (i.e., xTi xj → 0 for i 6= j as p→∞) and (ii) having asymptotically equal Euclidean norms (i.e., ‖xi‖2 → τ ), independently of the underlying class they belong to. As we shall see, this high-dimensional “concentration” of xTi xj → τ · δi=j plays a crucial role in “decoupling” the two dependent but asymptotically Gaussian random variables σ(wTxi) and σ(wTxj) in the definition (1). This, up to some careful control on the higher-order (but well “concentrated”) terms, leads to the following result on the asymptotic behavior of K, the proof of which is given in Section A.3 of the appendix.
Theorem 1 (Asymptotic equivalent of K) Under Assumption 1-3, for K defined in (2), as n→∞,
‖K− K̃‖ → 0,
2That is, when classification is neither too hard nor too easy.
almost surely with P = In − 1n1Tn/n,
K̃ = P ( d1 · ( Z + M
JT √ p
)T( Z + M
JT √ p
) + d2 ·VAVT + d0 · In ) P, (8)
and
V = [J/ √ p, φ] ∈ Rn×(K+1), A =
[ ttT + 2T t
tT 1
] ∈ R(K+1)×(K+1),
where, for z ∼ N (0, 1),
d0 = E[σ 2( √ τz)]− E[σ( √ τz)]2 − τE[σ′( √ τz)]2, d1 = E[σ ′( √ τz)]2, d2 = 1
4 E[σ′′(
√ τz)]2
and “first-order” random matrix Z = [z1, · · · , zn] ∈ Rp×n as defined in (7), “second-order” random (fluctuation) vector φ = { ‖zi‖2 − E[‖zi‖2] }n i=1 ∈ Rn, and GMM data statistics
M = [µ1, · · · ,µK ] ∈ Rp×K , t = {
tr(C◦a)√ p }K a=1 ∈ RK , T = { trCaCb p }K a,b=1 ∈ RK×K (9)
as well as the class label vectors J = [j1, · · · , jK ] ∈ Rn×K with [ja]i = δxi∈Ca .
First note that the “asymptotic equivalence” established in Theorem 1 holds for the expected kernel K, not on the empirical random feature kernel G defined in (6), and is thus independent of the number of random features m. On closer inspection of Theorem 1, we see that the data statistics for classification, i.e., the means (M) and covariances (ttT, T) are respectively weighted by the generalized Gaussian moments of first (d1) and second order (d2), while the coefficient d0 merely acts as a regularization term to shift all the eigenvalues,3 and has thus asymptotically no impact on the performance of, e.g., kernel spectral clustering for which only eigenvector structures are exploited.4
Theorem 1 unveils a surprising universal behavior of the expected kernel K in the large n, p regime. Specifically, the expression of K̃ (and thus the spectral behavior of K, see our discussion in the paragraph that follows) is universal with respect to the distribution of W and the activation function σ, when they are “normalized” to satisfy Assumptions 2 and 3. This technically improves previous efforts such as (Liao & Couillet, 2018b),5 and allows us to design computationally more efficient random features approach by wisely choosing W and σ.
As a direct consequence of Theorem 1, one has, by Weyl’s inequality and Davis–Kahan theorem that (i) the difference between each corresponding pair of eigenvalues and (ii) the distance between the “isolated” eigenvectors or the “angle” between the “isolated subspaces” of K and K̃ vanish asymptotically as n, p → ∞. This is numerically confirmed in Figure 2, where one observes a close match between the spectra (eigenvalue “bulk” and isolated eigenvalue-eigenvector pairs) of K and those of its asymptotic equivalent K̃ given in Theorem 1, already for n, p only in hundreds. In particular, we compare, in Figure 2, the eigenspectra of K and K̃ for Gaussian versus Student-t distributed W, on GMM data. A close match of the spectra and the isolated eigen-pairs is observed, irrespective of the distribution of W, as predicted by our theory.
In the following corollary, we exploit the universal result in Theorem 1 to design computationally efficient Ternary Random Features (TRFs) by specifying the distribution of (the entries of) W to symmetric Bernoulli and σ to a ternary function, so as to obtain a limiting kernel Kter that is asymptotically equivalent to any random features kernel matrix K of the form (2), for the highdimensional GMM data under study. This is proved in Section A.4 of the appendix.
3This can be seen as another manifestation of the implicit regularization in high-dimensional kernel and random features (Jacot et al., 2020; Derezinski et al., 2020; Liu et al., 2021b).
4We provide in Table 1 (Section A.5 of the appendix) the corresponding Gaussian moments d0, d1, d2 for various commonly used activation functions in random features and neural network contexts.
5 From a technical perspective, Theorem 1 improves (Liao & Couillet, 2018b, Theorem 1) in the fact that the latter relies on the explicit forms of the expectation K in (2), and is thus limited to (i) a few nonlinear σ for which K can be computed explicitly, see (Liao & Couillet, 2018b, Table 1); and (ii) Gaussian distributed W in which case the p-dimensional integral can be easily reduced to a two-dimensional one. Here, Theorem 1 holds for a much broader family of random W and nonlinear σ as long as Assumptions 2 and 3 hold.
Corollary 1 (Ternary Random Features) For a given random features kernel matrix K of the form (2) with W and nonlinear σ satisfying Assumptions 1-3, with associated generalized Gaussian moments d0, d1, d2 defined in Theorem 1, let σter be defined in (5) with s− = ŝ−, s+ = ŝ+, and ŝ−, ŝ+ satisfying the following equations
d1 = 1
π2
( e−ŝ 2 +/τ + e−ŝ 2 −/τ )2 , d2 = 1
2πτ3
( ŝ+e −ŝ2+/τ + ŝ−e −ŝ2−/τ )2 . (10)
Define the Ternary Random Features matrix Σter = σter(WterX) with Wter defined in (4) having sparsity level , the associated Gram matrix Gter = 1m (Σ ter)TΣter as in (6), and the limiting kernel
Kter , P{κter(xi,xj)}ni,j=1P (11)
for κter defined in (3). Then, there exists λ ∈ R such that6 ‖K−Kter − λP‖ → 0 almost surely as n→∞.
Note that the system of equations in (10) defining ŝ− and ŝ+ is a function of the key parameter τ = trC◦/p defined in Assumption 1, which can be consistently estimated from the data; see Algorithm 1 below and a proof in Lemma 1 of the appendix (the intuition of which follows from the fact that ‖xi‖2 = ‖µa‖2/p+ 2µTazi/ √ p+ ‖zi‖2 = E[tr(zizTi )] +O(p−1/2) according to (7) and Assumption 1). This makes the proposed TRFs and its limiting kernel Kter data-statistics-dependent. Yet, the system of equations (10) does not have a closed-form solution and might have multiple solutions on the real line; we practically solve it using a numerical least squares method, by gradually enlarging the search range (from say [−1, 1]) until a solution is found (the time complexity of which is independent of the dimension n, p and it is observed in the experiments in Section 4 to converge in a few iterations). The details of the proposed TRF approach are described in Algorithm 1.
Algorithm 1 Ternary Random Features Input: Data X and level of sparsity ∈ [0, 1). Output: Ternary Random Features Σter and Gram matrix Gter. Estimate τ as τ̂ = 1n ∑n i=1 ‖xi‖2.
Solve for thresholds ŝ−, ŝ+ using (10), which defines σter via (5). Construct a random matrix Wter ∈ Rm×p having i.i.d. entries distributed according to (4). Compute Σter = σter(WterX) and then TRFs Gram matrix Gter = 1m (Σ ter)TΣter as in (6).
6The parameter λ characterizes the possibly different d0 between K and Kter , see details in Appendix A.4.
Computational and storage complexity For W ∈ Rm×p a random matrix with i.i.d. N (0, 1) entries and Wter ∈ Rm×p with i.i.d. entries satisfying (4) with sparsity level ∈ [0, 1), let G = 1 mσ(WX)
Tσ(WX) for some given smooth function σ (e.g., sine and cosine in the case of random Fourier features) and Gter = 1mσ
ter(WterX)Tσter(WterX) with data matrix X ∈ Rp×n. It is then beneficial to use Gter instead of G as the computation of σ(WX) requires O(mnp) multiplications and O(mnp) additions, while the computation of σter(WterX) requires no multiplication and only O((1− )mnp) additions. In terms of storage, it requires a factor of b = 32 times more bits to store σ(WX) when computing G compared to storing σter(WterX) when computing Gter (assuming full precision numbers are stored using b = 32 bits).
The computationally and storage efficient TRFs σter(WterX) can then be used instead of the “expensive” random features σ(WX) and lead (asymptotically) to the same performance as the latter on downstream tasks, at least for GMM data, according to Theorem 1. In the following section, we provide empirical results showing that (i) this “performance match” between TRFs and random features of the form (1) is not limited to Gaussian data and empirically holds when applied on popular real-world datasets such as MNIST (LeCun et al., 1998), some UCI ML datasets, as well as DNNfeatures of CIFAR10 data in Section A.6 of the appendix; and (ii) due to the competitive performance of TRFs with respect to standard random features approaches, when compared to state-of-the-art random feature compression/quantization techniques (for which a “performance-complexity tradeoff” generally arises, that is, as one compresses more the original random features, the performance decays), TRFs yield significantly better performances for a given storage/computational budget.

heading: 4 EXPERIMENTS
text: The experiments in this section and Section A.6 of the appendix are performed on a Ubuntu 18.04 machine with Intel(R) Core(TM) i9-9900X CPU @ 3.50GHz and 64 GB of RAM.

heading: 4.1 TRFS MATCH THE PERFORMANCE OF POPULAR KERNELS WITH LESS BUDGET
text: We first consider ridge regression with random features Gram matrix G on MNIST data in Figure 3. We compare RFFs with σ(t) = [cos(t), sin(t)] and Gaussian Wij ∼ N (0, 1) to the proposed TRF method with σter(t) in (5) and ternary random projection matrix Wter defined in (4), with different sparsity levels . The thresholds s−, s+ of σter are tuned in such away that the generalized Gaussian moments d1 and d2 match with those of RFFs,7 as described in Corollary 1 and Algorithm 1. Figure 3 displays the test mean squared errors as a function of the regularization parameter γ for our TRF method with different sparsity levels compared to the RFF method, as well as to the baseline of Kernel Ridge Regression (KRR) using Gaussian kernel. Note that despite 90% sparsity in the projection matrix W, virtually no performance loss is incurred compared with RFFs. This is further confirmed in the right hand side of Figure 3 which shows the gains in running time8 when using TRFs. These experiments show that the proposed TRF approach yields similar performance as popular random features such as RFFs, with a significant reduction in terms of computation and storage.

heading: 4.2 COMPUTATIONAL AND STORAGE GAINS – COMPARISONS TO STATE-OF-THE-ART
text: In this section, we compare TRFs with state-of-the-art quantized or non-quantized random features methods, for both random features-based logistic and ridge regressions. Specifically, in Figure 1, we compare logistic regression performance of TRFs versus the Low Precision Random Fourier Features (LP-RFFs) proposed by Zhang et al. (2019), on the Cov-Type dataset from UCI ML repo. As in Section 4.1, the TRFs are tuned to match the limiting Gaussian kernel of RFFs. For a single datum x ∈ Rp, the associated TRFs σter(Wterx) ∈ Rm use m bits for storage while LP-RFFs use 8m bits. We follow the same protocol as in (Zhang et al., 2019) and use SGD with mini-batch size 250 in training logistic regressor and {1, 8} bits precision for the LP-RFF approach. Figure 1 compares the
7For given xi,xj , we use [cos(Wxi), sin(Wxj)] as the random Fourier features so that the (i, j) entry of the corresponding random features Gram matrix is cos(Wxi)T cos(Wxj)+sin(Wxi)T sin(Wxj) (Rahimi & Recht, 2008). The generalized Gaussian moments of the RFFs are thus the sum of the d1’s and d2’s corresponding to sin and cos functions. Note that the d0’s of RFFs and TRFs may be different.
8The running time is taken as the total clock-time of the whole ridge regression solver including the random features calculation.
Published as a conference paper at ICLR 2022
10−2 10−1 1 102 104 0
0.2
0.4
0.6
0.8
1
γ
M SE
s
Baseline (KRR) RFF
= 0.1 (TRF) = 0.5 (TRF) = 0.9 (TRF)
10−2 10−1 1 102 104
0
0.2
0.4
0.6
0.8
1
γ
M SE
s Baseline (KRR)
RFF = 0.1 (TRF) = 0.5 (TRF) = 0.9 (TRF)
0.1 0.3 0.5 0.7 0.9
1,000
2,000
R un
ni ng
tim e
(s )
Baseline (KRR) RFF TRF
Figure 3: Testing mean squared errors (MSEs with ±1 std, LEFT) and running time (RIGHT) of random features ridge regression as a function of regularization parameter γ, p = 512, n = 1024, ntest = 512,m = 5.10 4. With W distributed according to (4) and = [0.1, 0.5, 0.9], versus RFFs and KRR on a 2-class MNIST dataset – digits (7, 9). We find that d0 = 0.44 for TRFs and d0 = 0.39 for RFFs, making λ in Corollary 1 small. Results averaged over 5 independent runs.
logistic regression test accuracy as a function of the total memory budget for LP-RFF (1 bit and 8 bits) and the Nyström approximation of Gaussian kernel matrices (using 32 and 16 bits) (Williams & Seeger, 2001) (see also Table 1 in (Zhang et al., 2019)), versus the proposed TRFs approach. As seen from the shift in the x-axis (memory), by using 8× less memory than LP-RFFs and 32× or 16× less memory than the Nyström method, TRFs achieve a superior generalization performance on the logistic regression task. Similar comparisons are performed in Figure 4 on a random features ridge regression task for the Census dataset (Rahimi & Recht, 2008), confirming that TRFs outperform alternative approaches in terms of test mean square errors, with a significant save in memory.

heading: 5 CONCLUSION
text: Our large dimensional spectral analysis of the random features kernel matrix K reveals that its spectral properties only depend on the nonlinear activation through the corresponding generalized Gaussian moments and are universal with respect to zero-mean and unit-variance random projection vectors. This allows us to design the new TRF approach which turns both the random weights W and the activations σ(WX) into ternary integers, thereby allowing to only perform addition operations and to only store 1 bit for the activations. This drastically saves the storage and computation of random features while preserving the performances on downstream tasks with respect to their counterpart expensive kernels. Our article comes along with (Couillet et al., 2021) as first steps in re-designing machine learning algorithms using Random Matrix Theory, in order to be able to perform computations on massive data using desktop computers instead of relying on energy consuming giant servers.

heading: ACKNOWLEDGEMENT
text: Z L would like to acknowledge the National Natural Science Foundation of China (NSFC-12141107), the Fundamental Research Funds for the Central Universities of China (2021XXJS110), and CCFHikvision Open Fund (20210008) for providing partial support. R C would like to acknowledge the MIAI LargeDATA chair (ANR-19-P3IA-0003) at University Grenobles-Alpes, the HUAWEI LarDist project, as well as the ANR DARLING project for providing partial support of this work.

heading: A APPENDIX
text: We provide detailed discussions on our working assumptions in Section A.1, a consistent estimator of the key parameter τ in Section A.2, the proof of Theorem 1 in Section A.3 and that of Corollary 1 in Section A.4, as well as the generalized Gaussian moments for standard kernels arising from random features and neural net contexts in Section A.5. Additional numerical experiments are placed in Section A.6.

heading: A.1 NOTES ON THE WORKING ASSUMPTIONS
text: For completeness, let us restate the working settings of the paper.
Let x1, · · · ,xn ∈ Rp be n independent vectors belonging to one of K distributional classes C1, · · · , CK . Class a has cardinality na, and we assume that xi follows a Gaussian Mixture Model (GMM), i.e., for xi ∈ Ca,
xi = µa/ √ p+ zi (12)
with zi ∼ N (0p,Ca/p) for some mean µa ∈ Rp and covariance Ca ∈ Rp×p associated to class Ca. We position ourselves in the non-trivial regime of high-dimensional classification as described by the following growth rate conditions.
Assumption 4 (High-dimensional asymptotics) As n → ∞, we have (i) p/n → c ∈ (0,∞) and na/n → ca ∈ (0, 1); (ii) ‖µa‖ = O(1); (iii) for C◦ = ∑K a=1 na n Ca and C ◦ a = Ca − C◦, then
‖Ca‖ = O(1), tr(C◦a) = O( √ p) and tr(CaCb) = O(p) for a, b ∈ {1, · · · ,K}. We denote τ , tr(C◦)/p that is assumed to converge in (0,∞).
Beyond Gaussian mixtures While the theoretical results in this paper are derived for Gaussian mixture data under the non-trivial setting, we conjecture that they can be extended to a much broader family of data distributions beyond the Gaussian setting, e.g., to the so-called concentrated random vector family (Seddik et al., 2020), under similar non triviality assumptions.
Definition 1 (Concentrated vector) Given a normed vector space (X , | · |), and q > 0, a random vector x ∈ X is said to be q-exponentially concentrated if for any 1-Lipschitz function φ : X → R, there exists C > 0 independent of dim(X ) and σ > 0 such that for all t ≥ 0,
P (|φ(x)− E[φ(x)])| > t) ≤ Ce−(t/σ) q .
Multivariate Gaussian distributed x ∼ N (0, Ip) can be checked to belong to the family of concentrated random vectors. The major advantage of using concentrated random vectors for data modeling is that this concentration property is stable under Lipschitz transformation, that is, for any 1-Lipschitz map f : Rp → Rq, if x ∈ Rp is concentrated, so is x′ = f(x). Among the broad family of concentrated random vectors, it has been particularly shown in (Seddik et al., 2020) that artificial images generated by a Generative Adversarial Network (GAN) (Seddik et al., 2020), which look extremely close to real-world images (as they are designed to), are Lipshitz transformations of Gaussian vectors and can thus be, by definition, concentrated random vectors.
As a result, concentrated random vectors are more appropriate (than GMM for instance) in modeling realistic data, at least for those “close” to data generated by GANs.
The same was shown experimentally for CNN representations of real images (Seddik et al., 2019; 2020) as well as words embeddings in Natural Language Processing (Couillet et al., 2020, November).
As for the random projector matrix, we assume that the matrix W has i.i.d. entries of zero mean, unit variance and bounded fourth-order moment, with no restriction on their particular distribution as follows.
Assumption 5 (On random projection matrix) The random matrix W has i.i.d. entries such that E[Wij ] = 0, E[W 2ij ] = 1 and E[W 4 ij ] < λW for some constant λW <∞.
Remark 2 (On the i.i.d. assumption of entries of W) In Assumption 5 we consider the setting where the entries of W are independently and identically distributed: this is the case for the popular vanilla random Fourier features in (Rahimi & Recht, 2008) and the arc-cosine kernels in a neural network context; but not the case for, e.g., data dependent random features approaches such as leverage score based methods (Li et al., 2021).
We consider the family of activation functions σ(·) satisfying the following assumption.
Assumption 6 (On activation function) The function σ is at least twice differentiable (in the sense of distributions when applied to a random variable having a non-degenerate distribution function), with max{E|σ(z)|,E|σ2(z)|,E|σ′(z)|,E|σ′′(z)|} < λσ for some constant λσ < ∞ and z ∼ N (0, 1).
Remark 3 (Activation functions not differentiable everywhere) Some popular activation functions used in machine learning such as ReLu, Sign, Absolute value, etc., are not differentiable everywhere. For those functions, we will use a derivative in the sense of distributions (Friedlander et al., 1998). A distribution g is a continuous linear functional on the set D of infinitely differentiable functions with bounded support
g : D → R φ 7→ g(φ) = ∫ +∞ −∞ g(x)φ(x) dx.
The distributional derivative g′(φ) is defined such that g′(φ) = −g(φ′). In particular, we will be interested here in the expectation of the derivatives of the activation function with respect to the Gaussian measure i.e., ∫ +∞ −∞ σ ′(x)e−x 2/2 dx. Following the previous definition, we have in this particular case ∫ +∞ −∞ σ′(x)e−x 2/2 dx = ∫ +∞ −∞ xσ(x)e−x 2/2 dx
which can be evaluated by some integration by parts. This is also refereed to the as “weak derivative” in the functional analysis literature (Stein & Shakarchi, 2012).

heading: A.2 AUXILIARY RESULTS AND PROOFS
text: Lemma 1 (Consistent estimation of τ ) Let Assumption 4 hold and define τ , trC◦/p. Then as n→∞, with probability 1,
1
n n∑ i=1 ‖xi‖2 − τ → 0.
Proof 1 (Proof of Lemma 1) From equation 7, we have that
1
n n∑ i=1 ‖xi‖2 = 1 n K∑ a=1 n∑ i=1 1 p ‖µa‖2 − 2 √ p µTazi + 1 n n∑ i=1 ‖zi‖2. (13)
From Assumption 4, we have that 1n ∑K a=1 ∑n i=1 1 p‖µa‖
2 = O(p−1). The second term of equation 13 2√ pµ T azi is a weighted sum of independent zero mean random variables; it thus vanishes with probability 1 as n, p→∞ by a mere application of Chebyshev’s inequality and the Borell Cantelli lemma. Finally, using the strong law of large numbers on the last term of equation 13, we have almost surely,
1
n n∑ i=1 ‖zi‖2 = K∑ a=1 na n trCa + o(1)
= K∑ a=1 na n 1 p trC◦ + o(1)
where in the last line we use trC◦a = O( √ p) from Assumption 4, and thus 1n ∑n i=1 ‖zi‖2 − τ → 0 almost surely. This concludes the proof.

heading: A.3 PROOF OF THEOREM 1
text: In the sequel, we will make use of Bachmann–Landau notations but specific for random variables in the asymptotic regime where n → ∞. For x ≡ xn a random variable and un ≥ 0, we write xn = O(un) if for any η and D > 0, nDP(x ≥ nηun)→ 0 as n→∞. For a vector v or a diagonal matrix with random entries, v = O(un) means that the maximal entry of v in absolute value is O(un) in the sense defined above. When M is a square matrix, M = O(un) means that the operator norm of M is O(un). For x a random variable, x = o(un) means that there exists κ > 0 such that x = O(n−κun). And the same definition for the o(·) notation applies for vectors and matrices having random entries.
Let us define X = [x1, · · · ,xn] ∈ Rp×n and define
Σ = σ(WX)
where W = [w1, · · · ,wm]T ∈ Rm×p with W satisfying Assumption 5, and σ a function satisfying Assumption 6. We consider the gram matrix
G = 1
m ΣTΣ
whose (i, j) entry is given by
Gij = 1
m m∑ k=1 σ(wTkxi)σ(w T kxj).
We are interested in computing
κ(xi,xj) = Ew[σ(w Txi)σ(w Txj)] (14)
under Assumptions 4-6.
Under Assumption 4, we have xTi xj = z
T i zj︸︷︷︸
O(p− 1 2 )
+µTaµb/p+ µ T awj/ √ p+ µTbwi/ √ p︸ ︷︷ ︸
O(p−1)
(15)
and ‖xi‖2 = τ︸︷︷︸
O(1) + tr(C◦a)/p+ Φi︸ ︷︷ ︸ O(p− 1 2 ) + ‖µ2a‖/p+ 2µTazi/ √ p︸ ︷︷ ︸ O(p−1)
(16)
where φi , (‖zi‖2 − E[‖zi‖2]). It can be checked that, for random vector w ∈ Rp having i.i.d. entries with E[wi] = 0 and E[w2i ] = 1, we have, conditioned on xi, that Ew[(wTxi)2] = ‖xi‖2 and
Ew[(w Txi) 4] = (m4 − 3)‖xi‖2 + 2‖xi‖4. (17) with m4 = E[w4j ] < λW according to Assumption 5. From the trace lemma, (Bai & Silverstein, 2010, Lemma B.26), one has that
xTi xi − 1
p trCa → 0 (18)
almost surely as p → ∞, for xi ∈ Ca. Thus, under Assumption 4 we have in particular xTi xi − trC◦/p → 0 holds almost surely, regardless of the class of xi, see the proof of Lemma 1. It then follows from Lyapunov CLT (see, e.g., (Billingsley, 2012, Theorem 27.3)) that, under Assumption 4 and 5, (wTxi,wTxj) is asymptotically bivariate Gaussian. We can thus perform, in the large p limit, a Gram-Schmidt orthogonalization procedure for some standard Gaussian variables ξa, ξb ∼ N (0, 1). By construction, ξa, ξb are uncorrelated and thus independent by Gaussianity. Let us denote the shortcuts ζa = wTxi and ζb = wTxj . We thus have
ζa ≡ wTxi = uaξa + o(1) ζb ≡ wTxj = vbξa + ubξb + o(1)
with ua = ‖xi‖
vb = xTi xj ‖xi‖
ub = √ ‖xj‖2 − (xTi xj) 2
‖xi‖2 .
Since we have that ‖xi‖ = √ τ + o(1) and ‖xj‖ = √ τ + o(1) (from equation 16), we can perform
a Taylor expansion of σ(ζa) (respectively of σ(ζb)) around √ τξa (resp. √ τξb) giving
σ(ζa) = σ( √ τξa) + σ ′( √ τξa)(ζa − √ τξa) + 1
2 σ′′(ξa)(ζa −
√ τξa) 2 + o((ζa − √ τξa) 2)
σ(ζb) = σ( √ τξb) + σ ′( √ τξb)(ζb − √ τξb) + 1 2 σ′′( √ τξb)(ζb − √ τξb) 2 + o((ζb − √ τξb) 2).
We then have κ(xi,xj) = E[σ(ζa)σ(ζb)]
= E[σ( √ τξa)σ( √ τξb)] + E[σ( √ τξa)σ ′( √ τξb)(ζb − √ τξb)]
+ 1 2 E[σ( √ τξa)σ ′′( √ τξb)(ζb − √ τξb) 2] + E[σ′( √ τξa)(ζa − √ τξa)σ( √ τξb)] + E[σ ′( √ τξa)(ζa − √ τξa)σ ′( √ τξb)(ζb − √ τξb)]
+ 1
2 E[σ′(ξa)(ζa −
√ τξa)σ ′′(ξb)(ζb − √ τξb) 2] + 1
2 E[σ′′(ξa)(ζa −
√ τξa) 2σ( √ τξb)]
+ 1
2 E[σ′′(
√ τξa)(ζa − √ τξa) 2σ′( √ τξb)(ζb − √ τξb)]
+ 1
4 E[σ′′(
√ τξa)(ζa − √ τξa) 2σ′′( √ τξb)(ζb − √ τξb) 2] + o((ζa − √ τξa) 2).
where the expectation is with respect to w as in the definition in (14).
Using the independence between ξa, ξb along with the following
ζa − √ τξa = ( ua√ τ − 1) √ τξa (19) ζb − √ τξb = (
ub√ τ − 1) √ τξb + vb√ τ √ τξa (20)
we have for ξ ∼ N (0, 1),
E[σ( √ τξa)σ( √ τξb)] = ( E[σ( √ τξ)] )2 (21)
E[σ( √ τξa)σ ′( √ τξb)(ζb − √ τξb)] = ( E[σ( √ τξ)]E[ √ τξσ′( √ τξ)] )( ub√ τ − 1 )
+ ( E[σ′( √ τξ)]E[ √ τξσ( √ τξ)] ) vb√ τ
(22)
Similarly writing down the products (ζa − √ τξa) 2, (ζb − √ τξb) 2, (ζa − √ τξa)(ζb − √ τξb), (ζa −√
τξa)(ζb − √ τξb) 2, (ζa − √ τξa) 2(ζb − √ τξb), (ζa − √ τξa) 2(ζb − √ τξb)
2 as in Equations 21- 22, and using Equations 19- 20, we obtain a complete expression of all the terms in κ(xi,xj) as follows
κ(xi,xj) =
[( E[σ( √ τξ)] )2 + ( E[σ( √ τξ)]E[ √ τξσ′( √ τξ)] )( ub√ τ − 1 ) + ( E[σ′( √ τξ)]E[ √ τξσ( √ τξ)] ) vb√ τ
+ ( E[ √ τξσ′( √ τξ)]2 )( ( ua√ τ − 1)( ub√ τ − 1) ) + ( E[σ( √ τξ)]E[τξ2σ′′( √ τξ)] ) 2 ( (ub − 1)2 ) + ( E[ √ τξσ′( √ τξ)]E[τξ2σ′′( √ τξ)] ) 2 ( ua√ τ − 1 )( ub√ τ − 1 )2
+
( E[ √ τξσ′( √ τξ)]E[τξ2σ′′( √ τξ)] ) 2 ( ua√ τ − 1 )2( ub√ τ − 1 )
+
( E[τξ2σ′′( √ τξ)]2 ) 4 ( ua√ τ − 1 )2( ub√ τ − 1 )2 + ( E[ √ τξσ( √ τξ)]E[ √ τξσ′′( √ τξ)] ) vb√ τ ( ub√ τ − 1 )
+
( E[σ′′( √ τξ)]E[τξ2σ′′( √ τξ)] ) 2 ( vb√ τ )2 + ( E[σ( √ τξ)]E[ √ τξσ′( √ τξ)] )( ua√ τ − 1 )
+ ( E[σ′( √ τξ)]E[ξ2σ′( √ τξ)] ) vb√ τ ( ua√ τ − 1 ) + ( E[σ′′( √ τξ)]E[τ √ τξ3σ′( √ τξ)] ) 2 ( vb√ τ )2( ua√ τ − 1 )
+
( E[σ′( √ τξ)]E[τξ2σ′( √ τξ)] ) 2 vb√ τ ( ua√ τ − 1 )( ub√ τ − 1 ) + ( E[σ( √ τξ)]E[τξ2σ′′( √ τξ)] ) 2 ( ua√ τ − 1 )2
+
( E[σ′( √ τξ)]E[τ √ τξ3σ′′( √ τξ)] ) 2 vb√ τ ( ua√ τ − 1 )2
+
( E[σ′′( √ τξ)]E[τ2ξ4σ′′( √ τξ)] ) 4 ( vb√ τ )2( ua√ τ − 1 )2
+
( E[ √ τξσ′′( √ τξ)]E[τ √ τξ3σ′′( √ τξ)] ) 2 vb√ τ ( ub√ τ − 1 )( ua√ τ − 1 )2] +O(p− 3 2 ).
(23)
Since |xTi xj | ≤ for a sufficiently small , (following from equation 15), using a Taylor approximation of √ ‖xj‖2 − (xTi xj) 2
‖xi‖2 we obtain
ub√ τ − 1 = 1√ τ
√ ‖xj‖2 − (xTi xj) 2
‖xi‖2 − 1 = ( ‖xj‖√ τ − 1 ) − 1 2 √ τ‖xj‖ (xTi xj) 2 ‖xi‖2 + o ( (xTi xj) 2 ‖xi‖2 ) ( ub√ τ − 1 )2 = ( 1√ τ √ ‖xj‖2 − (xTi xj) 2 ‖xi‖2‖xj‖2 − 1 )2 = ( ‖xj‖√ τ − 1 )2 − ‖xj‖√ τ − 1 ‖xj‖ (xTi xj) 2 ‖xi‖2 + o ( (xTi xj) 2 ‖xi‖2
) ua√ τ − 1 = ‖xi‖√ τ − 1.
We thus have
vb√ τ = xTi xj√ τ‖xi‖ = 1 τ  zTi zj︸︷︷︸ O(p− 1 2 ) + ( µTaµb p + µTawj√ p + µTbwi√ p ) ︸ ︷︷ ︸
O(p−1) +O(p− 32 ) ( ua√ τ − 1 ) = ‖xi‖√ τ − 1 = 1 2τ ( tr(C◦a) p + Φi ) ︸ ︷︷ ︸
O(p− 1 2 )
+ 1
2τ
( ‖µa‖2
p + 2µTa zi√ p ) ︸ ︷︷ ︸
O(p−1)
+O(p− 3 2 )
( ua√ τ − 1 )( ub√ τ − 1 ) = ( ‖xi‖√ τ − 1 )( ‖xj‖√ τ − 1 ) =
( tr(C◦a) p + Φi )( tr(C◦b ) p + Φj ) 4τ2︸ ︷︷ ︸
O(p− 1 2 )
+O(p− 3 2 )
( vb√ τ )2 = 1 τ2 (zTi zj) 2 +O(p− 3 2 ). (24)
All other terms in equation 23 are of order at most O(p− 3 2 ) and thus vanish asymptotically. We thus get
κ(xi,xj) =
[( E[σ( √ τξ)] )2 + ( E[σ( √ τξ)]E[ √ τξσ′( √ τξ)] )( ub√ τ − 1 ) + ( E[σ′( √ τξ)]E[ √ τξσ( √ τξ)] ) vb√ τ +
+ ( E[ √ τξσ′( √ τξ)]2 )( ( ua√ τ − 1)( ub√ τ − 1) ) + ( E[σ′′( √ τξ)]E[τξ2σ( √ τξ)] ) 2 ( vb√ τ )2 + ( E[σ( √ τξ)]E[ √ τξσ′( √ τξ)] )( ua√ τ − 1 ) +O(p− 3 2 ) ] . (25)
Plugging the terms in equation 24 into equation 25 and rearranging in matrix form, we obtain for K = P{κ(xi,xj)}ni,j=1P and K̃ defined in Theorem 1, that ‖K− K̃‖ → 0, as expected, where we used the fact that ‖A‖ ≤ n‖A‖∞ for A ∈ Rn×n and ‖A‖∞ = maxni,j=1 |Aij |. This concludes the proof of Theorem 1.

heading: A.4 PROOF OF COROLLARY 1
text: Define the Ternary Random Features matrix Σter = σter(WterX) with Wter defined in (4) having sparsity level , the associated Gram matrix Gter = 1m (Σ ter)TΣter as in (6), and the limiting kernel
Kter , P{κter(xi,xj)}ni,j=1P (26)
for κter defined in (3). After calculation similar to (Liao & Couillet, 2019, Section 4) and (Liao et al., 2021), we obtain
E[ ( σter )′ ( √ τz)]2 = e− s2+τ + e− s2−τ π 2
E[ ( σter )′′ ( √ τz)]2 = s+e− s2+τ + s−e− s2−τ τ √ 2πτ 2 . (27) Thus, according to Theorem 1, for random features kernel matrix K of the form (2) with W and nonlinear σ satisfying Assumption 1-3, with associated generalized Gaussian moments d0, d1, d2 defined in Theorem 1, by choosing s−, and s+ such that d1 and d2 are respectively equal to the first and second term of equation 27, we have that
‖K−Kter − λP‖ → 0,
almost surely with λ = d0 − ( E[ ( σter )2 ( √ τz)]− E[ ( σter ) ( √ τz)]2 − τE[ ( σter )′ ( √ τz)]2 ) where
E[ ( σter )2 ( √ τz)]− E[ ( σter ) ( √ τz)]2 − τE[ ( σter )′ ( √ τz)]2
= ( 1− erf ( s+√ τ )) − 1 2 e− s2+τ + e− s2−τ π 2 .

heading: A.5 GAUSSIAN MOMENTS OF POPULAR ACTIVATION FUNCTIONS
text: We provide in Table 1 the calculation of the generalized Gaussian moments d0, d1, d2 for popular activation functions used in random features and neural network contexts. Note that our Table 1 matches the results in (Liao & Couillet, 2018b, Table 2).

heading: A.6 ADDITIONAL EXPERIMENTS
text: In this section, we complement Section 4 by providing additional experiments on ridge regression and support vector machine, to support the robustness our TRF method compared to state-of-the-art approaches.
Remark 4 (On the empirical and expected random features kernels) It is worth noting that our main results in Theorem 1 and Corollary 1 characterize the behavior of the expected/ limiting random features kernel K instead of the empirical Gram kernel matrix G defined in (6) obtained by averaging over m random features. The operator norm difference ‖K−G‖ is then bound to vanish as m, p, n→∞ for a sufficiently large m (for example with m/max(n, p)→∞).
Random features based Ridge regression To complement the experiments in Section 4.1, we provide in Figures 5-9, the test mean square error of random features kernel ridge regression with increasing number of random features m ∈ {512, 4096, 104} for GMM data in Figure 5-7, and m ∈ {512, 104} for MNIST data in Figure 8 and 9, as a function of the regularization parameter γ and for different choices of sparsity levels ∈ {0.1, 0.3, 0.5, 0.7, 0.9}. It is interesting to note that, for both datasets, the performance gap between RFFs and TRFs significantly decreases as m the number of random features grows large: this is in agreement with our Theorem 1 in which guarantee is provided only for the expected kernel matrix (that corresponds to m → ∞), not the empirical Gram matrix as the sample mean over m random features.
Random features based Support Vector Machine We empirically evaluate the classification performance of various random features approximation algorithms, on several benchmark datasets. We compared the different algorithms (RFF (Rahimi & Recht, 2008), ORF (Yu et al., 2016), SSF (Lyu, 2017)) on 3 datasets IJCNN1 (Prokhorov, 2001), Cov-Type, and a8a from the UCI ML repository considered in (Liu et al., 2021a), with our TRF method where we choose the thresholds coefficients s−, s+ according to Algorithm 1 (to match the generalized Gaussian moments of Gaussian kernel). Figure 10 shows the results for the a8a dataset, Figure 11 for the IJCNN1 dataset and Figure 12 for the Cov-Type dataset. The lower running time along with higher SVM test accuracy indicates the superiority of our RMT-inspired TRF method over the other Gaussian kernel approximation methods.
Comparison of TRF with equivalent 0th order Arc-cosine kernel (with ReLU function) We consider Support Vector Machine (SVM) classification with random features Gram matrix G on Fashion MNIST data (LeCun et al., 1998) and on VGG16 embeddings of CIFAR10 (Krizhevsky et al., 2009) and Imagenet (Deng et al., 2009) datasets in Figures 14 -13 -15 respectively. We use VGG16 with batch normalization (Ioffe & Szegedy, 2015) pre-trained on ImageNet (Deng et al., 2009) as a feature extractor. We fine-tune this model on the CIFAR10 dataset with 240 epochs and a mini-batch size 64 with a SGD optimizer with momentum 0.9 and an initial learning rate of 0.1. We then extract the output of the first fully connected layer of the classifier as our features. For Imagenet, we directly extract the features from the pretrained model. We compare (i) RF with σ(t) = max(t, 0) (ReLU) and Gaussian Wij ∼ N (0, 1) (known as the 0th order Arc-cosine kernel) to (ii) the proposed TRF
method with σter(t) in (3) and ternary random projection matrix Wter defined in (4). The thresholds s−, s+ of σter are tuned in such away that the generalized Gaussian moments d1 and d2 are matched with those of ReLU (see Table 1), as described in Corollary 1 and Algorithm 1. Figures 14 -13 -15 display the test SVM accuracy as a function of the number of random features, for our TRF method compared to the random features corresponding to the Arc-Cosine kernel. We observe similar test performances with the two kernels while having a computation and storage gains for the ternary kernel.

Review summaries:
The authors present a resource-efficient approach for approximating kernels using random feature transformations. The random feature approach was introduced in a seminal paper by Rahimi and Recht, where each random feature is extracted from the input $x$ as $\sigma(\langle w, x\rangle )$, for a random (appropriately distributed) vector $w$ and an activation function (e.g. sinusiod) $\sigma$. In this paper, the goal is to save memory and time by restricting $w$ and $\sigma$ to take values in {$-1,0,1$}, meaning that $\sigma(\langle w, x\rangle )$:  1. can be evaluated only with additions/subtractions (no multiplications) and 2. only takes $1$ bit of storage.   The authors prove that under some conditions, their kernel estimate approaches the true kernel at the limit where the number of features and the number of data points are comparable and approach infinity. The assumptions are roughly: 1. The data follows a Gaussian mixture model, 2. the input features are approximately pairwise orthogonal, 3. the entries of $w$ are chosen i.i.d. and have bounded fourth moment, and 4. some boundedness assumptions on the generalized Gaussian moments of $\sigma$.  The theoretical result is accompanied by numerical experiments on various tasks like kernel ridge regression and SVM on both real and synthetic datasets. The results demonstrate a significant advantage of the authors' approach compared to previous quantized random features approches, both in terms of runtime and memory, without a significant drop in accuracy.

This paper studies random features method, which randomly projects the data onto a low-dimensional space, while at the mean time approximates the original kernel structure.  Some contributions of this paper: - The authors show that, under the Gaussian mixture model, in the high-dimensional limit, in the spectrum sense, the random features-type kernel (defined in equation (1) of the paper) only depends on the first two generalized Gaussian moments of the activation function (with the assumption that the random projection weights are i.i.d. from some mean 0 and var 1 distribution). - Based on the above result, the authors propose a special choice of the random weights and the activation function, that could be efficiently computed and needs less memory to store. - The authors provide empirical evidence showing that the proposed method performs as well as other random features methods, and with advantage in computation and storage.

The paper proposes a way to approximate common kernel matrices in a way that is (1) Sparse (2) Low Bit-Complexity (i.e. the nonzero entries have few bits) (3) efficient to compute.  The theoretical justification from this model comes from a random matrix theory (RMT) analysis of kernel matrices for datasets of vectors drawn iid from a fairly general gaussian mixture model. The RMT analysis shows that the true kernel matrix is asymptotically distributed as the approximation described above. No non-asymptotic results are given.  The introduction mentions deep networks, but this paper really has no strong connection to deep learning. View it as a kernel paper.  A good amount of supporting experiments are given to demonstrate that the approximate kernel matrix has the same statistical performance as (eg) a Random Fourier Features approximation.

With Gaussian mixture assumption for input data, this paper proposed a sparse random features approach where the approximated kernel is independent of the iid weights and depends on Gaussian moments of the activation functions. This paper provided a theoretical guarantee for the asymptotical equivalence with the centered kernel. And then, the authors derived the computationally efficient random features and devised a simple algorithm. Finally, they validate the accuracy and efficiency of the proposed random features by several experiments.