PaperlayerView exact revision on arXiv

Source-linked AI summary

Sharp analysis of low-rank kernel matrix approximations

Francis Bach

arXiv:1208.2015v3cs.LGmath.ST

TL;DR

Kernel methods often require at least quadratic computation because they use an n-by-n kernel matrix, motivating low-rank approximations. The paper analyzes random column sampling for kernel ridge regression and shows that the rank can scale with degrees of freedom while preserving predictive performance with subquadratic algorithms.

  • Problem

    Kernel methods are computationally limited by kernel-matrix construction, and the rank needed for low-rank approximations remains the key practicality question.

  • Method

    The paper analyzes kernel ridge regression using low-rank approximations formed from a random subset of kernel-matrix columns.

  • Results

    The rank p can be chosen linear in the problem's degrees of freedom, yielding subquadratic algorithms with the same predictive performance as classical quadratic-time methods for square loss.

  • Takeaways & Limitations

    Degrees of freedom characterize both predictive performance and the computational rank required by the analyzed column-sampling approximation.

  • Takeaways & Limitations

    The analysis focuses on kernel least-squares regression and fixed-design settings, while greedy sampling has potentially better behavior but remains harder to analyze.

Abstract

from arXiv · show

We consider supervised learning problems within the positive-definite kernel framework, such as kernel ridge regression, kernel logistic regression or the support vector machine. With kernels leading to infinite-dimensional feature spaces, a common practical limiting difficulty is the necessity of computing the kernel matrix, which most frequently leads to algorithms with running time at least quadratic in the number of observations n, i.e., O(n^2). Low-rank approximations of the kernel matrix are often considered as they allow the reduction of running time complexities to O(p^2 n), where p is the rank of the approximation. The practicality of such methods thus depends on the required rank p. In this paper, we show that in the context of kernel ridge regression, for approximations based on a random subset of columns of the original kernel matrix, the rank p may 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, and is often seen as the implicit number of parameters of non-parametric estimators. This result enables simple algorithms that have sub-quadratic running time complexity, but provably exhibit the same predictive performance than existing algorithms, for any given problem instance, and not only for worst-case situations.

1. Introduction

The paper asks whether positive-definite kernel methods can achieve subquadratic computation without losing predictive performance. It studies random column approximations and relates the required rank to problem-specific degrees of freedom.

  • Kernel methods require at least quadratic time in n because computing the kernel matrix has that complexity.
  • The paper frames degrees of freedom as an implicit parameter count that supports problem-dependent rather than only worst-case analysis.
  • Random subsets of p kernel columns yield approximations with running time O(p^2 n) and computable error bounds.
  • In fixed-design least-squares regression, p can be chosen linear in the problem's degrees of freedom for any problem instance.
  • The resulting algorithms have subquadratic running time while matching the predictive performance of classical quadratic-time algorithms for square loss.
  • For rapidly decaying spectra, computational limits may prevent exploring relevant regularization paths and lead to underfitting.

2. Supervised learning with positive-definite kernels

This section places kernel learning in equivalent primal and dual formulations and reviews computational approaches to approximating or solving kernel problems. It emphasizes that existing approximation-error guarantees may not explain observed predictive performance.

  • 2.1. Equivalent formulations: Kernel supervised-learning objectives can be formulated in an RKHS and transformed into optimization problems over R^n.
  • 2.1. Equivalent formulations: The representer theorem expresses the solution as a linear combination of feature maps evaluated at the training points.
  • 2.1. Equivalent formulations: The kernel matrix K has entries K_ij = k(x_i, x_j), providing the finite-dimensional representation used in these formulations.
  • 2.1. Equivalent formulations: For convex losses, a dual formulation is available; for square and logistic loss, its solution is unique and also solves the primal formulation.
  • 2.2. Related work: Primal optimization is often inefficient because kernel problems are ill-conditioned, while dual methods generally have better condition numbers but still require the full kernel matrix.
  • 2.2. Related work: Column-sampling methods reduce approximation costs, but conventional norm-error bounds can be too pessimistic to match practical predictive performance.

3. Approximation from subset of columns

The section constructs a low-rank kernel approximation from randomly sampled columns, relates it to an explicit p-dimensional feature map, and contrasts random sampling with pivoted greedy selection.

  • Approximation from columns: A random subset I of p indices defines an approximation L using the columns K(V, I) of the kernel matrix.The construction uses the kernel matrix columns indexed by I as its available information.
  • Approximation from columns: L is the unique symmetric matrix whose column space is spanned by K(V, I) and whose sampled columns equal K(V, I).It can also be interpreted through orthogonal projections of the data points onto the subspace spanned by the sampled points.
  • Approximation from columns: The approximation creates an explicit feature map in R^p, enabling predictions on test points and improved test-time running performance.The feature map is constructed from K(I, I)^−1/2 and kernel evaluations against the sampled indices.
  • Approximation from columns: O(p^2n) running time suffices to obtain the feature map using incomplete Cholesky decomposition, with an online trace-norm error bound.Incomplete Cholesky is often interpreted as partial Gram–Schmidt orthonormalization.
  • Pivoting vs. random sampling: Pivoted incomplete Cholesky has the same complexity as random subsampling and may achieve better approximation performance through approximate greedy selection.The paper compares the two approaches but notes that analyzing the greedy algorithms remains an open problem.

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

Under fixed design and square loss, the analysis characterizes ridge-regression prediction error and shows that random column sampling can preserve it with a rank tied to degrees of freedom. The guarantees are instance-based and expectation-based, while the analysis identifies scope limits beyond square loss and fixed design.

  • The analysis assumes deterministic inputs, real-valued responses, square loss, and noise with zero mean and finite covariance, without requiring independent or identically distributed noise.
  • The ridge predictor is a smoothed estimate, obtained as ˆz = K(K + nλI)−1y, and its expected prediction error is analyzed through a bias/variance decomposition.
  • Degrees of freedom quantify the variance contribution like an implicit parameter count, and the required number of sampled columns is linear in this quantity up to logarithmic terms.
  • Theorem 1 gives instance-based guarantees: for given λ and K, column sampling can avoid average predictive-performance loss, unlike cited worst-case guarantees.
  • The relative bound can make sampled-kernel predictions perform as well as full-kernel predictions in expectation, but high-probability guarantees over the noise require additional assumptions.
  • The framework extends conceptually to smooth losses and random designs, but logistic-loss analysis becomes more complex and random designs require controlling empirical covariance deviations.

5. Simulations

The simulations examine how kernel spectrum decay, regularization, and degrees of freedom affect prediction and the rank required by low-rank approximations. Across datasets, sufficient ranks remain close to the degrees of freedom, supporting the theoretical result.

  • Synthetic examples: The experiments vary periodic smoothing-spline kernels and functions to study regularization parameters and degrees of freedom under different eigenvalue-decay behaviors.Synthetic problems use kernels with prescribed spectral sequences, including µ_i = i^-2β and ν_i = i^-2δ.
  • Optimal values of λ: For β = 1, the convergence rate n^(1/(4β+1)−1) is achieved but is slower than n^(1/(2δ)−1).
  • Optimal values of λ: For β = 4, the optimal convergence rate n^(1/(2δ)−1) is achieved as expected.
  • Performance of low-rank approximations: The sufficient rank is defined as the rank yielding predictive performance at most 1% worse than the full-rank case, normalized by degrees of freedom.Experiments use both marginal and average degrees-of-freedom definitions.
  • Performance of low-rank approximations: Across all datasets, sufficient-rank ratios stay relatively close to one, while pivoting generally changes results little and can reduce required columns by a constant factor.Sudden increases in some Figure 3 plots are attributed to instability in the 1%-worse-performance criterion.

6. Conclusion

The conclusion identifies degrees of freedom as both a statistical and computational guide for column-sampling approximations in kernel least-squares regression. It also outlines extensions to sampling schemes, designs, losses, and computational settings.

  • Conclusion: Column sampling analysis shows that the approximation rank may be chosen proportional to properly defined degrees of freedom.The result holds for any problem instance, rather than only in a worst-case regime.
  • Conclusion: Computing degrees of freedom exactly requires the full kernel matrix, although efficient approximation methods may be used instead.
  • Conclusion: Figure 2 compares RKHS size, predictive performance, and sufficient-rank ratios for random sampling and pivoted incomplete Cholesky.
  • Conclusion: Figure 3 reports sufficient-rank ratios for pumadyn datasets 32fh, 32nh, and 32nm under the same two column-selection methods.
  • Conclusion: Future work includes nonuniform column sampling, random-design analysis, other losses, and extensions beyond the batch setting.

Appendix B. Proof of Theorem 1

The proof of Theorem 1 begins by establishing a Bernstein-type concentration inequality for subsampled covariance matrices.

  • Appendix B. Proof of Theorem 1: The lemma follows prior matrix-concentration work and supplies the probabilistic tool used in the theorem’s proof.

B.1. Concentration of subsampled covariance matrices

This subsection proves concentration for a covariance matrix formed from a uniformly sampled subset of rows. The argument converts sampling without replacement into a form suitable for matrix Bernstein bounds.

  • B.1. Concentration of subsampled covariance matrices: For a matrix Ψ with bounded row norms, a uniformly sampled subset I of p rows defines the subsampled matrix Ψ_I.
  • B.1. Concentration of subsampled covariance matrices: The proof introduces a centered matrix difference Δ and establishes that its expectation is zero.
  • B.1. Concentration of subsampled covariance matrices: An auxiliary matrix Ξ is constructed by independently sampling p rows with replacement to support the concentration argument.
  • B.1. Concentration of subsampled covariance matrices: Applying matrix Bernstein’s inequality yields the required probability bound and completes the concentration result.

B.2. Proof of Theorem 1

The proof controls the approximation’s bias and variance through matrix inequalities and probabilistic bounds for a random subset of columns. Monotonicity and regularization arguments extend the result to the unregularized case, while the analysis acknowledges variance-bound looseness and possible predictive degradation.

  • Proof setup: The proof factorizes K as ΦΦ⊤ and analyzes the approximation through the matrix Nγ.For γ > 0, the argument later obtains the γ = 0 case by monotonicity.
  • Variance control: The variance term for Nγ is bounded by comparing Nγ with N0 and I through matrix monotonicity.Because 0 ≼ Nγ ≼ N0 ≼ I, the approximated variance term is no larger than the corresponding unapproximated comparison term.
  • Bias control: The bias term is controlled by lower-bounding Nγ, which reduces the problem to an upper bound on an associated matrix expression.The resulting bound compares the approximated bias with the original bias through a γ- and λ-dependent factor.
  • Probabilistic control: The probabilistic step bounds the largest eigenvalue of the discrepancy between the full empirical moment and a random subset moment.The argument assumes γ ⩽ λ and bounds the squared ℓ2-norms of rows of Ψ before applying the concentration lemma.
  • Caveats: The bound could be sharpened by explicitly reducing the variance term, and approximation performance may sometimes be smaller than the non-approximated version.These statements qualify the sharpness and predictive-performance scope of the proof.
  • Probabilistic control: The proof relies on a high-probability bound over the sampled columns rather than a deterministic guarantee.Its final result is presented after selecting t = 1/2 and γ = λδ.

C.2. Bias terms

This section derives bias bounds under polynomial, exponential, and mixed eigenvalue–signal decay patterns. The bounds depend jointly on kernel eigenvalues and signal components, with polynomial cases splitting according to the relative decay exponents.

  • Decay cases: The bias term depends on both the kernel eigenvalue sequence (µi) and the signal-component sequence (νi).The analysis therefore considers combinations of polynomial and exponential decay.
  • Polynomial decays: If 2δ − 4β > 1, the polynomial-decay bias is O(nλ2).The bound follows because the relevant integral is finite.
  • Polynomial decays: When 2δ − 4β < 1, the polynomial-decay bias is further bounded using a finite integral under the assumptions on β and δ.The supplied passage does not state the resulting order explicitly.
  • Exponential decays: For exponential decays of both sequences, the bias bound uses an integral that is finite and uniformly bounded in λ.This gives a λ-uniform control of the relevant bias expression.
  • Mixed decays: With polynomial kernel decay and exponential signal decay, the signal decays fast enough that the bias is nλ2.The argument transfers the earlier sufficiently-fast-signal result to this mixed case.

C.3. Optimal regularization parameters

The section optimizes λ by balancing the variance and bias orders for different kernel–signal decay regimes. The resulting rates vary with the decay exponents, with logarithmic factors appearing in exponential cases.

  • Optimization: The analysis computes optimal λ and regularization error across six combinations of kernel and signal eigenvalue decays.Each case minimizes the corresponding bias–variance upper bound.
  • Polynomial decays: λ ≈ n^−1/(2+1/(2β)) and the optimal value is n^(1/(4β+1))−1 when 2δ > 4β + 1.This is the polynomial-decay case with µi = i^−2β and νi = i^−2δ.
  • Mixed decays: For µi = i^−2β and νi = e^−κi, the optimal regularization calculation is the same as in the first exponential case.The section explicitly identifies this mixed case by reference to that computation.
  • Exponential decays: λ ≈ n^−1/2 and the optimal value is log n/n when µi = e^−ρi and νi = e^−κi with κ > 2ρ.The minimized objective is λ + λ2.
  • Exponential decays: λ ≈ n^−ρ/κ and the optimal value is log n/n when µi = e^−ρi, νi = e^−κi with κ < 2ρ.The minimized objective is λ + λ^κ/ρ.

Appendix D. Kernels on [0,1]

The appendix constructs periodic kernels on [0,1] with tractable spectral behavior and uses their eigenvalues and signal components in simulations. It also connects polynomial decay to Sobolev spaces and known estimation rates.

  • Kernels: The appendix studies kernels on X = [0,1] whose eigenvalues and signal components have closed-form expressions or asymptotic equivalents for simulations.The construction uses positive summable eigenvalue sequences.
  • Kernels: The kernel is 1-periodic, depends on x − y modulo one, and has k(x,x) independent of x.These properties follow from its Fourier-series construction.
  • Kernels: For µi = i^−2β, the kernel can be expressed using a Bernoulli polynomial.The appendix gives the corresponding periodic-kernel form and an example for B2(x).
  • Kernels: For µi = e^−ρi, the appendix gives an explicit trigonometric kernel expression.This provides a closed-form example for exponential eigenvalue decay.
  • Spectral quantities: The kernel matrix for equally spaced points is circulant and diagonalizable in the discrete Fourier basis.Its eigenvalues are the discrete Fourier transform of the first column.
  • Spectral quantities: For large n with µi tending to zero, the kernel-matrix eigenvalues have asymptotic equivalent nµi.This relation is used to connect the sequence-level decay assumptions to finite-sample matrices.
  • Signal components: The signal component along the i-th eigenvector has asymptotic equivalent (nνi)1/2.The appendix derives this for the specified Fourier-type signal construction.
  • Sobolev connection: For polynomial νi, the associated function has controlled differentiability and the minimax estimation rate is achieved up to logarithmic terms when β is sufficiently large.The passage states the best possible rate as O(n1/2δ).
Loading 1208.2015v3…