Sharp analysis of low-rank kernel matrix approximations

Introduction

Kernel methods, such as the support vector machine or kernel ridge regression, are now widely used in many areas of science and engineering, such as computer vision or bioinformatics (see, e.g., Schölkopf et al., 2004; Zhang et al., 2007). Their main attractive features are that (1) they allow non-linear predictions through the same algorithms than for linear predictions, owing to the kernel trick; (2) they allow the separation of the representation problem (designing good kernels for non-vectorial data) and the algorithmic/theoretical problems (given a kernel, how to design, run efficiently and analyze estimation algorithms). Moreover, (3) their applicability goes beyond supervised learning problems, through the kernelization of classical unsupervised learning techniques such as principal component analysis or K-means. Finally, (4) probabilistic Bayesian interpretations through Gaussian processes allow their simple use within larger probabilistic models. For more details, see, e.g., Rasmussen and Williams (2006); Schölkopf and Smola (2002); Shawe-Taylor and Cristianini (2004).

However, kernel methods typically suffer from at least quadratic running-time complexity in the number of observations $n$, as this is the complexity of computing the kernel matrix. In large-scale settings where $n$ may be large, this is usually not acceptable. In these situations where plain kernel methods cannot be run, practitioners would commonly (a) turn to methods such as boosting, decision trees or random forests, which have both good running time complexity and predictive performance. However, these methods are typically run on data coming as vectors and usually put a strong emphasis on a sequence of decisions based on single variables. Another common solution is (b) to stop using infinite-dimensional kernels and restrict the kernels to be essentially linear kernels (i.e., by choosing an explicit representation of the data whose size is independent of the number of observations) where the non-parametric kernel machinery (of adapting the complexity of the underlying predictor to the size of the dataset) is lost, and the methods may then underfit.

In this paper, we consider the traditional kernel set-up for supervised learning, where the input data are only known through (portions of) the kernel matrix. The main question we try to tackle is the following: Is it possible to run supervised learning methods with positive-definite kernels in time which is subquadratic in the number of observations without losing predictive performance? Of course, if adaptation is desired, linear complexity seems impossible, and therefore we should expect (hopefully slightly) super-linear algorithms. Statistically, a quantity that characterizes the non-parametric nature of kernel method is the degrees of freedom, which play the role of an implicit number of parameters and which we define and review in Section 4.1. This quantity allows to go beyond worst-case analyses which are common in statistical learning theory: our generalization bounds will then depend on problem-dependent quantities which may not be known at training time, but that characterize finely the behavior on any given problem instance, and not only for the worst case over a large class of problems. In this paper, we try to tackle the following specific question: Do the degrees of freedom play a role in the computational properties of kernel methods?

An important feature of kernel matrices is that they are positive-semidefinite, and thus they may well be approximated from a random subset of $p$ of their columns, in running-time complexity $O(p^{2}n)$ and with a computable bound on the error (see details in Section 3). This appears through different formulations within numerical linear algebra or machine learning, e.g., Nyström method (Williams and Seeger, 2001), sparse greedy approximations (Smola and Schölkopf, 2000), incomplete Cholesky decomposition (Fine and Scheinberg, 2001; Bach and Jordan, 2005), Gram-Schmidt orthonormalization (Shawe-Taylor and Cristianini, 2004) or CUR matrix decompositions (Mahoney and Drineas, 2009). It has been thoroughly analyzed in contexts where the goal is kernel matrix approximation or approximate eigenvalue decomposition (see, e.g., Boutsidis et al., 2009; Mahoney and Drineas, 2009; Kumar et al., 2012; Gittens, 2011). Such bounds have also been subsequently used to characterize the approximation of predictions made from these low-rank decompositions (Cortes et al., 2010; Jin et al., 2001), but these two-stage analyses do not lead to guarantees that reflect the good observed practical behavior. In this paper, our analysis aims at answering explicitly the simple question: how big should $p$ be to incur no loss of predictive performance compared to the full kernel matrix? The key insight of this paper is not to try to approximate the kernel matrix well, but to predict well from the approximation. This requires a sharper analysis of the approximation properties of the column sampling approach.

In the fixed design least-squares regression setting, we show in Section 4.2 that the rank $p$ can be chosen to be linear in the degrees of freedom associated with the problem, a quantity which is classically used in the statistical analysis of such methods. Note that our results hold for any problem instance, and not only in a worst-case regime.

We present in Section 4.4 simple algorithms that have sub-quadratic running time complexity, and, for the square loss, provably exhibit the same predictive performance as classical algorithms than run in quadratic time (or more).

We provide in Section 4.3 explicit examples of optimal values of the regularization parameters and the resulting degrees of freedom, as functions of the decay of the eigenvalues of the kernel matrix, shedding some light in the joint computational/statistical trade-offs for choosing a good kernel. In particular, we show that with kernels with fast spectrum decays (such as the Gaussian kernel), computational limitations may prevent exploring the relevant portions of the regularization paths, leading to underfitting.

Supervised learning with positive-definite kernels

In this section, we present the problem we try to solve, as well as several areas of the machine learning and statistics literatures our method relates to.

Let $(x_{i},y_{i})$, $i=1,\dots,n$, be $n$ pairs of points in $\mathcal{X}\times\mathcal{Y}$, where $\mathcal{X}$ is the input space, and $\mathcal{Y}$ is the set of outputs/labels. In this paper, we consider the problem of minimizing

First, using the representer theorem, the unique solution $f$ may be found as $f=\sum_{i=1}^{n}\alpha_{i}\phi(x_{i})$ (see, e.g., Wahba, 1990; Schölkopf and Smola, 2002; Shawe-Taylor and Cristianini, 2004). Thus, by replacing the expression of $f$ in Eq. (1), $\alpha$ is a solution of the following optimization problem:

Second, for convex losses only, an equivalent dual problem is classically obtained as (see proof in Appendix A):

2 Related work

In order to solve Eq. (1), algorithms typically consider a primal or a dual approach. Solving Eq. (2), i.e., the primal formulation after application of the representer theorem, is typically inefficient because the problem is ill-conditionedThe objective function in Eq. (2) is a function of $K^{1/2}\alpha$, with a kernel matrix $K$ which is often ill-conditioned, usually leading to ill-conditioning of the original problem (Chapelle, 2007). and thus second-order algorithms are typically used (Chapelle, 2007). Alternatively, $K$ is represented explicitly as $K=\Phi\Phi^{\top}$ and a change of variable $w=\Phi^{\top}\alpha$ is considered (note that when the kernel $k$ is linear, $\Phi$ is simply the design matrix, and we are solving directly a linear supervised learning problem). Then, the classical battery of convex optimization algorithms may be used, such as gradient descent, stochastic gradient descent (Shalev-Shwartz et al., 2007) or cutting-planes (Joachims et al., 2009). However, in a kernel setting where a small matrix $\Phi$ (i.e., with few columns) is not known a priori, then they all exhibit at least quadratic complexity in $n$, as the full kernel matrix is used.

The dual problem in Eq. (3) is usually better-behaved (it has a better condition number) (Chapelle, 2007), and algorithms such as coordinate descent and its variants such as sequential minimal optimization may be used (Platt, 1999). Again, in general, the full kernel matrix is needed.

Some algorithms operate online and do not need to compute the full kernel matrix, such as the “forgetron” (Dekel et al., 2005), the “projectron” (Orabona et al., 2008), BGSD (Wang et al., 2012), or LASVM (Bordes et al., 2005), with typically a fixed computational budget and good practical performance. They often come with theoretical approximation guarantees, which are either data-dependent or based on worst-case analysis; however, these do not characterize the required rank which is needed to achieve the same accuracy than the problem with a full kernel matrix. In fact, one of the main motivations for this work is to derive precise bounds for reduced-set stochastic gradient algorithms for supervised kernel problems.

Analysis of column sampling approximation.

Given a positive semi-definite matrix $K$ of size $n$, many methods exist for approximating it with a low-rank (typically also positive semidefinite) matrix $L$. While the optimal approximation is obtained from the eigenvalue decomposition, it is not computationally efficient as it has complexity at least quadratic in $n$ (since it requires the knowledge of $K$). In order to achieve linear complexity in $n$, approximations from subsets of columns are considered and appear under many names: Nyström method (Williams and Seeger, 2001), sparse greedy approximations (Smola and Schölkopf, 2000), incomplete Cholesky decomposition (Fine and Scheinberg, 2001), Gram-Schmidt orthonormalization or CUR matrix decompositions (Mahoney and Drineas, 2009). Note that reduced-set methods (see, e.g., Keerthi et al., 2006) typically consider using a subset of columns after the predictor has been estimated. These low-rank methods are described in Section 3 and have running time complexity $O(p^{2}n)$ for an approximation of rank $p$. Note that they may also be used in a Bayesian setting with Gaussian processes (see, e.g., Lawrence et al., 2002).

Column sampling has been analyzed a lot (Mahoney and Drineas, 2009; Cortes et al., 2010; Kumar et al., 2012; Talwalkar and Rostamizadeh, 2010; Gittens, 2011); however, typically the analysis provides a high-probability bound on the error $\|K-L\|$ for an appropriate norm (typically operator, Frobenius or trace norm), but this is too pessimistic and does not really match with good practical performance (see empirical evidence in Figure 1). Some works do consider prediction guarantees (Cortes et al., 2010; Jin et al., 2001), but as shown in Section 4.2, these are not sufficient to reach sharp results depending on the degrees of freedom. Moreover, many analyses consider situations where the matrix $K$ is close to low-rank, which is not the case with kernel matrices. In this paper, the control of $K-L$ is more precise and adapted to the use of $K$ within a supervised learning method.

Randomized dimension reduction.

The method presented in this paper, which considers random columns from the original kernel matrix, is also related to random projection techniques used for linear prediction problems (Mahoney, 2011; Maillard and Munos, 2009). These techniques are not kernel methods per se, as they require the knowledge of a matrix square root $\Phi$ (such that $K=\Phi\Phi^{\top}$), which leads to complexity greater than quadratic. For certain kernel functions that can be explicitly expressed as an expectation of dot-products in low-dimensional spaces, similar randomized dimensionality reduction may be performed (Rahimi and Recht, 2007). Note however that the dimension reduction is then independent of the particular distributions of the input data points, while the column sampling approach is; see Yang et al. (2012) for more discussion.

Theoretical analysis of predictive performance of kernel methods.

In order to assess the required precision in approximating the kernel matrix, it is key to understand the typical predictive performance of kernel methods. For the square loss, this is classically obtained from a bias-variance decomposition of this performance (see Section 4). A key quantity is the degrees of freedom, which play the role of an implicit number of parameters and is applicable to many non-parametric estimation methods which consists in “smoothing” the response vector by a linear operator (see, e.g., Wahba, 1990; Hastie and Tibshirani, 1990; Gu, 2002; Caponnetto and De Vito, 2007; Hsu et al., 2011). See precise definitions in Section 4.1.

Approximation from subset of columns

Given a random subset $I$ of $V=\{1,\dots,n\}$ of cardinality $p$, we simply consider the approximation of the kernel matrix $K$ from the knowledge of $K(V,I)$ (the columns of $K$ indexed by $I$), by the matrix

where $M^{\dagger}$ denotes the pseudo-inverse of $M$. As shown by Bach and Jordan (2005), $L$ is the only symmetric matrix with column space spanned by the columns of $K(V,I)$, and such that $L(V,I)=K(V,I)$. Alternatively, given that $K$ is the matrix of dot-products of points in a Hilbert space, it may be seen as the kernel matrix of the orthogonal projections of all points onto the affine subspace spanned by the points indexed by $I$ (Mahoney, 2011).

Such a feature map may be efficiently obtained in running time $O(p^{2}n)$ using incomplete Cholesky decomposition—often interpreted as partial Gram-Schmidt orthonormalization, with the possibility of having an explicit online bound on the trace norm of the approximation error (see, e.g., Shawe-Taylor and Cristianini, 2004).

Pivoting vs. random sampling.

While selecting a random subset is computationally efficient, it may not lead to the best performance. For the task of approximating the kernel matrix, algorithms such as the incomplete Cholesky decomposition with pivoting, provide an approximate greedy algorithm with the same complexity than random subsampling (Smola and Schölkopf, 2000; Fine and Scheinberg, 2001).

In Section 5, we provide comparisons between the two approaches, showing the potential advantage of the greedy method over random subsampling. However, the analysis of such algorithms is harder, and, to the best of our knowledge, still remains an open problem.

Fixed design analysis for least-square regression (ridge regression)

Following classical results from the statistics literature (see, e.g., Wahba, 1990; Hastie and Tibshirani, 1990; Gu, 2002), we obtain the following expected prediction error:

which may be classically decomposed in two terms:

Note that the bias term is a matrix-decreasing function of $K/\lambda$ (and thus an increasing function of $\lambda$), while the variance term is a matrix-increasing function of $K/\lambda$ and the noise covariance matrix $C$.

Degrees of freedom.

Note that an assumption which is usually made is $C=\sigma^{2}I$; the variance term then takes the form $\frac{\sigma^{2}}{n}\mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}$ and $\mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}$ is referred to as the degrees of freedom (Wahba, 1990; Hastie and Tibshirani, 1990; Gu, 2002; Hsu et al., 2011) (note that an alternative definition is often used, i.e., $\mathop{\rm tr}K(K+n\lambda I)^{-1}$, and that as shown in Appendix C, they often behave similarly). In ordinary least-squares estimation from $d$ variables, the variance term is equal to $\sigma^{2}d/n$, and thus the degrees of freedom play the role of an implicit number of parameters. In this paper, we show that a proxy to this statistical quantity also plays a role in optimization: the number of columns needed to approximate the kernel matrix precisely enough to incur no loss of performance is linear in the degrees of freedom.

More precisely, we define the maximal marginal degrees of freedom $d$ as

We have \mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}\leqslant\mathop{\rm tr}K(K+n\lambda I)^{-1}=\big{\|}\mathop{\rm diag}\big{(}K(K+n\lambda I)^{-1}\big{)}\big{\|}_{1}\leqslant d, and thus $d$ provides an upper-bound on the regular degrees of freedom. It may be significantly larger in situations where there may be outliers and the vector \mathop{\rm diag}\big{(}K(K+n\lambda I)^{-1}\big{)} is far from uniform—precise results are out of the scope of this paper. Moreover, the diagonal elements of $K(K+n\lambda I)^{-1}$ are related to statistical leverage scores introduced for best-rank approximations (Mahoney and Drineas, 2009); it would be interesting to see if this link could lead to non-uniform sampling schemes with better behavior.

In Section 4.3, we study in detail how the degrees of freedom vary as a function of $\lambda$ and $n$: in order to minimize predictive performance, the best choice of $\lambda$ depends on $n$ (as a decreasing function), typically smaller than a constant times $1/\sqrt{n}$, and the degrees of freedom typically grow as a slow function of $n$, reflecting the non-parametric nature of kernel methods.

2 Predictive performance of column sampling

We consider sampling $p$ columns (without replacement) from the original $n$ columns. We consider the column sampling approximation defined in Eq. (4) and provide sufficient conditions (a lower-bound on $p$) to obtain the same predictive performance than with the full kernel matrix.

Proof sketch. The proof relies on approximating the expected error directly, and not through bounding the error $\|K\!-\!L\|$. This is done by (a) considering a regularized version of $L$, i.e., L_{\gamma}=K(V,I)\big{[}K(I,I)+p\gamma I\big{]}^{-1}K(I,V), (b) using a Bernstein inequality for an appropriately rescaled covariance matrix and (c) using monotonicity arguments to obtain the required bound. See more details in Appendix B. \BlackBox

The bound in Eq. (7) provides a relative approximation guarantee: the predictions $\hat{z}_{L}$ are shown to perform as well as $\hat{z}_{K}$ (no kernel matrix approximation). Small values of $\delta$ impose no loss of performance, while $\delta=1/4$ impose that the prediction errors have a similar behavior (up to a factor of $2$). Note that relative bounds may be more easily obtained by eigenvalue thresholding of the kernel matrix, i.e., through replacing soft-shrinkage by hard-thresholding (Blanchard et al., 2004; Dhillon et al., 2011). However, these bounds do not allow the proportionality constant to go arbitrarily close to one, and they depend on the full knowledge of the kernel matrix.

Lower bounds.

The lower bound for the rank $p$ in Eq. (6) shows that the maximal marginal degrees of freedom provides a quantity which, up to logarithmic terms, is sufficient to scale with, in order to incur no loss of prediction performance. Note that the previous result also allows the derivation of an approximation guarantee $\delta$ given a rank $p$, by inverting Eq. (6). Moreover, Theorem 4.1 provides a sufficient lower-bound for the required rank $p$. Deriving precise necessary lower-bounds is outside the scope of this paper. However, given that with a reduced space of $p$ dimensions, we can achieve a prediction error of $O(p/n)$ from ordinary least-squares, we should expect $p$ to be larger than the known minimax rates of estimation for the problem at hand (Johnstone, 1994; Caponnetto and De Vito, 2007; Steinwart et al., 2009). In Section 4.3, we show that in some situations, it turns out that $d$ is of the order of the minimax rate; therefore, we could expect that in certain settings, $d$ is also a necessary lower-bound on $p$ (up to constants and logarithms).

High-probability results.

Avoiding terms in 1/λ1𝜆1/\lambda.

Theorem 4.1 focuses on average predictive performance; this is different from achieving a good approximation of the kernel matrix (Mahoney and Drineas, 2009). Previous work (Cortes et al., 2010; Jin et al., 2001) considers explicitly the use of kernel matrix approximation bounds within classifiers or regressors, but obtains bounds that involve multiplicative terms of the form $1/\lambda$ or $1/\lambda^{2}$, which, as we show in Section 4.3, would grow as $n$ grows. More precisely, the bound from Cortes et al. (2010, Eq. (5)) has a term of the form $1/\lambda_{\min}(K+\lambda I)^{2}$; however, for the non-parametric problems we are considering in this paper, the lowest eigenvalue of $K$ is often below machine precision and hence the bound behaves as $1/\lambda^{2}$. Moreover, the subsequent bound of Cortes et al. (2010, Theorem 2) has a term of the form $1/(\lambda^{2}p)$ which can only be small if $p$ is larger than $1/\lambda^{2}$, which, according to our analysis in Section 4.3, is typically larger than $n$ since $\lambda$ should typically decrease at least as $1/\sqrt{n}$ (however, note that their bound has a stronger nature than the one in Theorem 4.1, as it states a guarantee in high probability).

Our proof technique, that focuses directly on prediction performance and side-steps the explicit approximation of the kernel matrix, avoids these terms in $1/\lambda$, and, beyond the dependence on $\lambda$ through the degrees of freedom (which we cannot avoid), our dependence is only logarithmic in $\lambda$ (see details in the proof in Appendix B).

Instance-based guarantees.

Theorem 4.1 shows that in the specific instance that we are faced with, we do not lose any average predictive performance. As opposed to Jin et al. (2001), the bound is not on the worst-case predictive performance (obtained from optimizing over $\lambda$, and with worst-case analysis over $K$), but for given $\lambda$ and $K$ (however, the bound of Jin et al. (2001) is a high-probability result while ours is only in expectation). Moreover, even in this worst-case regime, Jin et al. (2001) state that for polynomial decays of the eigenvalues of $K$ such that the optimal prediction performance is of the form $O(n^{-1}n^{1/(\gamma+1)})$ (and for binary classification rather than regression), the rank to achieve this optimal prediction is $p=n^{2\gamma/(\gamma^{2}-1)}$, which may be larger than $n$ and is significantly higher than $n^{1/(\gamma+1)}$ (which corresponds to our result since the degrees of freedom are then equal to $n^{1/(\gamma+1)}$).

Link with eigenvalues.

In the existing analysis of sampling techniques for kernel methods, another source of inefficiency which makes our result sharper is the proof technique for bounding $\|K-L\|$. Indeed, most analyses use a linear algebra lemma from Mahoney and Drineas (2009); Boutsidis et al. (2009), that relies on the $(p+1)$-th eigenvalue to be small; hence it is adapted to matrices with sharp eigenvalue decrease, which is not the case for kernel matrices (see an illustrative example in Figure 1). We provide a new proof technique based on regularizing the column sampling approximation and optimizing the extra regularization parameter using a monotonicity argument.

Additional regularization effect.

In our experiments, we have noticed that the low-rank approximation may have an additional regularizing effect leading to a better prediction performance than with the full kernel matrix.

Beyond square loss.

The notion of degrees of freedom can be extended to smooth losses like the logistic loss. However, the simple bias/variance decomposition only holds asymptotically, forcing a control of the two terms, which would lead to significant added complexity.

Beyond fixed design.

In order to extend our analysis to random design settings, we would need to additionally control the deviation between covariance operators and empirical covariance operators, with quantities like the degrees of freedom that depend on the decay of non-zero eigenvalues and not on their number. This could be done using tools from Hsu et al. (2011, 2012).

3 Optimal choice of the regularization parameter

As seen in Sections 4.1 and 4.2, the computational and statistical properties of kernel ridge regression depend heavily on the choice of the regularization parameter $\lambda$ as $n$ increases, which we now tackle.

For simplicity, in this section, we assume that the noise variables $\varepsilon$ are i.i.d. (i.e., $C=\sigma^{2}I$). Our goal is to study simplified situations, where we can derive explicit formulas for the bias, the variance, and the optimal regularization parameter. Throughout this section, we will consider specific decays of certain sequences, which we characterize with the notation $u_{n}=\Theta(v_{n})$, which means that there exist strictly positive constants $A$ and $B$ such that $Au_{n}\leqslant v_{n}\leqslant Bv_{n}$ for all $n$.

We assume that the kernel matrix $K$ has eigenvalues of the form $\Theta(n\mu_{i})$, $i=1,\dots,n$, for some summable sequence $(\mu_{i})$—so that $\mathop{\rm tr}K=\Theta(n)$, and that the coordinates of $z$ on the eigenbasis of $K$ have the asymptotic behavior $\Theta(\sqrt{n\nu_{i}})$ for a summable sequence $(\nu_{i})$—so that $\frac{1}{n}z^{\top}z=\Theta(1)$. In Table 1, we provide asymptotic equivalents of all quantities for several pairs of sequences $(\mu_{i})$ and $(\nu_{i})$ (see proofs in Appendix C), with polynomial or exponential decays.

Note that for decays of $\nu_{i}$ which are polynomial, i.e., $\nu_{i}=O(i^{-2\delta})$, then the best possible prediction performance is known to be $O(n^{1/2\delta-1})$ (Johnstone, 1994; Caponnetto and De Vito, 2007) and is achieved if the RKHS is large enough (lines 2 and 4 in Table 1). For exponential decay, the best performance is $O(\log n/n)$. See also Steinwart et al. (2009).

The RKHS is too large, for lines 1 and 3 in Table 1: the eigenvalues of $K$, which depend linearly on $\mu_{i}$, do not decay fast enough. In other words, the functions in the RKHS are not smooth enough. In this situation, the prediction performance is suboptimal (do not attain the best possible rate).

The RKHS is too small, for lines 2, 4, and 6 in Table 1: the eigenvalues of $K$ decay fast enough to get an optimal prediction performance. In other words, the functions in the RKHS are potentially smoother than what is necessary. In this situation however, the required value of $\lambda$ may be very small (much smaller than $O(n^{-1})$), leading to potentially harder optimization problems (since the condition number that depends on $1/\lambda$ may be very large).

There is thus a computational/statistical trade-off: if the RKHS is chosen too large, then the prediction performance is suboptimal (i.e., even with the best possible regularization parameter, the resulting error is not optimal); if the RKHS is chosen too small, the prediction performance could be optimal, but the optimization problems are harder, and sometimes cannot be solved with the classical precision of numerical techniques (see examples of such behavior in Section 5). Indeed, for least-squares regression, where a positive semi-definite linear system has to be solved, its condition number is proportional to $1/\lambda$ and for the best choice of $\lambda$ in line 4 of Table 1, it grows exponentially fast in $n$ (if $\lambda$ is chosen larger, then the bias term will lead to sub-optimal behavior).

4 Optimization algorithms with column sampling

Given a rank $p$ and a regularization parameter $\lambda$, we consider the following algorithm to solve Eq. (1) for twice differentiable convex losses:

Select at random $p$ columns of $K$ (without replacement).

The complexity of step 2 is already $O(p^{2}n)$, therefore using faster techniques for step 3 (e.g., accelerated gradient descent) does not change the overall complexity, which is thus $O(p^{2}n)$. Moreover, since we use a second-order method for step 3, we are robust to ill-conditioning and in particular to small values of $\lambda$ (though not below machine precision as seen in Section 5). This is not the case for algorithms that relies on the strong convexity of the objective function, whose convergence is much slower when $\lambda$ is small (as seen in Section 4.3, when $n$ grows, the optimal value of $\lambda$ can decay very rapidly, making these traditional methods non robust).

According to Theorem 4.1, at least for the square loss, the dimension $p$ may be chosen to be linear in the maximal marginal degrees of freedom $d$ defined in Eq. (5); it is in practice close to the traditional degrees of freedom $d_{ave}$, which, as illustrated in Table 1, is typically smaller than $n^{1/2}$. Therefore, if $p$ is properly chosen, the complexity is subquadratic. Given $\lambda$, $d$ (and thus $p$) can be estimated from a low-rank approximation of $K$. However, our current analysis assumes that $\lambda$ is given. Selecting the rank $p$ and the regularization parameter $\lambda$ in a data-driven way would make the prediction method more robust, but this would require extra assumptions (see, e.g., Arlot and Bach, 2009, and references therein).

Simulations

In order to study various behaviors of the regularization parameters $\lambda$ and the degrees of freedom $d$, we consider periodic smoothing splines on $$and points$x_{1},\dots,x_{n}$uniformly spread over$$, either deterministically or randomly. In order to generate problems with given sequences$(\mu_{i})$and$(\nu_{i})$, it suffices to choose$k(x,y)=\sum_{i=1}^{\infty}2\mu_{i}\cos 2i\pi(x-y)$, and a function$f(x)=\sum_{i=1}^{\infty}2\nu_{i}^{1/2}\cos 2i\pi x$. For$\mu_{i}=i^{-2\beta}$, we have$k(x,y)=\frac{1}{(2\beta)!}B_{2\beta}(x-y-\lfloor x-y\rfloor)$, where$B_{2\beta}$is the$(2\beta)$-th Bernoulli polynomial (see details in Appendix D).

Optimal values of λ𝜆\lambda.

In a first experiment, we illustrate the results from Section 4.3, and compute in Figure 2 the best value of the regularization parameter (left) and the obtained predictive performance (middle), for a problem with $\nu_{i}=i^{-2\delta}$ for $\delta=8$, and for which we considered several kernels, for which $\mu_{i}=i^{-2\beta}$, for $\beta=1$, $\beta=4$ and $\beta=8$.

For $\beta=1$, the rate of convergence of $n^{1/(4\beta+1)-1}$ happens to be achieved (line 1 in Table 1), with a certain asymptotic decay of the regularization parameter, and it is slower than $n^{1/(2\delta)-1}$.

For $\beta=4$, the optimal rate of $n^{1/(2\delta)-1}$ is achieved (line 2 in Table 1), as expected.

For $\beta=8$, the rate of convergence should be $n^{1/(2\delta)-1}$ (line 2 in Table 1), however, as seen in the left plot, the regularization parameter saturates as $n$ grows at the machine precision, leading, because of numerical errors, to worse prediction performance. The problem is so ill-conditioned that the matrix inversion cannot be algorithmically robust enough.

Performance of low-rank approximations.

In this series of experiments, we compute the rank $p$ which is necessary to achieve a predictive performance at most $1\%$ worse than with $p=n$, and computeNote that in practice, computing the degrees of freedom exactly requires to know the full matrix. However, it could also be approximated efficiently, following for example Drineas et al. (2012). the ratio with the marginal degrees of freedom d=n\big{\|}\mathop{\rm diag}\big{(}K(K+n\lambda I)^{-1}\big{)}\big{\|}_{\infty} and the traditional degrees of freedom $d_{ave}=\mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}$. In the right plot of Figure 2, we consider data randomly distributed in $$with the same kernels and functions than above, while in Figure 3, we considered three of the pumadyn datasets from the UCI machine learning repository (there we compute the classical generalization performance on unseen data points and estimate$\lambda$ by cross-validation). For further experimental evaluations, see Cortes et al. (2010); Talwalkar and Rostamizadeh (2010).

On all datasets, the ratios stay relatively close to one, illustrating the results from Theorem 4.1. Moreover, the two different versions $d$ and $d_{ave}$ of degrees of freedom are within a factor of 2. Moreover, using pivoting to select the columns does not change significantly the results, but may sometimes reduce the number of required columns by a constant factor. Note that the sudden increase (of magnitude less than 2 in the middle and right plots of Figure 3) is due to the chosen criterion (ratio of the sufficient rank to obtain 1% worse predictive performance), which may be unstable.

Conclusion

In this paper, we have provided an analysis of column sampling for kernel least-squares regression that shows that the rank may be chosen proportional to properly defined degrees of freedom of the problem, showing that the statistical quantity characterizing prediction performance also plays a computational role. The current analysis could be extended in various ways: First, other column sampling schemes beyond uniform, such as presented by Boutsidis et al. (2009); Kumar et al. (2012), could be considered with potentially better behavior. Moreover, while we have focused on a fixed design, it is of clear interest to extend our results to random design settings using tools from Hsu et al. (2011). The analysis may also be extended to other losses than the square loss, such as the logistic loss, using self-concordant analysis (Bach, 2010) or the hinge loss, using eigenvalue-based criteria (Blanchard et al., 2008) or tighter approaches to sample complexity analysis (Sabato et al., 2010). Finally, in this paper, we have considered a batch setting and extending our results to online settings, in the line of Tarrès and Yao (2011), is of significant practical and theoretical interest.

This work was supported by the European Research Council (SIERRA Project). The author would like to thank the reviewers for suggestions that have helped improve the clarity of the paper.

Appendix A Duality for kernel supervised learning

In this section, we review classical duality results for kernel-based supervised learning, which extends the dual problem of the support vector machine to all losses. For more details and examples, see Rifkin and Lippert (2007). We consider the following problem, where $\mathcal{F}$ is an RKHS with feature map $\phi:\mathcal{X}\to\mathcal{F}$:

which may be rewritten with the feature map $\phi$ as:

Minimizing with respect to $(f,u)$, we get $f=\sum_{i=1}^{n}\alpha_{i}\phi(x_{i})$ and the dual problem:

Appendix B Proof of Theorem 4.1

We first prove a lemma that provides a Bernstein-type inequality for subsampled covariance matrices. The proof follows Tropp (2011, 2012) and Gittens (2011).

where $\Xi$ is obtained by sampling independently $p$ rows with replacement, i.e., is equal to

We can then apply the matrix Bernstein inequality of Tropp (2012, Theorem 6.1) to obtain the probability bound:

B.2 Proof of Theorem 4.1

We consider the regularized low-rank approximation $L_{\gamma}=\Phi N_{\gamma}\Phi^{\top}$, with

(obtained using the matrix inversion lemma). We have $L=L_{0}$ but we will consider $L_{\gamma}$ for $\gamma>0$ to obtain a bound for $\gamma=0$, using a monotonicity argument.

The function $\gamma\mapsto N_{\gamma}$ is matrix-non-increasing (i.e., if $\gamma\geqslant\gamma^{\prime}$, then $N_{\gamma}\preccurlyeq N_{\gamma^{\prime}}$). Therefore, we have $0\preccurlyeq N_{\gamma}\preccurlyeq N_{0}\preccurlyeq I$. Since the variance term {\rm variance}(L_{\gamma})=\frac{1}{n}\mathop{\rm tr}C\big{[}\Phi N_{\gamma}\Phi^{\top}(\Phi N_{\gamma}\Phi^{\top}+n\lambda I)^{-1}\big{]}^{2} is non-decreasing in $N_{\gamma}$, this implies that the variance term with $N_{\gamma}$ is smaller than the one with $N_{0}$ and then less then the one with $N_{\gamma}$ replaced by $I$ (which corresponds to the variance term without any approximation). For the bias term we have:

which is a non-decreasing function of $\gamma$. Therefore, if we prove an upper-bound on the bias term for any $\gamma>0$, we have a bound for $\gamma=0$. This requires lower-bounding $N_{\gamma}$.

Thus, in order to obtain a lower-bound on $N_{\gamma}$, it suffices to have an upper-bound of the form

Assume $\frac{\gamma/\lambda}{1-t}\leqslant 1$. We then have, using the previous inequality:

Probabilistic control.

Using the bound from Eq. (11), we get, given $\delta\in(0,1)$, $t=1/2$, and $\gamma=\frac{\lambda\delta}{4}$ (which satisfy $\gamma\leqslant\lambda$ and $\frac{\gamma/\lambda}{1-t}=\delta/2\leqslant 1$),

Thus, if p\geqslant\big{(}\frac{32d}{\delta}+2\big{)}\log\frac{nR^{2}}{\delta\lambda}, we obtain that $B\leqslant 1+4\delta$.

We could improve the bound by expliciting the reduction of the variance term.

In some situations, the prediction performance for the approximated version may in fact be smaller than the non-approximated version.

The proof technique relies on a high-probability bound with respect to the sampling of columns. More precisely, from Eq. (11) and Eq. (12), with $t=1/2$ and $\gamma=\frac{\lambda\delta}{4}$, we get:

Appendix C Asymptotics of bias and variance terms

In this appendix, we consider various decays of eigenvalues $n\mu_{i}$ of $K$ and components $\sqrt{n\nu_{i}}$ (in magnitude) of the signal $z$ to estimate. We follow the reasoning of Harchaoui et al. (2008), i.e., replacing sums by integrals. Given our assumptions, we have:

For all cases we need to consider, for simplicity, we only provide an upper-bound for $\mu_{i}$ exactly equal to its asymptotic equivalent. Considering lower-bounds and a constant times the asymptotic equivalent may be done in a similar way.

We consider the only two possible cases (the variance term only depends on $(\mu_{i})$). Moreover we show that the two traditional definitions of the degrees of freedom, $\mathop{\rm tr}K(K+n\lambda I)^{-1}$ and $\mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}$, have the same asymptotically equivalents.

The renormalized variance term is less than

With the same reasoning, we have $\mathop{\rm tr}K(K+n\lambda I)^{-1}\leqslant\int_{0}^{\lambda n^{2\beta}}\frac{1}{(1+u)}{\lambda}^{-1/2\beta}u^{1/2\beta-1}\frac{1}{2\beta}du=O({\lambda}^{-1/2\beta}).$

The renormalized variance term is less than

We the same technique, we get bounds on $\mathop{\rm tr}K(K+n\lambda I)^{-1}$ in the same way we just did for $\mathop{\rm tr}K^{2}(K+n\lambda I)^{-2}$.

C.2 Bias terms

The bias terms depend on both $(\mu_{i})$ and $(\nu_{i})$ and we consider all combinations.

If $2\delta-4\beta>1$, then we have an upper bound of $2n\lambda^{2}\int_{1}^{\infty}t^{4\beta-2\delta}dt=O(n\lambda^{2})$, because the integral is finite.

If $2\delta-4\beta<1$, then we can further bound Eq. (14) as

because the integral is finite (due to the assumptions made on $\beta$ and $\delta$).

because the integral is finite and uniformly bounded in $\lambda$.

Mixed decays.

For $\mu_{i}$ with polynomial decays and $\nu_{i}$ with exponential decays, we are in a situation where $\nu_{i}$ is decaying fast enough (faster than $i^{-2\delta}$ for any $\delta>1/2$) so that, given previous results, the bias is $n\lambda^{2}$.

The only remaining result to show is $\mu_{i}=e^{-\rho i}$ and $\nu_{i}=i^{-2\delta}$, $\delta>1/2$, which we now consider. The bias term is equal to

C.3 Optimal regularization parameters

We can now take all six cases, and compute the optimal $\lambda$ and the resulting optimal regularization error.

$\mu_{i}=i^{-2\beta}$, $\nu_{i}=i^{-2\delta}$ ($2\delta>4\beta+1$): we need to minimize with respect to $\lambda$ the function $n^{-1}\lambda^{-1/2\beta}+\lambda^{2}$, which leads to $\lambda\approx n^{-1/(2+1/2\beta)}$ and an optimal value of $n^{1/(4\beta+1)-1}$.

$\mu_{i}=i^{-2\beta}$, $\nu_{i}=i^{-2\delta}$ ($2\delta<4\beta+1$): we need to minimize with respect to $\lambda$ the function $n^{-1}\lambda^{-1/2\beta}+\lambda^{(2\delta-1)/2\beta}$, which leads to $\lambda\approx n^{-\beta/\delta}$ and an optimal value of $n^{1/(2\delta)-1}$.

$\mu_{i}=i^{-2\beta}$, $\nu_{i}=e^{-\kappa i}$: same computation as the first one.

$\mu_{i}=e^{-\rho i}$, $\nu_{i}=i^{-2\delta}$: we need to minimize with respect to $\lambda$ the function $n^{-1}\log\frac{1}{\lambda}+(\log\frac{1}{\lambda})^{1-2\delta}$, which leads to $\log\frac{1}{\lambda}\approx n^{1/2\delta}$ and an optimal value of $n^{1/2\delta-1}$.

$\mu_{i}=e^{-\rho i}$, $\nu_{i}=e^{-\kappa i}$ ($\kappa>2\rho$): we need to minimize with respect to $\lambda$ the function $n^{-1}\log\frac{1}{\lambda}+\lambda^{2}$, which leads to $\lambda\approx n^{-1/2}$ and an optimal value of $\log n/n$.

$\mu_{i}=e^{-\rho i}$, $\nu_{i}=e^{-\kappa i}$ ($\kappa<2\rho$): we need to minimize with respect to $\lambda$ the function $n^{-1}\log\frac{1}{\lambda}+\lambda^{\kappa/\rho}$, which leads to $\lambda\approx n^{-\rho/\kappa}$ and an optimal value of $\log n/n$.

Appendix D Kernels on [0,1]

In this appendix, we consider kernels on $\mathcal{X}=$ that lead to closed-form expressions (or asymptotic equivalents) for eigenvalues of $K$ and components of $z$. These are used in simulations.

For a positive summable sequence $(\mu_{i})_{i\geqslant 1}$, we consider $k(x,y)=\sum_{i=1}^{\infty}2\mu_{i}\cos 2i\pi(x-y)$. It is defined for any $(x,y)\in^{2}$ and is 1-periodic in $x$ and $y$. It is a function $g$ of $x-y-\lfloor x-y\rfloor$, i.e., $k(x,y)=g(x-y-\lfloor x-y\rfloor)$. Moreover $k(x,x)$ is independent of $x$.

For $\mu_{i}=\frac{1}{i^{2\beta}}$, we have $k(x,y)=\frac{1}{(2\beta)!}B_{2\beta}(x-y-\lfloor x-y\rfloor)$, where $B_{2\beta}$ is the $(2\beta)$-th Bernoulli polynomial (Wahba, 1990). For example, we have $B_{2}(x)=\textstyle x^{2}-x+\frac{1}{6}$ and $B_{6}(x)=\textstyle x^{6}-3x^{5}+\frac{5}{2}x^{4}-\frac{1}{2}x^{2}+\frac{1}{42}$.

For $\mu_{i}=e^{-\rho i}$, we have, $k(x,y)=2\frac{e^{\rho}\cos 2\pi(x-y)-1}{e^{2\rho}-2e^{\rho}\cos 2\pi(x-y)+1}$. Indeed, we have

Data and eigenvectors.

If we consider $n$ data points $x_{i}=\frac{i-1}{n}$, $i=1,\dots,n$, then the kernel matrix $K$ has components $K_{ij}=k(\frac{i-1}{n},\frac{j-1}{n})$. It is a circulant matrix, thus it is diagonalizable in the discrete Fourier basis (Gray, 2006), with eigenvalues equal to the discrete Fourier transform of the first column of the matrix, i.e., $(g(0),g(1/n),\dots,g(1-1/n))^{\top}$.

Thus, the $i$-th eigenvector has $j$-th component $\frac{1}{\sqrt{n}}e^{2\omega i(j-1)\pi/n}$ (with $\omega^{2}=-1$) and the $i$-th eigenvalue is

If $n$ is large and $\mu_{i}$ tends to zero when $i$ tends to $+\infty$, then an asymptotic equivalent for $\lambda_{i}$ is $n\mu_{i}$.

For data sampled from the uniform distribution in $$, then similar equivalents hold (see, e.g., Harchaoui et al., 2008).

Functions.

Let $f(x)=\sum_{i=1}^{\infty}2\nu_{i}^{1/2}\cos 2i\pi x$, for $\nu_{i}$ a non-negative summable sequence. We consider $z_{i}=f(x_{i})=f((i-1)/n)$. The component of $z$ on the $i$-th eigenvector of $K$ is (following the same reasoning as above):

and the asymptotic equivalent is $(n\nu_{i})^{1/2}$.

Link with Sobolev spaces.

The kernel $k(x,y)$ defined above corresponds for $\mu_{i}=i^{-2\beta}$ to certain Sobolev spaces (Wahba, 1990; Gu, 2002). Indeed, for $\beta$ integer, the associated RKHS is the Sobolev space of periodic functions which are $\beta$-times differentiable.

Moreover, when $\nu_{i}=i^{-2\delta}$, then for $\delta>\delta_{0}$, then the corresponding function is $(\delta_{0}-1/2)$-times differentiable, and the minimax rate of estimation is known to be exactly $O(n^{1/2\delta_{0}})$ (Speckman, 1985; Johnstone, 1994). Thus, up to logarithmic terms, the best possible rate is $O(n^{1/2\delta})$, and is achieved if $\beta$ is large enough (Section 4.3).

References