Random Feature Mapping with Signed Circulant Matrix Projection

1 Introduction

Support Vector Machine (SVM) is one of the most popular classification tools, which is based on statistical learning theory and delivers excellent results for non-linear classification in machine learning [Vapnik, 1998; Sch¨olkopf and Smola, 2002]. The underlying training problem can be formulated as a Quadratic Programming (QP) problem that can be solved by standard optimization algorithms in O(l3 ) time and O(l2 ) space complexity, where l is the number of training data [Tsang et al., 2005]. This is computationally infeasible for training on very large datasets.

There is a lot of work proposed for scaling up non-linear kernel SVMs on large scale datasets, such as solving the QP problem exactly through decomposition methods [Platt, 1999; Chang and Lin, 2011], approximation algorithms using core vector set [Tsang et al., 2005], and parallel interior point solver [Chang et al., 2007]. However, the increasing volumes of datasets would lead to the curse of support [Kar and Karnick, 2012], where the number of Support Vectors (SVs) that have to be explicitly maintained grows linearly with the sample size on noisy data [Steinwart, 2003]. According to Representer Theorem [Kimeldorf and Wahba, 1970], and KarushKuhn-Tucker conditions [Boyd and Vandenberghe, 2004], one can typically represent the decision function f(x) via the kernel trick with SVs, f(x) = ∑xi∈SVs αik(xi, x), αi > 0. Obviously, it requires huge computational time and space to train a model as SVs increase on large scale datasets. Moreover, to evaluate a learned model with hundreds of thousands of SVs will bring in additional time and space burden in the predicting stage.

∗Corresponding author

Random feature mapping methods, such as Random Kitchen Sinks (RKS) [Rahimi and Recht, 2007; 2008] and Random Maclaurin Feature Maps [Kar and Karnick, 2012], are proposed for addressing the curse of support, which embed the implicit kernel feature space into a relatively lowdimensional explicit Euclidean space. In the embedded feature space, the kernel value of any two points is well approximated by their inner product and one can apply existing fast linear algorithms, such as linear SVMs that run in time linear with sample size [Joachims, 2006; Fan et al., 2008], to abstract data relations corresponding to non-linear kernel methods. After learning a hyperplane with a linear classifier, one can predict an input in O(dD) time complexity (independent on the number of training data and mainly used for computing random feature mapping), where d represents the dimensionality of the input and D the dimensionality of the random feature space. Therefore, random feature mapping methods with liner classifiers can inherit high efficiency of linear learning algorithms and good generalization performance of nonlinear kernel methods, providing a promising way to the curse of support. However, the random feature mapping itself is a bottleneck when dD is not small.

Fastfood uses an approximation of the unstructured Gaussian matrix of RKS with several special matrices to accelerate the computation of random feature mapping [Le et al., 2013]. Because the special matrices are easy to store and multiply, Fastfood computes the random feature mapping in O(D log d) time and O(D) space complexity, a significant improvement from O(dD) computation and storage. However, Fastfood brings great increasing variance of approximating the kernel, which will cause inaccurate approximation and loose concentration bound.

In this paper, we propose a novel random feature mapping method, called Signed Circulant Random Feature mapping (SCRF), for approximating non-linear kernel. Our proposed SCRF uses the structured random matrix Π ∈ RD×d , a stacking of D/d signed circulant Gaussian matrices, to project input data, followed by a non-linear transform, i.e., φ(Πx). We prove that the approximation using SCRF is an unbiased estimate of the corresponding non-linear kernel and does not increase the variance of that using RKS. Therefore, the approximation using SCRF concentrates with the same rate as RKS and faster than that using Fastfood. Because Π has the circulant structure, instead of saving it directly, we can store the first columns of the corresponding circulant Gaussian matrices with a random sequence in O(D) space complexity and then easily recover it. For the structured matrix projection, we can apply the Fast Fourier Transform (FFT) to compute the feature mapping in O(D log d) time complexity. Experimental results demonstrate that our proposed SCRF is much more time and space efficient than both RKS and Fastfood while retaining similar accuracies with state-of-the-art random feature mapping methods.

1.1 RelatedWork

Random feature mappings have been successfully used for approximating non-linear kernels to address the curse of support. As the first proposed random feature mapping method, Random Kitchen Sinks (RKS) [Rahimi and Recht, 2007] focuses on approximating translation invariant kernels (e.g., Gaussian kernel, Laplacian kernel). There have been several approaches proposed to approximate other kernels as well, including intersection [Maji and Berg, 2009], group invariant [Li et al., 2010], additive [Vedaldi and Zisserman, 2012] and dot product kernels [Kar and Karnick, 2012]. Moreover, there has been work aiming at improving quality of the existing random feature mapping methods, some of which include a more compact representation of accurately approximating polynomial kernels [Hamid et al., 2014] and a more effective Quasi-Monte Carlo feature mapping for translation invariant kernels [Yang et al., 2014]. Yen et al. [Yen et al., 2014] propose a sparse random feature algorithm so that the resulting model doesn't grow linearly with the number of random features.

Recently, two promising quasilinear kernel approximation techniques have been proposed to accelerate the existing random feature mappings. Tensor sketching applies recent results in tensor sketch convolution to deliver approximations for polynomial kernels in O(d + D logD) time [Pham and Pagh, 2013]. Fastfood uses an approximation of the unstructured Gaussian matrix of RKS with several special matrices to project input data [Le et al., 2013]. Because the special matrices are inexpensive to multiply and store, Fastfood computes the random feature mapping in O(D log d) time and O(D) space complexity.

We summarize the contributions of this paper as follows:

• We propose a novel scheme for random feature mapping that uses structured randommatrix (signed circulant ran-
• dom matrix) instead of unstructured random matrix to project data.
• We prove that the approximating Gaussian kernel using our proposed SCRF is unbiased and does not increase the variance of that using RKS.
• We save the random parameters in O(D) space complexity and implement our random feature mapping in O(D log d) time complexity by using FFT.
• Empirically, our random feature mapping method is time and space efficient, making kernel methods scalable for large scale training and predicting problems.

2 Preliminaries

In this section, we first review two well-known random feature mapping methods, Random Kitchen Sinks (RKS) [Rahimi and Recht, 2007] and Fastfood [Le et al., 2013], for kernel approximation, and then introduce a structured matrix, circulant matrix [Davis, 1979; Gray, 2006].

2.1 RandomFeature Mapping

The starting point of RKS is a celebrated result that characterizes the class of positive definite functions.

Theorem 1 (Bochner's theorem [Rudin, 2011]). For any normalized continuous positive definite function f : Rd → C, there exists a finite non-negative Borel measure µ on Rd such that

i.e. f is the Fourier transform of a finite non-negative Borel measure µ on Rd .

Without loss of generality, we assume henceforth that the kernel κ(δ) is properly scaled and µ(·) is a probability measure with associated probability density function p(·).

Corollary 1. For shift-invariant kernel k(x, y) = κ(x-y),

where p(w) is a probability density function and can be calculated through the inverse Fourier transform of κ.

For Gaussian kernel, k(x, y) = exp(-‖x - y‖2/(2σ2 )), we calculate p(w) through the inverse Fourier transform of k(x, y) and obtain w ∼ N(0, I/σ2 ), where I is an identity matrix. In addition, we have Ew[sin(wT(x - y))] = 0 and Ew,b[cos(wT(x + y) + 2b)] = 0, where b is drawn from [-π, π] uniformly. Note that

so 〈Zw,b(x),Zw,b(y)〉 is an unbiased estimate of the Gaussian kernel. Through a standard Monte Carlo (MC) approximation to the integral representation of the kernel, we can lower the variance of 〈Zw,b(x),Zw,b(y)〉 by concatenating D randomly chosen Zw,b into a column feature mapping vector and normalizing each component by √D,

where Z ∈ RD×d is a Gaussian matrix with each entry drawn i.i.d. from N(0, 1/σ2 ), b ∈ RD is a random vector drawn i.i.d. from [-π, π] uniformly and cos(·) is an element-wise function.

As derived above, the associated feature mapping converges in expectation to the Gaussian kernel. In fact, convergence occurs with high probability and at the rate of independent empirical averages. In the explicit random feature space, we can use primal space methods for training, which delivers a potential solution to the curse of support. However, such approach is still limited by the fact that we need to store the unstructured Gaussian matrix Z and, more importantly, we need to compute Zx for each x. Because the unstructured matrix has the disadvantage that no fast matrix multiplication is available, large scale problems are not practicable with such unstructured matrix.

Fastfood finds that Hadamard matrix H when combined with binary scaling matrix B, permutation matrix Q, diagonal Gaussian matrix G and scaling matrix S exhibits properties similar to the unstructured Gaussian matrix Z, i.e., V ≈ Z, where

The decomposed matrices are inexpensive to multiply and store, which speeds up the random feature mapping. However, approximating Gaussian kernel using Fastfood brings increasing variance, which causes inaccurate approximation and loose concentration bound. Different from Fastfood, in this paper, we make use of a structured random matrix to take the place of Z and propose a novel scheme for random feature mapping.

2.2 Circulant Matrix

A circulant matrix C is an m × m Toeplitz matrix with the form

where each column is a cyclic shift of its left one [Davis, 1979]. The structure can also be characterized as follows: For the (k, j) entry of C, Ckj ,

A circulant matrix is fully determined by its first column, so we rewrite the circulant matrix of order m as follows

Therefore, it only needs to store the first column vector so that we can reconstruct the whole matrix, which saves a lot of storage.

Definition 1 (Circulant RandomMatrix, CRM). A circulant matrix C[m] = circ [cj : j ∈ {0, 1, . . . ,m- 1}] is called a circulant random matrix if its first column is a random sequence with each entry drawn i.i.d according to some distribution probability.

Definition 2 (Signed Circulant Random Matrix, Signed CRM). If a matrix P = [σ0C0·; σ1C1·; . . . ; σm-1C(m-1)·] satisfies that Ci· is the i-th row vector of a circulant random matrix C[m] and σi is a Rademacher variable (P[σi = 1] = P[σi = -1] = 1/2), where i = 0, 1, . . . , m - 1, we call P a signed circulant random matrix.

The following lemma provides an equivalent form of a circulant matrix with the Discrete Fourier Transform (DFT) [Davis, 1979; Tyrtyshnikov, 1996].

Lemma 1. Suppose C is a matrix with the first column c = [c0, c1, . . . , cm-1]T. Then C is a circulant matrix, i.e.

if and only if

where

is the discrete Fourier matrix of order m (i = √-1), diag(Fc) signifies the diagonal matrix whose diagonal comprises the entries of the vector Fc, and F ∗ represents the conjugate transpose of F .

Obviously, we could realize the calculation of Cx efficiently via the Fast Fourier Transform algorithm (FFT). In the following, we introduce our structured feature mapping method and apply such computational advantage to accelerate the feature mapping.

3 Signed Circulant RandomFeature Mapping

In this section, we propose a random feature mapping method with signed circulant matrix projection, called Signed Circulant Random Feature mapping (SCRF), and analyse its approximation error.

3.1 Our RandomFeature Mapping

Without loss of generality, we assume that d | D (D is divisible by d). And then we generate t = D/d signed circulant Gaussian matrices P (1) ,P (2) , . . . ,P (t) and stack them by row to take the place of the unstructured Gaussian matrix Z in Eq. (4), i.e., replacing Z with projection matrix Π = [P (1) ;P (2) ; · · · ;P (t) ]. The first column of the corresponding circulant random matrix C(i) [d] of P (i) (i = 1, 2, . . . , t) is randomly drawn i.i.d. from a Gaussian distribution N(0, I/σ2 ). Therefore, we obtain a novel random feature mapping, called Signed Circulant Random Feature mapping (SCRF),

If d ∤ D, we can generate t = ⌈D/d⌉ signed circulant Gaussian matrices. By combining these t signed circulant Gaussian matrices in a column, after random feature mapping we take the {1, 2, . . . , d, . . . , (t-1)d, (t-1)d+1, . . . ,D}-th entries of Eq. (9) as the final random feature mapping. In this paper, without any specification, we assume that d | D.

With the structured random matrix projection, we can implement the SCRF efficiently by using FFT.

Lemma 2 (Computational Complexity). The feature mappings of SCRF (Eq.(9)) can be computed in O(D log d) time and O(D) space complexity. To predict an input data is in O(D log d) time complexity.

Proof. For the i-th block (Lemma 1), the corresponding feature representation x(i) can be calculated as follows,

To compute x(i) could proceed efficiently by using FFT and the inverse FFT, which is in O(d log d) time complexity. Changing the signs of rows will be in O(d) time complexity. x(i) +bi will cost O(d) time complexity. Only σi, ci and bi need to be stored, which is in O(3d) space complexity. To sum up, the total time complexity is O(D log d) and the space complexity is O(3D).

Obviously, the predicting function has the form f(x) = 〈w, ΦSCRF(x)〉, which will be computed in O(D log d) time complexity.

3.2 Theoretical Analysis

The goal of this subsection is to analyse the approximation quality of SCRF. Lemma 3 states that the approximation using our proposed SCRF is an unbiased estimate of the Gaussian kernel. Corollary 2 shows that the variance of the approximation using SCRF is the same with that using RKS and much smaller than that using Fastfood. SCRF has smaller concentration error than Fastfood as shows in Theorem 3.

Lemma 3. The expected feature mapping (Eq.(9)) recovers the Gaussian kernel, i.e.,

Proof. Any given row of a circulant Gaussian matrix is a Gaussian vector with distribution N (0, I/σ2 ).

Obviously, Eq. (10) holds.

Remark 1. Any two rows of circulant Gaussian matrix are considerably more correlated, which causes increasing variance of approximating Gaussian kernel. We change signs of the rows randomly to eliminate the correlations and define the signed circulant random matrix. Theoretical analysis states that our proposed random feature mapping with signed circulant matrix projection is unbiased and does not increase the variance of the approximation. Fastfood generates correlated random features while SCRF generates independent random features.

Theorem 2. Assume that P is a signed circulant Gaussian matrix with the associated circulant Gaussian matrix C = Circ(c0, c1, . . . , cd-1), where c0, c1, . . . , cd-1 are randomly drawn i.i.d. from N(0, 1/σ2 ). Let v = x- y and φj(v) = cos([Pv]j), the estimate of the kernel value that comes from the k-th pair of random features for each k ∈ {0, 1, . . . , d - 1}. Then for each k we have

and

Proof. Pk· = σkCk·, k ∈ {0, 1, . . . d - 1}. Sign changes retain Gaussianity so that Pk· is also Gaussian vector with i.i.d entries, where E[Pkj ] = E[σkCkj ] = 0 and

so zk ∼ N(0, ‖v‖2/σ2 ). Therefore, we have

Hence, we can obtain the following equation

For j ̸= k,

Therefore, the correlation coefficient ρjk between zj and zk equals to zero,

Obviously, Eq. (12) holds.

The proof completes.

Corollary 2. Let v = x- y. Then for the feature mapping of SCRF (Eq.(9)), Φ : Rd → RD , obtained by stacking D/d signed circulant Gaussian matrices, we have that

Proof. 〈ΦSCRF(x),ΦSCRF(y)〉 is the average of D/d independent esitmates. Following the proof of Theorem 2, if the two rows Pj· and Pk· are in distinct signed circulant Gaussian matrix blocks, the corresponding correlation coefficient ρjk between zj and zk is also zero. Obviously, this corollary holds.

All of the three approximations are unbiased, but lower variance means tighter concentration error. Corollary 2 states that SCRF approximates kernel with the same variance as RKS and smaller variance than Fastfood. Combining Corollary 2 and Chebyshev's inequality, we can easily obtain the concentration error of SCRF, which will converge quickly as D increases with probability 1 - δ.

Theorem 3 (Concentration Error). For the feature mapping of SCRF (Eq.(9)), Φ : Rd → RD , obtained by stacking D/d signed circulant Gaussian matrices, and any δ ∈ (0, 1), the following inequality holds with probability 1 - δ,

4 Experiments

We implement random feature mappings in R 3.1.1 and conduct experiments on a public SUSE Linux enterprise server 10 SP2 platform with 2.2GHz AMD Opteron Processor 6174 CPU and 48GB RAM. We compare the performance of random feature mappings, including RKS, Fastfood and our SCRF in terms of kernel approximation, efficiency of kernel expansion and generalization performance.

4.1 Kernel Approximation

First, we evaluate the kernel estimates from RKS, Fastfood and SCRF. We uniformly sample l = 100 vectors from [0, 1]16 , and set D = 512 and kernel parameter γ = 0.25 (γ = 1/(2σ2 )). Figure 1(a)-1(c) show kernel estimates from the three methods plotted against exact kernel values respectively. Each point represents a combination of two vectors x, y. The coordinate corresponds to (k(x, y), 〈Φ(x),Φ(y)〉). A perfect mapping would manifest as a narrow 45-degree line. As we can see, both RKS and SCRF perform a little better than Fastfood in terms of kernel estimates, which coincides with the fact that the variance of approximation using SCRF is the same with that using RKS and smaller than that using Fastfood.

Next, we investigate the convergence of kernel approximation using SCRF. We uniformly sample l = 500 vectors from [0, 1]16 and compare RKS, Fastfood and SCRF with the exact kernel values. Figure 1(d) shows the relative kernel approximation error (‖̂ K-K‖F/‖K‖F) w.r.t. D. From Figure 1(d), all the three approaches converge quickly to the exact kernel values as D increases, where both RKS and SCRF converge faster than Fastfood.

4.2 Efficiency of Kernel Expansion

In this subsection, we compare the CPU time of computing random feature mappings of the three approaches. We uniformly sample l = 5000 vectors from [0, 1]d and set D = 8192. As analysed above, both Fastfood and SCRF share O(lD log d) time complexity while RKS requires O(ldD) time. Obviously, the running time of both Fastfood and SCRF are less dependent on d, a very promising property since random feature mapping often contributes a significant computational cost in predicting.

Table 1 shows the CPU time in seconds of the three methods with speedup and memory improvement. The running time of both Fastfood and SCRF is almost independent on d. It costs more time in generating random parameters for Fastfood than SCRF, thus SCRF runs a little faster than Fastfood. Both Fastfood and SCRF save much storage compared with RKS. However, SCRF uses around 1.5x less storage than Fastfood. This is because Fastfood needs to store 5D parameters while SCRF only stores 3D parameters.

For higher-dimensional problems, we need to increase D to boost the accuracy, D = O(d). Figure 2 shows that the CPU time of computing D = d random feature mapping with RKS is quadratic with d, a bottleneck of kernel methods on high-dimensional datasets. However, both Fastfood and SCRF scale well in this case. Especially, SCRF runs faster and uses less memory than Fastfood.

4.3 Generalization Performance

In this subsection, we compare RKS, Fastfood and SCRF with non-linear SVMs on 9 well-known classification benchmark datasets of size ranging from 1,000 to 60,000 with dimensionality ranging from 22 to 3,072. For mushrooms, we select 4,062 training data randomly and the left as test data. We use LIBSVM [Chang and Lin, 2011] for non-linear kernels and LIBLINEAR [Fan et al., 2008] for random feature mappings for classification task. All averages and standard deviations are over 5 runs of the algorithms except on cifar10. We select the kernel parameter γ and regularization coefficient C of Gaussian kernel SVM by using 5-fold cross validation with LIBSVM.

Figure 3 depicts the test accuracies of RKS, Fastfood and SCRF with LIBLINEAR w.r.t. D on dna and ijcnn1. Their generalization performances are not obviously distinguishable as D increases. Table 2 shows the results of the comparison. There is virtually no difference among them in terms of test accuracy except on cifar10. All of the results of the three approaches on cifar10, which come from one random trial respectively, are worse than LIBSVM. This is because the selected γ and C by LIBSVM may not be optimal for the other three approaches. Except on mnist and cifar10, Fastfood runs slower than RKS and SCRF in both training and predicting. This is because Fastfood needs to format d = 2k , k ∈ N, through padding the vector with zeros, which demands more time for preprocessing input data and random feature mapping. In addition, if d is not so large that there will be no speedup of computing random feature mappings. For larger dimensional dataset mnist/cifar10, Fastfood can obtain efficiency gains. In general, SCRF is more time and space efficient than both Fastfood and RKS.

5 Conclusion

In this paper, we have proposed a random feature mapping method for approximating Gaussian kernel using signed circulant matrix projection that makes the approximation unbiased and have lower variance. The adoption of circulant matrix projection guarantees a quasilinear random feature mapping and promotes scalable and practical kernel methods for large scale machine learning.

Acknowledgement

This work was supported in part by National Natural Foundation of China (No. 61170019), 973 Program (2013CB329304), National Natural Foundation of China (No. 61222210) and Key Program of National Natural Science Foundation of China (No. 61432011).